Int J Fract DOI 10.1007/s10704-015-0051-0

ORIGINAL PAPER

Modeling the interaction between fluid-driven fracture and natural fault using an enriched-FEM technique A. R. Khoei · M. Vahab · M. Hirmand

Received: 25 March 2015 / Accepted: 3 November 2015 © Springer Science+Business Media Dordrecht 2015

Abstract In this paper, the interaction between the fluid-driven fracture and frictional natural fault is modeled using an enriched-FEM technique based on the partition of unity method. The intersection between two discontinuities is modeled by introducing a junction enrichment function. In order to model the fluid effect within the fracture, the fluid pressure is assumed to be constant throughout the propagation process. The frictional contact behavior along the fault faces is modeled using an X-FEM penalty method within the context of the plasticity theory of friction. Finally, several numerical examples are solved to illustrate the accuracy and robustness of proposed computational algorithm as well as to investigate the mechanism of interaction between the fluid-driven fracture and the natural fault. Keywords Multiple cracks · X-FEM method · Fluid-driven fracture · Frictional natural faults

1 Introduction Hydraulic fracturing is known as an appealing technique in oil industry for stimulation of gas and oil reservoirs. In a hydraulic fracturing operation, a highly A. R. Khoei (B) · M. Vahab · M. Hirmand Center of Excellence in Structures and Earthquake Engineering, Department of Civil Engineering, Sharif University of Technology, P.O. Box 11365-9313, Tehran, Iran e-mail: [email protected]

pressurized fluid is pumped into the borehole to create fractures that form a path through which the oil stored in the reservoir flows toward the borehole (Dong and De Pater 2001). The behavior of hydraulic fracture near a natural fault or discontinuity is of great importance for an efficient reservoir simulation, as natural discontinuity can significantly influence the hydraulic fracturing process (Zhang and Ghassemi 2011). However, the primary models proposed for hydraulic fracturing process were performed based on a neat crack condition within the domain (Réthoré et al. 2007; Mohammadnejad and Khoei 2013; Khoei et al. 2014; Khoei and Vahab 2014), which is proven to be over idealized as nearly all reservoirs contain natural discontinuities substantially affecting the hydraulic fracturing process (Zhang and Ghassemi 2011; Dahi and Olson 2011; Khoei et al. 2015). Dahm (2000) demonstrated that the implementation of a single-crack model in the simulation of hydraulic fracturing process results in considerable loss of accuracy. In hydraulic fracturing process, the natural fault may experience the mixed mode opening/closing conditions, which in turn affects the hydraulic fracture behavior when the hydraulic fracture merges with the natural fault (Zhang et al. 2007), as shown in Fig. 1a. In the case that a junction takes place, different scenarios are conceivable for the hydro-fracture/natural-fault system, including the penetration, arrest, diversion and offset, as shown in Fig. 1b. It has been shown through a wide range of laboratory and field observations that hydro-fracture branch-

123

A. R. Khoei et al.

(b) Natural fault

Natural fault

Hydraulic fracture

Hydraulic fracture

(a)

σv

Region of contact

Penetration

σh

σh

Hydraulic fracture

Natural fault

Arrest Natural fault

Natural fault Hydraulic fracture

Hydraulic fracture

σv Diversion

Offset

Fig. 1 a The opened and closed zones along the natural fault, b schematic view of possible scenarios for a hydraulic fracture merged with a natural fault

ing, arrest or diversion is not only possible but also fairly common. The changes in rock properties and insitu stresses substantially affect the hydraulic fracture growth occurring in layered sedimentary rocks (Price and Cosgrove 1990). Several investigations on the complex behavior of hydraulic fracturing subjected to geologic discontinuities such as faults, joints and bedding planes have been documented from mineback data and laboratory tests by Warpinski and Teufel (1987) and Jeffrey and Settari (1995). Daneshy (1974) performed a set of experiments to study the hydraulic fracturing in the presence of natural flaws; it was observed that the hydro-fractures appear to cross the natural faults when they are closed; otherwise the hydro-fracture would be arrested by the natural fault. Anderson (1981) carried out experimental studies on hydraulically-driven fractures in the vicinity of unbounded interfaces in rocks. Based on experimental tests performed on Nugget sandstone from Utah and Indiana Limestone under uniaxial loading condition, it was shown that the hydrofracture growth is arrested below a threshold of confining normal stress which is proportional to the friction between the natural fracture surfaces. Blanton (1982) extended the works of Anderson (1981) by

123

performing laboratory experiments on naturally fractured blocks of Devonian shale and hydrostone under triaxial state of stress. It was shown that the hydrofractures were mostly either diverted or arrested by the natural faults unless high differential stresses and high angles of approach are provided. Murphy and Fehler (1986) claimed that the shear slippage along the natural discontinuities can be activated much faster than the conventional tensile failure, especially in the presence of high differences between the minimum and maximum in-situ stresses. Based on their observations, the occurrence of slippage along the natural fault faces leads to the hydro-fracture branching, or dendritic evolution patterns, which are in agreement with microearthquake locations. The numerical study of the interaction between hydraulically driven fracture and pre-existing discontinuities has been the issue of interest by many researchers. As a matter of fact, since a long list of parameters, such as the fluid viscosity, orientation of discontinuities, crack lengths, frictional resistance and farfield stresses, are determinant in the overall behavior of interacting discontinuities, a comprehensive study of the interaction problem is quite cumbersome and

Modeling the interaction between fluid-driven fracture and natural fault

complicated tasks. Basically, the research works studied on the interaction of hydraulic fracturing and preexisting discontinuities can be categorized in two main groups. In the first category, the study is mainly focused on the interaction between hydraulic fracture and natural faults; Dong and De Pater (2001) proposed an algorithm for simulation of hydraulic fracturing and its interaction with faults that was based on the displacement discontinuity method. Akulich and Zvyagin (2008) performed a numerical simulation to study the interaction between the fracture and natural fault where the fracture is propagated by an incompressible Newtonian fluid injected from a source located at its center. It was shown that the lower value of fluid viscosity results in a rapid hydraulic fracture propagation; it was also demonstrated that for higher fluid viscosity the interaction begins from a greater distance between the hydraulic fracture and natural fault. Zhou et al. (2008) studied the hydraulic fracture propagation behavior in naturally fractured reservoirs through a series of servo-controlled tri-axial fracturing experiments, and observed three types of interaction between the induced fracture and pre-existing fracture by operating different values of horizontal stress, angle of approach, and shear strength of pre-existing fracture. Dyskin and Caballero (2009) investigated the interaction between the hydraulically driven fracture and frictionless natural fault using the finite element method, and illustrated that a relatively long frictionless and cohesionless fault is capable of arresting the hydraulic fracture propagation. Dahi and Olson (2011) presented a numerical model to study fracture intersections by tracking fluid fronts in the network of reactivated fissures, where the hydraulic fracture was arrested by pre-existing natural fractures, and/or was controlled by shear strength and potential slippage at the fracture intersections. Zhang and Ghassemi (2011) performed a comprehensive study on the interaction between the hydraulic fracture and natural fault, and concluded that the fault influence is conditioned by its shear stiffness, its inclination, and its distance from the hydraulic fracture; it was also highlighted that the hydraulic fracture always tends to propagate along the maximum compressive stress direction. On the other hand, the second category concerns with the interaction between the hydraulically-driven fracture and material interfaces, such as layered media or bedding planes. Li (2000) investigated the interface debonding when a crack perpendicularly approaches an

interface between two dissimilar elastic materials, and emphasized that the most significant factor in the interface debonding is the interface toughness while the relative bimaterial stiffness is less important. Gudmundsson and Brenner (2001) and Brenner and Gudmundsson (2004) studied the hydro-fracture arrest and aperture variation in layered reservoirs, and highlighted that the growth of a hydro-fracture depends on the mechanical properties of the host rock and the over-pressure of the hydro-fracture; in this regard, it was noted that if the internal fluid pressure is the only loading condition, by approaching the hydro-fracture to a stiff layer, the tip becomes sharp and narrow, and continue its propagation through that layer since the high tensile stresses concentrate at the stiff layers, however—the soft layers suppress the tensile stresses associated with the hydrofracture tip and blunt the tip itself, so the hydro-fracture tends to become arrested on the intersection with a soft layer. Zhou et al. (2008) performed a numerical fracture analysis to study the propagation of hydraulic fracturing deflected at bedding interfaces in layered rocks using a two-dimensional boundary element model, and concluded that the fracture escape from the interface is likely under the conditions of fracture growth from stiff to soft rocks, small layer-to-layer far-field stress contrasts, and moderately low fluid viscosity, and small parent fracture lengths and offset distances; as a result, the change of the fracture propagation direction at the bedding plane leads to a variety of fracture and fluid flow patterns. It has been extensively shown that modeling of crack propagation with the classical finite element method is quite a cumbersome task since the evolving nature of the problem requires one to update the domain topology to conform the FE mesh to the faces of the crack as it propagates. Thus, the extended finite element method has been developed to eliminate the requirement of remeshing by incorporating additional discontinuous functions into the standard finite element approximation on the basis of the partition of unity method (Melenk and Babuška 1996). In this method, the discontinuity is modeled by introducing enrichment degrees of freedom for the nodal points which their supports are cut by the discontinuity (Belytschko and Black 1999; Moës et al. 1999; Stolarska et al. 2001). Hence, the geometry of the discontinuity can be updated independent of the FE mesh simply by updating the enrichment scheme. Moreover, by introducing the asymptotic crack-tip enrichment functions, the con-

123

A. R. Khoei et al.

vergence of the solution can be achieved even with a uniform course mesh. The X-FEM introduces a great simplicity while dealing with multiple crack growth problems where the probable junction and/or intersection between the discontinuities can be simply modeled by introducing appropriate junction enrichment functions. The X-FEM was performed by Daux et al. (2000) to model the intersecting cracks, in which a technique was proposed to construct the enriched approximation space for multiple discontinuities. The implementation of X-FEM in modeling of multiple cracks was also carried out by Chopp and Sukumar (2003) and Zi et al. (2004) for fatigue crack growth, and by Budyn et al. (2004) for brittle crack growth. The X-FEM has also been well developed to model frictional discontinuities on the basis of both penalty and Lagrangian approaches (Khoei and Nikbakht 2007; Liu and Borja 2008; Hirmand et al. 2015). A comprehensive overview of the X-FEM method and its applications in continuum mechanics were presented by Khoei (2015). In the present study, the interaction between the fluid-driven fracture and pre-existing natural fault is simulated on the basis of an extended finite element method, where a junction enrichment function is introduced to model the intersection of hydraulic fracture with the natural fault. The effect of hydraulic fracturing fluid within the propagating crack is modeled by imposing a uniform constant pressure along the crack faces throughout the propagation process. During the hydraulic fracturing process, the natural fault may experience the mixed mode opening/closing conditions. In order to model the frictional contact behavior on the faces of natural fault, an X-FEM penalty method is employed within the context of the plasticity theory of friction. Several numerical examples are analyzed to verify the robustness and versatility of the computational algorithm, and to investigate the interaction between the hydraulic fracture and natural fault.

t

n t

n

fault

θ Q

HF

n

x2 HF

x1

fault

u

u Fig. 2 Problem configuration and the boundary conditions

here is restricted to infinitesimal deformation regime with the linear elastic behavior for the bulk system. Consider a body  bounded by the external boundary , as shown in Fig. 2, in which the boundary  is subjected to the prescribed traction t on t and the displacement boundary condition u¯ on u , such that t ∪ u =  and t ∩ u = ∅. Moreover, the domain contains the internal discontinuities denoted by d that generally consist of two parts, one denoted by HF represents the hydraulically-driven fracture which is subjected to the fracturing fluid pressure pHF , and the other denoted by fault represents the frictional natural fault that can experience the frictional contact condition and is subjected to the contact traction tcont . The equilibrium equation of the fractured body can then be written using the well-known momentum balance equation expressed in the quasi-static condition as ∇ · σ + ρb = 0

(1)

2 Governing equations of a fractured domain

where b is the body force per unit mass vector, ρ is the density of the medium, and σ is the Cauchy stress tensor that can be related to the strain tensor ε using an appropriate constitutive relation. The boundary conditions associated with the external and internal boundaries of the domain are represented as σn = t on t u = u˜ on u (2) σnd = p on d

In order to model the crack growth simulation in a naturally fractured domain, the governing equations of a cracked body are first derived; the weak form of governing equations together with the X-FEM discretization equations are then obtained to investigate the interaction between the hydraulic fracture and frictional natural fault. It must be noted that the discussion given

where t is the prescribed traction on t , u˜ is the prescribed displacement on u , p is the traction imposed on the faces of the discontinuities d and, n and nd are the unit normal vectors on the external and internal boundaries, respectively. Note that the traction p imposed on HF and fault are respectively denoted by pHF and tcont defined as

123

Modeling the interaction between fluid-driven fracture and natural fault

σnHF = pHF on HF ⊂ d σnfault = tcont on fault ⊂ d

(3)

The constitutive equation relating the stress tensor to the strain tensor can be expressed by σ = Dε, where D is the elasticity material property matrix of the continuum body, and ε is the infinitesimal strain tensor defined by ε(u) = ∇ s u, with ∇ s denoting the symmetric spatial gradient operator, i.e. ∇ s = 21 (∇ + ∇ T ). Basically, the distribution of fracturing fluid pressure along the hydraulic fracture faces must be determined through the solution of continuity equation of the crack inflow defined based on the Darcy law, coupled with the above mentioned boundary value problem to account for the effect of fluid viscosity on the fluid pressure distribution. It is important to note that the effect of fracturing fluid viscosity may result in formation of the so called pressure-lag at the vicinity of hydraulic fracture tip (Réthoré et al. 2007; Mohammadnejad and Khoei 2013). However in the current study, the effect of fluid viscosity is neglected for the sake of simplicity, and the crack inflow is simply modeled by imposing a constant uniform pressure along the faces of hydraulically-driven fracture HF . Since the region of pressure-lag is quite small compared to the fracture length, the assumption of constant pressure along the hydro-fracture faces induces negligible error in the model (Boone and Ingraffea 1990), particularly if a low value of fluid viscosity is imposed. Moreover, the constant pressure approach introduces considerable simplicity in the numerical simulation, and has been extensively used by various researchers (Dong and De Pater 2001; Dyskin and Caballero 2009; Zhang and Ghassemi 2011). Consider the crack inflow pressure is denoted by p, the traction imposed by the fluid flow along the faces of hydraulic fracture HF can be expressed as pHF = − pnHF . In order to model the frictional contact behavior along the natural fault fault , the so called Kuhn–Tucker inequalities are imposed to account the boundary value problem for the contact condition at fault , as demonstrated in the next section. Moreover, to derive the weak form of mentioned equilibrium boundary value problem, i.e. Eqs. (1)–(3), the trial and test functions of the displacement field are introduced. The trial function of the displacement field u(x, t) is assumed to satisfy all essential boundary conditions on u while being smooth enough to define derivatives of the formulation. Moreover, the test func-

tion δu(x, t) is required to have the properties of trial function except to be zero on the boundaries where the essential boundary conditions are imposed. Following the standard procedure of derivation of the weak form of the equilibrium equation, the weak form of the governing boundary value problem can be obtained as 

  δε:σ d + δu · pHF d + δu · t cont d  HF fault   δu · ρb d + δu · t d (4) = 



which must be satisfied for all arbitrary test functions δu. Note that the notation δu stands for the jump value of δu across the discontinuities HF and fault . It is noteworthy to highlight that the second and third integral terms are derived to define the tractions imposed on the faces of HF and fault due to the fluid injection pressure along the hydraulic fracture and the contact constraint condition at the natural fault, respectively.

3 Modeling frictional contact along the natural fault In order to model the frictional contact condition along the natural fault faces, the frictional contact algorithm is incorporated into the X-FEM on the basis of penalty method within the context of the plasticity theory of friction (Khoei and Nikbakht 2007; Liu and Borja 2008; Khoei et al. 2009). In the X-FEM penalty method, the hypothetical springs are embedded along the contact interfaces such that the contact traction tcont is related to the relative displacement of contacting boundaries without introducing any new unknowns. The relative displacement of master and slave points located on the faces of natural fault fault can be evaluated using the normal and tangential gap functions, as shown in Fig. 3a, by   (5) on fault g N = nfault ⊗ nfault u   gT = I − nfault ⊗ nfault u on fault (6)   in which nfault ⊗ nfault is the normal projection tensor. The unilateral contact condition is imposed on the surface of natural fault fault through the standard Kuhn-Tucker conditions as g N ≥ 0, t Ncont ≤ 0, g N t Ncont = 0

(7)

123

A. R. Khoei et al. Fig. 3 a Definition of the normal and tangential gap functions, b the normal and tangential contact tractions

(a)

(b) Slave

t

n

cont

t cont N

S u

gT

gN

M

Master

in which t Ncont = tcont N · nfault denotes the normal component of the contact traction, and g N = u · nfault is the normal opening of the interface, as shown in Fig. 3b. In order to resolve the unilateral contact condition, the is related to the penetranormal contact traction tcont N tion of contacting boundaries using a normal penalty constant k N as tcont N = k N g N . It must be noted that the parameter k N is analogous to the normal stiffness of the interface. In the tangential direction, two conditions of stick and slip are distinguished. Considering the stick condition, the procedure is similar to that employed for the normal contact traction, i.e. the relative tangential movement along the interface is constrained by the hypothetical tangential springs with a tangential penalty parameter k T as tTcont = k T gT . However, the perfect stick condition is quite exceptional in contact problem, and the plasticity theory of friction is used to describe the stick-slip motion of contacting boundaries in the tangential direction. To this aim, a slip criterion is defined based on the Coulomb friction law as   cont   cont  = 0 Slip     − μ f tN (8) F f = tT < 0 Stick in which μ f is the Coulomb friction coefficient. Following this, a constitutive relation that accounts for both the stick and slip conditions can be achieved for the frictional contact tractions as   (9) dtTcont = k T dgT − dλmfault where dλ = |dgT | is the magnitude of sliding increment and mfault is the unit vector in the tangential (sliding)  direction  of the natural fault, defined as mfault = tTcont tTcont . A criterion for the tangential stick/slip condition can be defined based on an alternative definition of Kuhn–Tucker relations in the tangential direc-

123

fault

tTcont

Slave

fault

m M

S

fault

fault

Master

tion by dλ > 0, F f < 0 and dλ · F f = 0, in which the stick condition pertains to the case of dλ = 0 and F f < 0, and the slip condition pertains to the case dλ > 0 and F f = 0 (Hirmand et al. 2015). A return mapping algorithm can be developed analogous to the predictor–corrector scheme of the classical theory of plasticity. The details of return mapping algorithm can be found in Liu and Borja (2008) and Khoei and Vahab (2014). Accordingly, a constitutive relation that relates the contact traction to the relative displacement along the interface can be defined as dtcont ≡ Dcont du, in which Dcont is the contact constitutive matrix defined for the stick and stick-slip conditions as Dcont      k N nfault ⊗ nfault  +k T mfault ⊗ mfault  for stick ≡ for stick-slip k N nfault ⊗ nfault +μ f k N mfault ⊗ nfault

(10) 4 X-FEM discretization of governing equations In order to derive the system of X-FEM equations, the extended finite element approximation is employed to discretize the weak form of governing equation (4). In the classical finite element framework, the field variables are approximated using the standard continuous shape functions. However in the X-FEM, additional degrees of freedom are introduced to enrich the standard continuous field with discontinuous functions in the context of partition of unity method. To account for the displacement jump across the natural fault and hydro-fracture faces, the displacement field is enriched by discontinuous functions along fault and HF . To this aim, the trial and test functions of the displacement field, i.e. u(x, t) and δu(x, t) respectively, are approximated using the extended finite element discretization, where the discontinuity interface is embed-

Modeling the interaction between fluid-driven fracture and natural fault

ded within the finite elements by incorporating the discontinuous interpolation functions into the standard finite element approximation. The resulting systems of nonlinear algebraic equations are then solved using the Newton–Raphson iterative procedure.

and

4.1 Approximation of the displacement field

where (r, θ ) is the polar coordinates with the origin located at the crack tip (Fig. 2). In relation (12), ϕ(x) is the signed distance function defined based on the absolute value of level set function; this is a common procedure in the extended finite element method to define the configuration of discontinuities. In this manner, a singed distance function is defined to attribute a distance and a sign to each point of the domain with respect to the interface as

In order to introduce an enriched displacement field for a domain with the hydro-fracture and natural fault, the Heaviside function is employed as an enrichment function to model the discontinuous displacement field induced by the hydro-fracture and natural fault; moreover, the crack tip asymptotic functions that are available from the analytical solutions of LEFM, are utilized to enrich the displacement field near the tips of hydro-fracture and natural fault. Hence, the extended finite element approximation of the displacement field uh (x, t) can be expressed as  Ni (x)u¯ i (t) u(x, t) ≈ uh (x, t) = +



i∈N std

Ni (x) (H(ϕ(x)) − H(ϕ(x i ))) a¯ i (t)

i∈N dis

+



⎛ Ni (x) ⎝

i∈N ti p

4  

B j (x) − B j (x i )



⎞ j b¯ i (t)⎠

j=1

(11a) or ¯ + Nadis (x)¯a(t) uh (x, t) = Nustd (x)u(t) ti p ¯ + N (x)b(t)

(11b) b where Ni (x) are the standard shape functions of node i, N std is the set of standard nodes, N dis is the set of nodes whose support is cut by the hydro-fracture and natural fault discontinuities, and N ti p is the set of nodes that contain the crack tip in their supports; u¯ i (t) are the j standard nodal DOFs and, a¯ i (t) and b¯ i (t) are additional nodal DOFs corresponding to the Heaviside and crack tip enrichment functions, respectively. In rela¯ is the vector of standard displacement tion (11b), u(t) ¯ DOFs and, a¯ (t) and b(t) are the vectors of enriched displacement DOFs associated with the Heaviside and crack tip enrichment functions, respectively. The crack tip asymptotic functions B j (x) and the Heaviside function H(ϕ(x)) are defined as  H (ϕ(x)) =

+1 0

ϕ(x) ≥ 0 ϕ(x) < 0

(12)

B(r, θ) = {B1 , B2 , B3 , B4 }   √ √ θ √ θ √ θ θ = r sin , r cos , r sin sin θ, r cos sin θ 2 2 2 2

(13)

    ϕ(x) =  x−x d  sign (x − x d ) · nd

(14)

where x d is the normal projection of point x on the interface, and nd is unit normal vector at x d . It must be noted that if the domain contains multiple cracks/interfaces, it is necessary to introduce different level set functions corresponding to each specified interface. It follows from the enriched displacement approximation that the displacement jump across the hydro-fracture and natural fault can be expressed by uh (x, t) = Nustd (x)¯a(t). Accordingly, the strain vector ε(x, t) corresponding to the approximate displacement field (11) can be written in terms of the standard and enriched nodal values as ti p ¯ εh (x, t) = Bstd ¯ + Badis (x)¯a(t) + Bb (x)b(t) u (x)u(t)

(15) std dis dis where Bstd u (x) = SNu (x), Ba (x) = SNa (x) and ti p ti p Bb (x) = SNb (x) are respectively the spatial derivatives of the standard, Heaviside enrichment and cracktip enrichment shape functions, and S is the well-known strain differential operator.

4.2 Discretization of the equilibrium equation The discrete form of the integral equation (4) can be obtained in the framework of X-FEM formulation using the test and trial functions of the enriched displacement and strain fields. In this manner, the test functions of

123

A. R. Khoei et al.

displacement and strain fields δuh (x, t) and δεh (x, t), given in the weak form Eq. (4), can be approximated within the same approximating space as uh (x, t) and εh (x, t) by ¯ + Nadis (x)δ a¯ (t) δuh (x, t) = Nustd (x)δ u(t) ti p ¯ + N (x)δ b(t) δε (x, t) = h

b ¯ + Badis (x)δ a¯ (t) Bstd u (x)δ u(t) ti p ¯ + Bb (x)δ b(t)

(16) (17)

and the jump of displacement variation can be obtained as δuh (x, t) = Nustd (x)δ a¯ (t). By substituting the enriched FE approximations into the weak form of equilibrium equation (4), rearranging the integral terms, and satisfying the necessity that the weak form should hold for all kinematically admissible test functions, it is straight forward to obtain the discretized system of nonlinear equations as ⎫ ⎫ ⎧ ⎧ Fuint − Fuext ⎪ ⎬ ⎨ ⎨ 1 ⎬ ⎪ int int int ext  =  2 = Fa − Fa + F HF + F fault = 0 ⎪ ⎭ ⎪ ⎩ ⎭ ⎩ 3 Fint − Fext b

b

(18) where

  T σ d Bstd u     T T ext std Nu Nustd t d ρb d + Fu =   T int dis σ d Ba Fa =      T T Nadis ρb d + Nadis t d Faext =   T ti p σ d Bb Fbint =     ti p T ti p T Nb Nb ρb d + t d Fbext =     T std N = tcont d F int u fault fault   T int Nustd pHF d F HF =

Fuint =

behavior of system is considered to be linearly elastic, the system of Eq. (18) is not linear due to the nonlinear constitutive law considered for the frictional contact behavior, as described in Sect. 3. However, it must be emphasized that even for the frictionless condition the problem is still kinematically nonlinear since the active contact zone along the natural fault with non-zero contact tractions is not known a priori. Hence, the Newton– Raphson iterative procedure needs to be employed to linearize the system of Eq. (18) at each iteration k of time interval (n, n + 1) as ⎡

⎤k ⎧ ⎫k+1 Kua Kua ⎨ d u¯ ⎬ d a¯ Kaa + Kcont Kab ⎦ ⎩ ⎭ Kba Kbb n+1 d b¯ n+1 ⎧ ⎫k Fuint − Fuext ⎨ ⎬ int = − Faint − Faext + F int + F HF fault ⎭ ⎩ Fbint − Fbext n+1

Kuu ⎣ Kau Kbu

in which the vector given in the right hand side of Eq. (20) stands for the residual force vector evaluated from ¯ d a¯ and d b¯ the known values of iteration k and, d u, denote the incremental values of unknown standard and enriched DOFs at iteration k + 1. Finally, the components of stiffness matrix in (20) can be derived as

Kua Kub Kau Kaa

Kab

In relations (18), the internal force vectors F int fault and F int HF are respectively expressed over the surfaces fault and HF , and take into account the internal forces emanated from the contact condition along the natural fault faces and the fluid pressure along the hydraulic fracture faces. It must also be noted that although the

Kbu

123

  T Bstd DBstd u u d    T = Bstd DBadis d u    T ti p = Bstd DBb d u   T = Badis DBstd u d    T ∂ 2 = Badis DBadis d + Kcont = ∂ a¯    T Nustd Dcont Nustd d (21) +   fault T ∂ 2 ti p Badis DBb d = = ∂ b¯   ∂ 3 ti p T = Bb = DBstd u d ∂ u¯    ∂ 3 ti p T = Bb = DBadis d ∂ a¯    ∂ 3 ti p T ti p = Bb = DBb d ∂ b¯  ∂ 1 ∂ u¯ ∂ 1 = ∂ a¯ ∂ 1 = ∂ b¯ ∂ 2 = ∂ u¯

Kuu =

(19)

HF

(20)

Kba Kbb

=

Modeling the interaction between fluid-driven fracture and natural fault

4.3 Modeling multiple cracks with junction enrichment function

+



H (ϕHF (x)) H (ϕfault (x))

for x ∈ A1 for x ∈ A2

(22)

in which A1 is the one side of the natural fault that contains the hydraulic fracture and A2 is the other side of the natural fault within an intersected element; ϕfault (x) and ϕHF (x) are the level set functions associated with the natural fault and hydraulic fracture respectively, as shown in Fig. 4. Thus, the extended finite element approximation of the displacement field (11) for an element containing the junction of two discontinuities of the natural fault and hydraulic fracture can be rewritten as uh (x, t) =



Ni (x)u¯ i (t) +

i∈N std



Ni (x) (H(ϕ(x))

i∈N dis

− H(ϕ(x i ))) a¯ i (t)

Fig. 4 An element containing the junction of hydraulic fracture and natural fault; description of the level set functions ϕfault and ϕHF

⎛ Ni (x) ⎝

i∈N ti p

In order to investigate the interaction between the hydraulic fracture and natural fault when the hydraulic fracture merges with the natural fault, a junction enrichment function is introduced to enrich all nodes whose support of the finite element shape functions contains the junction point. The implementation of the X-FEM for modeling of arbitrary branched and intersecting cracks with multiple branches, multiple holes and cracks emanating from holes was presented by Daux et al. (2000), Zi et al. (2004) and Budyn et al. (2004) by introducing a junction enrichment function. It must be noted that when a junction occurs within an element, as shown in Fig. 4, the crack tip enrichment strategy that provides the stress singularity at the tip of hydraulic fracture must be replaced by a junction enrichment function defined as

J (x) =



fault

0

+



4  

⎞ j B (x) − B (x i ) b¯ i (t)⎠ j

j

j=1

Ni (x) (J (x) − J (x i )) c¯ i (t)

in which N junc is the set of nodes that contain the junction in their supports, and c¯ i (t) are additional nodal DOFs corresponding to the junction enrichment function. Note that N dis is the set of nodes whose support is cut by both the hydro-fracture and natural fault discontinuities, and N ti p is the set of nodes that contain the tips of the natural fault in their supports. A schematic view of the enrichment strategy is illustrated in Fig. 5 at the junction of the hydro-fracture and natural fault. In order to model the interaction between the fluid-driven fracture and the natural fault, it is assumed that the tips of the natural fault are fixed and cannot proceed during the simulation. However, the tip of hydraulic fracture propagates with a constant length of a single element size at each step of crack propagation. In this manner, an appropriate criterion is required to detect the intersection of two discontinuities as the fluid-driven fracture propagates. An efficient criterion is chosen here for detection of the junction based on the procedure proposed by Budyn et al. (2004), in which the intersection of hydraulic fracture and natural fault is determined first by obtaining the closest segment of the natural fault, i.e. CB, from the tip of the hydraulic fracture at each propagation step, as depicted in Fig. 6. If the distance AD between the tip of hydraulic fracture and the segment CB is less than the size of a crack tip element, the hydraulic fracture is then connected to the fault. In such case, the crack tip asymptotic functions of the element containing the junction of hydro-fracture and natural fault is replaced by the junction enrichment function defined in (22) to follow the solution at the next increment.

0

HF

Hydraulic fracture

1

fault

0

(23)

i∈N junc

HF

0

2

Natural fault

123

A. R. Khoei et al.

Old hydro-fracture New hydro-fracture

Natural fault Standard nodes

Hydraulic fracture Heaviside enriched nodes Natural fault Heaviside enriched nodes Crack tip enriched nodes Double enriched nodes with junction and fault Heaviside

After junction

Before junction

Fig. 5 A schematic view of the enrichment strategy at the junction of hydro-fracture and natural fault

C Closest segment to hydro-fracture tip

New hydraulic fracture segment

D

A New junction

B

Natural fault

Old crack tip

Hydraulic fracture

Fig. 6 Detection of the junction between the natural fault and the hydraulic fracture

It must be noted that the junction enrichment function is a truncated form of the Heaviside function. In accordance with the junction enrichment function proposed by Daux et al. (2000), a modified junction function was proposed by Budyn et al. (2004) as  J (x) =

H (ϕHF (x)) 0

for x ∈ A1 for x ∈ A2

(24)

5 Numerical simulation results In order to illustrate the accuracy and robustness of the proposed computational algorithm, and to study the interaction between the fluid-driven fracture and the natural fault, several examples are solved numerically. Moreover, to provide a meaningful discus-

123

sion of the numerical simulation results, comparisons are performed between the X-FEM analyzes and those reported in the literature. For all numerical simulations, the elements located at the vicinity of discontinuities are partitioned into uniform subrectangles in which the standard Gauss quadrature rule is employed within each sub-rectangle. The surface integrals associated with the contact interfaces and hydro-fracture faces are performed using the onedimensional Gauss integration rule with two integration points along each interface segment, generated by intersection of the discontinuity with the FE mesh. For each Newton iteration, the contact stiffness matrix Kcont and the contact force vector F int fault are incorporated into the system of Eq. (20) if the set of Gauss points of the contact segment has a negative opening value, i.e. g N < 0, so that the inequality of g N ≥ 0 is fulfilled over the natural fault faces within each iteration until the convergence is achieved. In order to determine the direction of crack propagation, the maximum hoop stress criterion is employed by introducing the crack propagation angle in a direction normal to the maximum hoop stress. In this way, the angle of crack growth θcr is defined in terms of the stress intensity factors K I and K II as ⎞ ⎛    KI 2 1 ⎝ KI θcr = 2arctan − sign(K I I ) + 8⎠ 4 K II KI I

(25) which is computed with respect to x1 -axis (Fig. 2). Moreover, the pure mode I fracture toughness K I C is

Modeling the interaction between fluid-driven fracture and natural fault Fig. 7 An orthogonal hydro-fracture approaching a natural fault; the problem configuration and boundary conditions

B = 1m

Bi

E = 2 ×105 MPa ν = 0.3

A

H = 2.2 m

p B O

natural fault

y Bl Lc Li

utilized to determine the occurrence of crack growth at each increment when K I ≥ K I C . Note that the crack growth length is limited to the crack tip element size at each step of the propagation. Finally, the evaluation of stress intensity factors K I and K I I is carried out using a contour J -integral method.

5.1 Orthogonal hydro-fracture approaching a frictionless natural fault The first example is chosen to illustrate the performance of proposed computational algorithm for X-FEM modeling of a fluid-driven fracture approaching a frictionless natural fault. In order to investigate the interaction

Br

between the hydraulic fracture and a natural fault, a block containing a vertical natural fault of 2000 mm length and a horizontal propagating hydro-fracture is modeled using the proposed X-FEM computational model. The specimen is chosen 1000 × 2200 mm with a horizontal hydro-fracture at its center, as shown in Fig. 7. This example was originally simulated by Dyskin and Caballero (2009) using the finite element method in conjunction with zero thickness contact elements. In order to perform the numerical simulation, the hydraulic fracture propagation is carried out by a series of discrete simulation with different hydrofracture lengths. It is assumed that the position of the left hydro-fracture tip (point A) is fixed at a distance of Bl = 400 mm from the left edge, while the right

123

A. R. Khoei et al.

(a)

(b)

1.22

1.19

1.16 Li Li Li Li Li Li Li Li

1.13

1.1 0.0E+00

1.5E-04

= 70 - XFEM = 40 - XFEM = 10 - XFEM = 0 - XFEM = 70 - Dyskin and Caballero (2009) = 40 - Dyskin and Caballero (2009) = 10 - Dyskin and Caballero (2009) = 0 - Dyskin and Caballero (2009)

3.0E-04

4.5E-04

6.0E-04

Normal opening at the natural fault (m)

Y - coordinate (m)

Y - coordinate (m)

1.19

1.22 Li = 70 - XFEM Li = 40 - XFEM Li = 10 - XFEM Li = 0 - XFEM Li = 70 - Dyskin and Caballero (2009) Li = 40 - Dyskin and Caballero (2009) Li = 10 - Dyskin and Caballero (2009) Li = 0 - Dyskin and Caballero (2009)

1.16

1.13

1.1 -1.2

-0.9

-0.6

-0.3

0.0

Normalc ontact stress on the natural falut (MPa)

Fig. 8 An orthogonal hydro-fracture approaching the natural fault; a the profiles of normal opening and, b the profiles of normal contact stress along the natural fault

hydro-fracture tip (point B) approaches the natural fault considering four values of L i = 70, 40, 10 mm, and zero, that makes possible to investigate the interaction between the natural fault and the hydraulic fracture, as the hydro-fracture tip propagates and merges with the natural fault. To remain consistent with the work given by Dyskin and Caballero (2009), the hydro-fracture inflow is modeled by imposing a uniform constant pressure of p = 1 MPa on the faces of hydro-fracture through the propagation process. Moreover, the contact behavior along the natural fault is modeled using the penalty X-FEM method described in Sect. 3, employing a normal penalty parameter k N = 1010 MPa/mm. The X-FEM analysis is performed employing a quadrilateral structured mesh of 3025 elements with 3136 nodal points, which is concentrated over the central region of the domain containing the hydro-fracture and the natural fault. In Fig. 8, the profiles of normal opening g N and normal contact stress t Ncont are plotted along the natural fault for different values of L i , and the results are compared with those reported by Dyskin and Caballero (2009). Note that the results are only plotted for one half of the natural fault at the upper part due to the symmetry of the specimen. A great agreement can be observed between the current X-FEM simulation and those reported in the literature. Obviously, the interaction between the hydraulic fracture and the nat-

123

ural fault results in inducing both opening and closing modes along the natural fault. As the hydro-fracture tip approaches the natural fault, the central opening zone induced in front of the hydro-fracture tip decreases and the normal opening displacement increases at the tips of natural fault. When the hydro-fracture tip merges with the fault, i.e. L i = 0, the overall behavior of the system changes dramatically; the central opening zone disappears and the fault experiences the closing mode in the vicinity of an intersection point O leaving the remaining parts of the interface in opening mode. Note that the magnitude of the normal contact stress is almost in the same order of the hydraulic fracturing fluid pressure. However, the length of the natural fault that was affected by the fracturing fluid pressure, which is equal to the hydro-fracture opening at the intersection point O, is virtually negligible compared to the length of contacting faces of the natural fault. In this manner, it can be predicted that the natural fault is, at least momentarily, capable to prevent the flow of fracturing fluid within the natural fault once the hydro-fracture tip merges with the fault. This prediction is in accordance with the criteria reported by Potluri et al. (2005) based upon which the natural fault dilates only if the fracturing fluid pressure is greater than the normal contact tractions on the natural fault faces. In Fig. 9, the profiles of normal stress σ y are plotted at the right side of the natural fault for different values of L i , and the results are compared with

Modeling the interaction between fluid-driven fracture and natural fault

(b)

2.1 Li = 70 - XFEM Li = 40 - XFEM Li = 10 - XFEM Li = 0 - XFEM Li = 70 - Dyskin and Caballero (2009) Li = 40 - Dyskin and Caballero (2009) Li = 10 - Dyskin and Caballero (2009) Li = 0 - Dyskin and Caballero (2009)

Y - coordinate (m)

1.9

1.7

1.5

1.3

1.1 -1.2

1.22 Li Li Li Li Li Li Li Li

1.2

Y - coordinate (m)

(a)

1.18

= 70 - XFEM = 40 - XFEM = 10 - XFEM = 0 - XFEM = 70 - Dyskin and Caballero (2009) = 40 - Dyskin and Caballero (2009) = 10 - Dyskin and Caballero (2009) = 0 - Dyskin and Caballero (2009)

1.16

1.14

1.12

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

0.2

0.4

σ y on the right side of the natural fault (MPa)

1.1 -1.2

-1.0

-0.8

-0.6

-0.4

-0.2

0.0

0.2

σ y on the right side of the natural fault (MPa)

Fig. 9 The distributions of normal stress σ y at the right side of natural fault; a the profiles of stress distribution for the upper part of the specimen, b the graphs zoomed in the vicinity of point O

those given by Dyskin and Caballero (2009). It can be seen from the stress curves that once the hydro-fracture tip merges with the natural fault, the stress component σ y becomes compressive in front of point O at the right side of natural fault so that it can prevent the hydrofracture to cross the natural fault. Moreover, the stress σ y becomes tensile in the upper part of natural fault and reaches its maximum value at the tip of the fault representing a singular behavior at this point. As a result, an offset crack can be initiated at the tips of natural fault. This implies that a sufficiently large frictionless natural fault is capable of arresting the fluid-driven fracture and that the offset crack can propagate from the ends of the fault. In Fig. 10, the contours of horizontal stress σx are shown for two hydro-fracture configurations, i.e. L i = 40 mm and L i = 0. The stress contours clearly illustrate the mechanism of fault closure in front of the hydro-fracture tip when a junction takes place at point O. In the case of L i = 40 mm, where the hydrofracture is not connected to the natural fault, the horizontal stress σx induced in front of the hydro-fracture tip is in tensile. This tensile stress in normal direction to the natural fault results in an opening mode on the natural fault in the vicinity of an intersection point O. In contrast, when the hydro-fracture intersects with the natural fault, i.e. L i = 0, the horizontal stress in front of point O becomes compressive leading to fault closure in the vicinity of this point. It must be noted that

although the results reported here provide an insight into the mechanism of interaction between the fluiddriven fracture and a long frictionless natural fault, the validity of conclusion is limited to current configuration of the natural fault and the imposed boundary conditions.

5.2 Hydraulic fracture propagation subjected to a frictional natural fault The next example is chosen to study the mechanism of interaction between the fluid-driven fracture and a frictional natural fault in an infinite domain, as shown in Fig. 11. In this study, the effects of various fault lengths and different frictional resistance along the natural fault are investigated on the mechanism of interaction between two discontinuities. The material parameters of the domain are assumed as the Young modulus E = 15.96 GPa and Poisson ratio ν = 0.33. The geometry and boundary conditions of the domain together with the finite element mesh are shown in Fig. 11. Note that the chosen material properties, geometry and boundary conditions are consistent with a benchmark example studied originally by Schrefler et al. (2006) and Khoei et al. (2015), where the hydraulic fracture propagation was investigated in an infinite medium. Due to symmetry of the problem, it is evident that the hydro-fracture propagates in a horizontal direction

123

A. R. Khoei et al.

x

Li

40 mm

Li

0

(at the intersecon of hydrofracture and natural fault)

Fig. 10 The contours of horizontal stress σx for an orthogonal crack approaching a natural fault at L i = 40 mm and L i = 0, where the junction occurs (color bar in Pa)

under pure mode I crack propagation. The hydraulic fracture propagation is performed by applying a uniform constant pressure p = 1 MPa on the faces of hydro-fracture through the propagation process. The hydro-fracture propagates in a horizontal direction for the total length of HF = 5 m until it merges with a pre-existing vertical natural fault. In order to investigate the effect of various fault lengths on the mechanism of interaction between two discontinuities, four lengths of fault = 10, 20, 30, and 40 m are assumed for the natural fault. Moreover, in order to study the effect of frictional behavior along the natural fault, three frictional coefficients of μ f = 0, 0.2, and 0.4 are employed. For all numerical simulations, the normal and tangential penalty parameters are assumed as k N = k T = 109 MPa/m. The numerical simulation is performed employing an X-FEM mesh of 2420 quadrilateral elements with 2448 nodes, in which most of the elements are clustered in the vicinity of fracture propagation trajectory, as shown in Fig. 11. In Fig. 12, the profiles of normal opening g N are plotted along the natural fault at different hydrofracture growths HF for various fault lengths fault and frictional coefficients μ f . Obviously, as the hydraulic fracture advances through the domain, an opening mode is induced along the natural fault with a maximum opening value at its center. The magnitude of this opening increases for larger values of fault length. Once

123

the hydro-fracture merges with the natural fault, i.e. HF = 5 m, the overall behavior of the system changes significantly, such that the normal opening displacement induced along the natural fault vanishes, leaving the entire length of the natural fault in contact conditions irrespective of the length of natural fault. In Fig. 13, the profiles of normal stress σ y are plotted along the natural fault at the right side of the fault for various hydro-fracture growths HF using different values of fault length fault and frictional coefficients μ f . In contrast to the previous example where the vertical stress σ y was in compression ahead of intersection point of the hydro-fracture and the natural fault, it can be seen here the development of both tensile and compressive vertical stresses ahead of intersection point, depending on the fault length and frictional resistance along the natural fault. Obviously, for the shortest length of natural fault, i.e. fault = 10 m, the stress σ y is in tensile ahead of intersection point irrespective of frictional resistance along the natural fault. However, the compressive vertical stress can be observed at the intersection area for larger values of fault length, i.e. fault = 20, 30, and 40 m. It is noteworthy to highlight that the larger value of frictional resistance results in a larger tensile stress ahead of intersection point of two discontinuities. Accordingly, it can be concluded that the lower value of fault length together with a larger frictional resistance along the natural fault produces a larger vertical

Modeling the interaction between fluid-driven fracture and natural fault Fig. 11 Hydraulic fracturing subjected to a frictional natural fault; the problem definition, boundary conditions and X-FEM mesh

10 m

80 m

ux = 0

160 m

p

fault

HF

ux = 0 uy = 0

5m

tensile stress ahead of intersection point of two discontinuities, and makes the possibility of penetration of the hydro-fracture through the natural fault. In contrast, a larger value of fault length together with the lower frictional resistance along the natural fault has the potential of producing the compressive vertical stress ahead of intersection point of two discontinuities that makes possible to arrest the hydro-fracture propagation. The above argument is in general in accordance with the results reported in the literature based on experimental observations. Daneshy (1978) and Hanson et al. (1980) pointed out that the shear strength along the interfaces affect the mechanism of interaction between the hydro-fracture and frictional natural fault. Their experiments illustrated that the more the frictional resistance of the natural interface, the higher the possibility of the hydro-fracture penetration through the natural fault.

In Fig. 14, the contours of vertical stress σ y are presented for various fault lengths fault and frictional coefficients μ f along the natural fault when the hydrofracture merges with the natural fault, i.e. HF = 5 m. These contours clearly present the distribution of compressive and tensile vertical stresses at the right side of the natural fault, which are generally in accordance with those plotted in Fig. 13. Moreover, it can be observed that the tensile vertical stress can be produced at the tips of natural fault; it illustrates how two tips of the natural fault can be considered as an obvious candidate for the offset crack propagation. 5.3 Mixed-mode hydro-fracture propagation in a naturally fractured domain In the last example, the mixed-mode hydro-fracture propagation is simulated in a naturally fractured domain

123

A. R. Khoei et al. 10

5

Y- coordinate (m)

= 10 m

0

-2.5

-5 0.00

fault

5

Y-coordinate (m)

fault

2.5

= 20 m

0

HF

= 3 m (μ = 0.0)

HF

= 3 m (μ = 0.4)

HF

= 4 m (μ = 0.0)

HF

= 4 m (μ = 0.4)

HF

= 5 m (μ = 0.0)

HF

= 5 m (μ = 0.4)

= 3 m (μ = 0.2)

HF

= 3 m (μ = 0.2)

HF

= 4 m (μ = 0.2)

HF

= 4 m (μ = 0.2)

HF

= 5 m (μ = 0.2)

HF

= 5 m (μ = 0.2)

HF

= 3 m (μ = 0.0)

HF

= 3 m (μ = 0.4)

HF

= 4 m (μ = 0.0)

HF

= 4 m (μ = 0.4)

HF

= 5 m (μ = 0.0)

HF

= 5 m (μ = 0.4)

HF

0.05

0.10

0.15

0.20

-5

-10

0.25

0.00

0.05

0.10

0.15

0.20

0.25

Normal opening at the natural fault (mm)

Normal opening at the natural fault (mm) 15

20

10

fault

= 30 m

fault

= 40 m

5

0

-5

-10

-15 0.00

0.05

Y-coordinate (m)

Y-coordinate (m)

10

0

HF

= 3 m (μ = 0.0)

HF

= 3 m (μ = 0.4)

HF

= 4 m (μ = 0.0)

HF

= 4 m (μ = 0.4)

HF

= 5 m (μ = 0.0)

HF

= 5 m (μ = 0.4)

= 3 m (μ = 0.2)

HF

= 3 m (μ = 0.2)

HF

= 4 m (μ = 0.2)

HF

= 4 m (μ = 0.2)

HF

= 5 m (μ = 0.2)

HF

= 5 m (μ = 0.2)

HF

= 3 m (μ = 0.0)

HF

= 3 m (μ = 0.4)

HF

= 4 m (μ = 0.0)

HF

= 4 m (μ = 0.4)

HF

= 5 m (μ = 0.0)

HF

= 5 m (μ = 0.4)

HF

0.10

0.15

0.20

0.25

Normal opening at the natural fault (mm)

-10

-20 0.00

0.05

0.10

0.15

0.20

0.25

Normal opening at the natural fault (mm)

Fig. 12 Hydraulic fracturing subjected to a frictional natural fault; the profiles of normal opening along the natural fault at different hydro-fracture growths HF for various lengths fault and frictional coefficients μ f of the natural fault

where different fault orientations are employed. The problem configuration together with the geometry and boundary conditions are presented in Fig. 15. A block of 2.5 × 5.0 m containing a vertical hydro-fracture of initial length HF = 0.5 m and a frictional natural fault of fixed length fault = 2.0 m with the frictional coefficient μ f = 0.2 is modeled using the Young modulus E = 2 × 104 MPa and Poisson ratio ν = 0.2. The hydro-fracture faces are subjected to a uniform internal pressure of 3 p = 29.1 MPa through the propagation process. In this study, the effect of fault orientations is investigated for three values of fault angle β = 30◦ , 60◦ and 75◦ . Moreover, in order to study the effect of

123

far-field stress conditions on the mechanism of interaction between the hydro-fracture and natural fault, four different conditions of far-field horizontal and vertical stresses σh and σv are investigated, as shown in Fig. 15. The X-FEM analyzes are performed using a structured FE mesh consisting of 2700 bilinear elements with 2828 nodal points, which are concentrated over the center of the block with respect to the horizontal direction. The frictional contact condition is modeled along the natural fault using the normal and tangential penalty parameters k N = k T = 1012 MPa/m. The problem configuration and material properties of this example are in accordance with Dong and De Pater (2001).

Modeling the interaction between fluid-driven fracture and natural fault 10

5

= 10 m

Y-coordinate (m)

2.5

0

-2.5

HF

= 3 m (μ = 0.0)

HF

= 4 m (μ = 0.0)

HF

= 5 m (μ = 0.0)

HF

= 3 m (μ = 0.2)

HF

= 4 m (μ = 0.2)

HF

= 5 m (μ = 0.2)

HF

= 3 m (μ = 0.4)

HF

= 4 m (μ = 0.4)

HF

= 5 m (μ = 0.4)

-5 -2.0

-1.0

fault

5

Li==3 3mm(μ( =u0.0) =

0

-5

0.0

HF

= 4 m (μ = 0.0)

HF

= 5 m (μ = 0.0)

HF

= 3 m (μ = 0.2)

HF

= 4 m (μ = 0.2)

HF

= 5 m (μ = 0.2)

HF

= 3 m (μ = 0.4)

HF

= 4 m (μ = 0.4)

HF

= 5 m (μ = 0.4)

-10 -2.0

1.0

σ y on the right side of the natural fault (MPa)

-1.0

0.0

1.0

σ y on the right hand side of the natural fault (MPa)

15

20 fault

10

= 30 m

fault

10

Li ==33mm(μ( =u0.0) =

0

-5

HF

= 4 m (μ = 0.0)

HF

= 5 m (μ = 0.0)

HF

= 3 m (μ = 0.2)

HF

= 4 m (μ = 0.2)

HF

= 5 m (μ = 0.2)

HF

= 3 m (μ = 0.4)

HF

= 4 m (μ = 0.4)

HF

= 5 m (μ = 0.4)

HF

Y-coordinate (m)

5

= 40 m Li==33mm(μ( =u0.0) =

HF

Y-coordinate (m)

= 20 m HF

Y-coordinate (m)

fault

0

-10

HF

= 4 m (μ = 0.0)

HF

= 5 m (μ = 0.0)

HF

= 3 m (μ = 0.2)

HF

= 4 m (μ = 0.2)

HF

= 5 m (μ = 0.2)

HF

= 3 m (μ = 0.4)

HF

= 4 m (μ = 0.4)

HF

= 5 m (μ = 0.4)

-10

-15 -2.0

-1.0

0.0

1.0

σ y on the right hand side of the natural fault (MPa)

-20 -2.0

-1.0

0.0

1.0

σ y on the right hand side of the natural fault (MPa)

Fig. 13 Hydraulic fracturing subjected to a frictional natural fault; the profiles of normal stress σ y at the right hand side of the natural fault at different hydro-fracture growths HF for various lengths fault and frictional coefficients μ f of the natural fault

In Fig. 16, the crack trajectories are presented for the fluid-driven fracture propagation at different fault angles using various loading conditions. Obviously, the fault does not affect the crack trajectory as long as the hydro-fracture tip is quite far from the natural fault such that the hydraulic fracture propagates in an almost straight path with respect to its initial orientation. The effect of natural fault, however, becomes significant when the hydraulic fracture tip approaches the natural fault and moves toward a normal direction to the fault. This is in accordance with experimental observations reported by Lamont and Jessen (1963), in

which a greater inclination angle of the natural fault with respect to the hydro-fracture leads to a further deviation of the hydro-fracture path from the center line of the model. It can be clearly observed from these figures that a larger value of far-field horizontal stress causes a greater deviation from the straight path of crack propagation. It is obvious that the farfield horizontal stress has a significant effect on the performance of internal pressure imposed on the hydrofracture faces, and plays an important role on the mechanism of interaction between the hydro-fracture and natural fault. In fact, the overall behavior of the sys-

123

A. R. Khoei et al.

μ f = 0.0

μ f = 0.2

μ f = 0.4

μ f = 0.0

μ f = 0.2

(a)

fault

= 10 m

(b)

fault

= 20 m

(c)

fault

= 30 m

(d)

fault

= 40 m

-1.5E+06

-7.5E+05

0.0E+00

7.5E+05

μ f = 0.4

1.5E+06

Fig. 14 The contours of vertical stress σ y for various fault length fault and frictional coefficients μ f along the natural fault when the hydro-fracture merges with the natural fault HF = 5 m (color bar in Pa)

tem is mainly governed by the difference between the hydro-fracture internal pressure and far-field horizontal stress. In Fig. 17, the contours of vertical stress σ y are shown for different loading conditions using three values of fault angles when the hydro-fracture merges with the natural fault. It can be seen from this figure that the stress contours are almost similar in loading cases 1 and 2 where the same far-field horizontal stress is employed. Similarly, the stress contours of loading conditions 3 and 4 display almost the same behavior. It can therefore be highlighted that the crack trajectory for a fractured domain with farfield stress conditions is mainly affected by the difference between the hydro-fracture internal pressure and far-field stress conditions. It was reported through experimental observations by Blanton (1982) that the hydro-fracture penetration across the natural fault can be happened under the high differential stresses and

123

high inclination angles. This is again in accordance with the results presented in Fig. 17 that shows the vertical stress σ y is always in tensile for the greatest inclination angle of β = 75◦ for all loading conditions, indicative of the possibility of hydro-fracture penetration in such case. In Fig. 18, the profiles of fault opening displacement are plotted along the natural fault at final stage of hydrofracture propagation using various loading conditions for three values of fault angles. Also plotted in Fig. 19 are the profiles of fault normal contact stress along the natural fault at final stage of hydro-fracture propagation for three values of fault angles. Obviously, an opening displacement mode is formed along the natural fault for all three cases at the intersection point of two discontinuities. Hence, it can be concluded that the fluid flow within the hydro-fracture can permeate through the frictional natural fault, proposing the feasibility of

Modeling the interaction between fluid-driven fracture and natural fault Fig. 15 The mixed mode hydro-fracture propagation in a naturally fractured domain; the problem configuration and boundary conditions

v

fault

= 2m

= 0.5 m HF

H = 5m

h

y

d =2m

O

s

h 3p

Loading case

Horizontal stress

Vertical stress

h

p

v

p

Case 2.

h

p

v

2p

Case 3.

h

2p

v

p

h

2p

v

2p

Case 1.

Case 4.

v B = 2.5 m Fig. 16 The crack trajectories for different fault angles at various loading conditions

fluid diversion condition in this particular mixed mode hydraulic fracturing problem. It is evident from this example that the results obtained here do not necessarily follow the results observed in previous examples. This is because there are a wide range of parameters that may affect the overall behavior of the interaction mechanism, including the hydraulic fracture/natural fault configuration, the fault inclination angle, far-field stress conditions and the frictional resistance along the natural fault, which clearly indicates the necessity of employing a robust computational algorithm for mod-

eling the interaction problem in hydraulic fracturing processes.

6 Conclusion In the present paper, an extended finite element method was presented for modeling the hydro-fracture propagation in naturally fractured domain, with a particular attention to the mechanism of interaction between the hydro-fracture and natural fault. The discontinu-

123

A. R. Khoei et al. Fig. 17 The contours of stress σ y for the mixed mode hydro-fracture propagation in a naturally fractured domain for different loading conditions; a β = 30◦ , b β = 60◦ , c β = 75◦

ities within the domain were modeled by introducing additional degrees of freedom at the nodal points of the FE mesh circumventing the need for remeshing during the hydro-fracturing process. The X-FEM method was utilized to discretize the governing equations of a fractured domain containing a frictional natural fault and a hydro-fracture interface subjected to a uniform constant pressure along its faces. In order to model the frictional contact behavior along the natural fault during the hydro-fracturing process, the X-FEM penalty method was employed by imposing the normal contact constraint to prevent the penetration between two faces of the natural fault in the normal direction, and applying a plasticity theory of friction to capture the stick-slip behavior along the contact faces in the tangential direc-

123

tion. The crack tip asymptotic functions were utilized to enrich the displacement fields at the tips of hydrofracture and natural fault; the Heaviside jump function was used to model the displacement jump along the hydro-fracture and natural fault faces; and a junction enrichment function was employed to capture the intersection between the hydro-fracture and natural fault. In order to illustrate the robustness and versatility of proposed computational algorithm several examples were solved numerically. The first example was chosen to illustrate the performance of proposed computational model for hydraulic fracturing process, in which the hydro-fracture faces are subjected to a uniform constant pressure. A good agreement was shown between the proposed X-FEM modeling and those reported in

Modeling the interaction between fluid-driven fracture and natural fault

(b)

8

Normal opning at the natural fault (mm)

Normal opning at the natural fault (mm)

(a)

Loading case 1 Loading case 2 Loading case 3 Loading case 4

6

4

2

0

M3 M4

0

M1

0.5

M2

8

Loading case 1 Loading case 2 Loading case 3 Loading case 4

6

4

2

0 1

1.5

2

0

M3 M4

M2

0.5

M1

S - coordinate

Normal opning at the natural fault (mm)

(c)

1

1.5

2

S - coordinate

8

Loading case 1 Loading case 2 Loading case 3 Loading case 4

6

4

2

0 0

M3 M4

M2 M1

0.5

1

1.5

2

S - coordinate Fig. 18 Profiles of the fault opening displacement at final stage of hydro-fracture propagation using various loading conditions; a β = 30◦ , b β = 60◦ , c β = 75◦ (Mi denotes the location of junction along S-direction)

the literature. In this example, the interaction between the fluid-driven fracture and a long frictionless natural fault was investigated. It was shown that a relatively long frictionless and cohesionless fault is capable of arresting the fluid-driven fracture; it was also observed that the offset cracks may propagate from the tips of natural fault. The second example was chosen to investigate the effects of fault length and frictional resistance along the natural fault on the mechanism of interaction between the fluid-driven fracture and nat-

ural fault. It was highlighted that the lower value of fault length together with a larger frictional resistance along the natural fault produces a larger vertical tensile stress ahead of intersection point of two discontinuities, and makes the possibility of penetration of the hydro-fracture through the natural fault. Furthermore, it was illustrated that a larger value of fault length together with the lower frictional resistance along the natural fault has the potential of producing the compressive vertical stress ahead of intersection point of

123

A. R. Khoei et al.

(b)

50 Loading case 1 Loading case 2 Loading case 3 Loading case 4

40

30

20

10

0

0

0.5

1

1.5

50

Normal contact stress on the fault (MPa)

Normal contact stress on the fault (MPa)

(a)

Loading case 3 Loading case 4

30

20

10

0

2

Loading case 1 Loading case 2

40

0

0.5

Normal contact stress on the fault (MPa)

(c)

1

1.5

2

S - coordinate

S - coordinate

50 Loading case 1 Loading case 2 Loading case 3 Loading case 4

40

30

20

10

0

0

0.5

1

1.5

2

S - coordinate Fig. 19 Profiles of the fault normal contact stress at final stage of hydro-fracture propagation using various loading conditions; a β = 30◦ , b β = 60◦ , c β = 75◦

two discontinuities that makes possible to arrest the hydro-fracture propagation. Finally, the last example was chosen to illustrate the performance of proposed computational algorithm for X-FEM modeling of the mixed mode hydro-fracture propagation in a naturally fractured domain, in which the effects of fault orientation and far-field stress conditions were investigated on the mechanism of interaction between the fluid-driven fracture and natural fault. It was observed that the farfield stress conditions have a significant effect on the performance of internal pressure imposed on the hydro-

123

fracture faces, and plays an important role on the mechanism of interaction between the hydro-fracture and natural fault. Moreover, it was concluded that there are a wide range of parameters that may affect the overall behavior of the interaction mechanism, including the hydraulic fracture/natural fault configuration, the fault inclination angle, far-field stress conditions, and the frictional resistance along the natural fault. Based on the results presented here, it has been shown that the proposed X-FEM computational model is a desirable candidate in tackling hydraulic fracturing problems in

Modeling the interaction between fluid-driven fracture and natural fault

even more complex conditions. However, the proposed model can be further improved by incorporating several complicated mechanisms involved in the interaction between the hydro-fracture and natural fault. For instance, one may investigate the after-junction mechanisms by considering the fracturing fluid diversion into the natural fault after coalescence of the hydraulic fracture with the natural fault, and/or considering the effect of fluid leakage flow from the fracture faces into surrounding medium. Indeed, simulation of such phenomena is not considered here and is left for a future study. Acknowledgments The senior author is grateful for the research support of the Iran National Science Foundation (INSF).

References Akulich AV, Zvyagin AV (2008) Interaction between hydraulic and natural fractures. Fluid Dyn 43:428–435 Anderson GD (1981) Effects of friction on hydraulic fracture growth near unbonded interfaces in rocks. Soc Petrol Eng 21:21–29 Belytschko T, Black T (1999) Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 45:601–620 Blanton TL (1982) An experimental study of interaction between hydraulically induced and pre-existing fractures. Soc Petrol Eng. doi:10.2118/10847-MS Boone TJ, Ingraffea AR (1990) A numerical procedure for simulation of hydraulically-driven fracture propagation in poroelastic media. Int J Numer Anal Meth Geomech 14:27– 47 Brenner SL, Gudmundsson A (2004) Arrest and aperture variation of hydrofractures in layered reservoirs. Geo Soc Lond 231:117–128 Budyn E, Zi G, Moës N, Belytschko T (2004) A method for multiple crack growth in brittle materials without remeshing. Int J Numer Meth Eng 61:1741–1770 Chopp DL, Sukumar N (2003) Fatigue crack propagation of multiple coplanar cracks with the coupled extended finite element/fast marching method. Int J Eng Sci 41:845–869 Dahi A, Olson JE (2011) Numerical modeling of multistrandedhydraulic-fracture propagation: accounting for the interaction between induced and natural fractures. Soc Petrol Eng 16:575–581 Dahm T (2000) Numerical simulations of the propagation path and the arrest of fluid-filled fractures in the earth. Geophys J Int 141:623–638 Daneshy AA (1974) Hydraulic fracture propagation in the presence of planes of weakness. Soc Petrol Eng. Paper no 4852 Daneshy AA (1978) Hydraulic fracture propagation in layered formations. Soc Petrol Eng 18:33–41 Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T (2000) Arbitrary branched and intersecting cracks with the extended finite element method. Int J Numer Meth Eng 48:1741–1760

Dong CY, De Pater CJ (2001) Numerical implementation of displacement discontinuity method and its application in hydraulic fracturing. Comp Meth Appl Mech Eng 191:745– 760 Dyskin AV, Caballero A (2009) Orthogonal crack approaching an interface. Eng Fract Mech 76:2476–2485 Gudmundsson A, Brenner SL (2001) How hydrofractures become arrested. Terra Nova 13:456–462 Hanson ME, Anderson GD, Shaffer RJ (1980) Theoretical and experimental research on hydraulic fracturing. J Energ Resour ASME 102:92–98 Hirmand M, Vahab M, Khoei AR (2015) An augmented Lagrangian contact formulation for frictional discontinuities with the extended finite element method. Finite Elem Anal Des. doi:10.1016/j.finel.2015.08.003 Jeffrey RG, Settari A (1995) A comparison of hydraulic fracture field experiments, including mineback geometry data, with numerical fracture model simulations. Soc Petrol Eng. doi:10.2118/30508-MS Khoei AR (2015) Extended finite element method, theory and applications. Wiley, New York Khoei AR, Biabanaki SOR, Anahid M (2009) A Lagrangianextended finite-element method in modeling large-plasticity deformations and contact problems. Int J Mech Sci 51:384– 401 Khoei AR, Hirmand M, Vahab M, Bazargan M (2015) An enriched FEM technique for modeling hydraulically driven cohesive fracture propagation in impermeable media with frictional natural faults; Numerical and experimental investigations. Int J Numer Meth Eng. doi:10.1002/nme.4944 Khoei AR, Nikbakht M (2007) An enriched finite element algorithm for numerical computation of contact friction problems. Int J Mech Sci 49:183–199 Khoei AR, Vahab M (2014) A numerical contact algorithm in saturated porous media with the extended finite element method. Comput Mech 54:1089–1110 Khoei AR, Vahab M, Haghighat E, Moallemi S (2014) A meshindependent finite element formulation for modeling crack growth in saturated porous media based on an enriched FEM technique. Int J Fract 188:79–108 Lamont N, Jessen FW (1963) The effects of existing fractures in rocks on the extension of hydraulic fractures. J Pet Technol 15:203–209 Li J (2000) Debonding of the interface as ‘crack arrestor’. Int J Fract 105:57–79 Liu F, Borja RI (2008) A contact algorithm for frictional crack propagation with the extended finite element method. Int J Numer Meth Eng 76:1489–1512 Melenk JM, Babuška I (1996) The partition of unity finite element method: basic theory and applications. Comp Meth Appl Mech Eng 139:289–314 Moës N, Dolbow J, Belytschko T (1999) A finite element method for crack growth without remeshing. Int J Numer Meth Eng 46:131–150 Mohammadnejad T, Khoei AR (2013) An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model. Finite Elem Meth Anal Des 73:77–95 Murphy HD, Fehler MC (1986) Hydraulic fracturing of jointed formations. Soc Petrol Eng. doi:10.2118/14088-MS

123

A. R. Khoei et al. Potluri N, Zhu D, Hill A (2005) The effect of natural fractures on hydraulic fracture propagation. Soc Petrol Eng 94568:1–6 Price NJ, Cosgrove JW (1990) Analysis of geological structures. Cambridge University Press, Cambridge Réthoré J, de Borst R, Abellan MA (2007) A two-scale approach for fluid flow in fractured porous media. Int J Numer Meth Eng 71:780–800 Schrefler BA, Secchi S, Simoni L (2006) On adaptive refinement techniques in multi-field problems including cohesive fracture. Comp Meth Appl Mech Eng 195:444–461 Stolarska M, Chopp DL, Moës N, Belytschko T (2001) Modelling crack growth by level sets in the extended finite element method. Int J Numer Meth Eng 51:943–960 Warpinski NR, Teufel LW (1987) Influence of geologic discontinuities on hydraulic fracture propagation. J Petrol Technol 39:209–220

123

Zhang X, Jeffrey RG, Thiercelin M (2007) Deflection and propagation of fluid-driven fractures at frictional bedding interfaces: a numerical investigation. J Struct Geol 29:396–410 Zhang Z, Ghassemi A (2011) Simulation of hydraulic fracture propagation near a natural fracture using virtual multidimensional internal bonds. Int J Numer Anal Meth Geomech 35:480–495 Zhou J, Chen M, Jin Y, Zhang G (2008) Analysis of fracture propagation behavior and fracture geometry using a triaxial fracturing system in naturally fractured reservoirs. Int J Rock Mech Mining Sci 45:1143–1152 Zi G, Song JH, Budyn E, Lee SH, Belytschko T (2004) A method for growing multiple cracks without remeshing and its application to fatigue crack growth. Model Simul Mat Sci Eng 12:901–915