International Journal of Fracture 96: 127–147, 1999. © 1999 Kluwer Academic Publishers. Printed in the Netherlands.
The effective fracture toughness in hydraulic fracturing PANOS PAPANASTASIOU Schlumberger Cambridge Research, High-Cross, Madingley Road, Cambridge, CB3 0EL, United Kingdom. e-mail:
[email protected] Received 19 February 1998; accepted in revised form 13 November 1998 Abstract. This paper examines the effective fracture toughness approach which is used in hydraulic fracturing in order to explain the high net-pressures that are often observed in field operations. The effective fracture toughness is calculated using a fully deterministic elasto-plastic hydraulic fracturing model. Rock is modelled by Mohr– Coulomb flow theory of plasticity for cohesive-frictional dilatant material. Fluid flow is modelled by lubrication theory. A cohesive crack model which takes into account the softening behaviour of rocks is employed as the propagation criterion. The fully coupled model is solved numerically by the finite element method and the effective fracture toughness is calculated using the path independent J -integral. The results show that plastic yielding near the tip of a propagating fracture provides an effective shielding, resulting in an increase in the rock effective fracture toughness by more than an order of magnitude. It is demonstrated that an elastic model based on the concept of effective fracture toughness matches the results of plasticity quite well. The effective fracture toughness increases with formation yielding, which is influenced by the deviator of the in-situ stresses, the rock strength, the elastic modulus and the pumping parameters. Tables of effective fracture toughness for a representative set of physical parameters are presented. Key words: Hydraulic fracturing, fracture toughness, net pressure, plasticity, coupled model, finite elements.
1. Introduction The problem of a hydraulically fluid-driven fracture has attracted numerous contributions since the early fifties (Kristianovitch and Zheltov, 1955; Geertsma and DeKlerk, 1969). Interest in such studies rises mainly from hydraulic fracturing, a technique for stimulating hydrocarbon reservoirs but also from other applications such as magma driven fracture (Spence and Turcotte, 1985). New applications of hydraulic fracturing have been recently found in geotechnical engineering for ground reinforcement, in petroleum engineering for reinjection of drilling cuttings and in environmental engineering for solids waste disposal. Modelling the propagation of hydraulic fractures is usually carried out prior to fracturing in order to optimize the treatment. In practice, a lot of attention is focussed on the prediction of wellbore pressure; wellbore pressure is normally measured during the treatment and is usually the only parameter available to evaluate the operation. Classical hydraulic fracturing simulators are based on linear elasticity and often underestimate the down-hole pressures which are measured in field operations. A recent world-wide survey on net-pressures (difference between fracturing pressure and far-field confining stress) by the Delft Fracturing Consortium indicated that net-pressures encountered in the field are 50 to 100 percent higher than the netpressures predicted by the conventional hydraulic fracturing simulator. The difference is even higher in weak formations (Van Dam et al., 1997). Several hypotheses have been proposed to explain the discrepancy. Among them the most consistent with observations are (1) high
128 Panos Papanastasiou friction losses in constrictions within a fracture and (2) effective fracture toughness effects (e.g. due to micro-cracking or plasticity in the process zone). Two main approaches have been proposed in order to explain the discrepancies between models and field observations. According to the first approach, the observed high net pressure is related to a sharp drop of fluid-pressure and the existence of a dry region near the crack-tip (Cleary et al., 1991). Cleary et al. (1991) claimed that elasticity can not provide the required high pressure drop to match the observed net pressure. They stated that other mechanisms must be introduced in the models, and proposed that rock dilation may be the source of high pressure drop near the tip. According to the dilation hypothesis, the rock dilation behind the advancing fracture would constrain the opening which may lead to the sharp pressure gradients. Barr (1991) and Gardner (1992) investigated the possible effect of dilatancy on the net-pressures using a fully-coupled elastic simulator, in which the dilatancy was modelled rather arbitrary by artificially constricting the fracture opening over a region near the tip. These studies completely ignored the fact that dilation is associated with plastic yielding. Papanastasiou and Thiercelin (1993) and Papanastasiou (1997a) carried out an elastoplastic analysis fully coupled with fluid flow and a robust fracture propagation criterion. They showed that the propagating elastoplastic fractures are wider than the elastic fractures. The wider openings of the elastoplastic fractures result in smaller fluid-lag regions than the elastic fractures. This result was in a complete contradiction with the dilation hypothesis. The second approach says that the value of the fracture toughness measured in the laboratory underestimates the in-situ value (Shlyapobersky, 1985). Estimates based on hydraulic fracturing field data have shown that typical toughness values are one to two orders of magnitude higher than those determined from conventional laboratory tests (Shlyapobersky, 1985). The high value of the in-situ fracture toughness could be attributed to a scale effect, to the influence of confining in-situ stress, to micro-cracking and to plasticity in the process zone near the crack-tip. In an attempt to explain the observed high net-pressures, Valco and Economides (1994) introduced a continuum-damage model. The damage introduced was a function of the fracture length. Experimental results show however that, after some crack propagation, the damage or plastic zones near the crack-tip develops fully; i.e. becomes independent of fracture length. Furthermore, they allowed sub-critical crack growth, i.e the fracture can propagate at a value of stress intensity less than the fracture toughness, in a time-dependent damage model based on a stress-corrosion cracking. Stress corrosion cracking is however a very slow process compared to the propagation of hydraulic fractures and therefore must not affect the response. Finally, their model was not fully coupled with the fracture propagation model to properly take into account the stress field ahead of the crack tip. The first results of work on the influence of inelastic rock behaviour in hydraulic fracturing were presented in references by Papanastasiou and Thiercelin (1993) and Papanastasiou (1997a). It was shown that plastic yielding results in a shorter and wider fracture than the elastic fracture and that wider openings of the elasto-plastic fractures result in smaller fluidlag regions than the elastic fractures. The pressure from the smaller fluid-lag region partially compensates the dissipated energy in the plastic zones. The difference strongly depends on the size of the plastic zones. The size of the plastic zones increases with the contrast of the in-situ stresses and is strongly influenced by the strength of the rock, the effective Young’s modulus and the pumping parameters (fluid viscosity and propagation velocity). This paper focusses on evaluating the effective fracture toughness (EFT) approach using a fully deterministic model based on finite element analysis. The fully coupled fluid-flow,
The effective fracture toughness in hydraulic fracturing 129
Figure 1. Geometry for a plane strain hydraulic fracture.
rock deformation and fracturing process model is summarized in the next section. The EFT determination is presented in Section 3. Section 4 presents three examples of elastoplastic fracture propagation under different conditions, with determination of the EFT using the path independent J -integral. Furthermore, results of the elastic model with EFT are compared with results of the elastoplastic model. Section 5 presents the results in tabular form and discusses the dependence of the EFT on the more important physical parameters. The main results and conclusions are outlined in Section 6. 2. Elasto-plastic hydraulic fracturing model This study examines the importance of plasticity and of the resulting EFT in the near tip processes. The classical plane strain fracture geometry (Figure 1) is analyzed because it is the simplest way to appropriately describe the behaviour near the crack tip. Although a meshing/re-meshing scheme is employed in order to carry out longer propagations, it is not claimed that the analysis is conducted under conditions representative of a massive fracture treatment. For example, when the fracture is properly contained, other fracture geometries are more appropriate than the plane strain geometry. However, it is suggested that the results from this analysis can be reasonably useful to other fractured geometries since the area near the fracture tip can be modelled locally as a plane strain geometry (Figure 2). The main physical processes which govern hydraulic fracturing in a weak formation are: • viscous fluid flow in the fracture, • elasto-plastic deformation caused by the stress concentration due to in-situ stresses and the action of fluid pressure, and • fracture propagation into the rock formation. Fluid leak-off from the fracture into the formation is not taken into consideration here, assuming that it is minimized by fluid additives forming a filtercake on the fracture walls. Leak-off from the fracture into the formation and pore pressure diffusion may change the effective stresses in the formation and hence the size of the plastic zones. However, the emphasis in this
130 Panos Papanastasiou
Figure 2. Geometry for a vertical hydraulic fracture.
study is put on determining the increase of EFT due to inelastic rock deformation. The effect of leak-off is left to a future study. Results on the influence of fluid leak off and pore pressure in a hydraulic fracture in a poro-elastic medium can be found in reference by Detournay and Cheng (1991). 2.1. F LUID
FLOW
It is assumed that the fracturing fluid is an incompressible, uniform, Newtonian viscous fluid. The continuity equation which imposes conservation of mass in one-dimensional flow, is ∂w ∂q + = 0, ∂t ∂x
(1)
where w is the local fracture width and q is the local flow rate. It is assumed that motion in the fracture occurs only along the x-direction. The first term in (1) represents the storage capacity of the fracture that arises from changes in the fracture width and the second term represents the mass flow in a cross-section perpendicular to the fracture axis. Equation (1) ignores any leak-off from the fracture surface into the rock formation. A second equation is derived from conservation of momentum. The lubrication equation, which relates the pressure gradient to the fracture width for a Newtonian fluid of viscosity µ degenerates to q=−
w 3 ∂p , 12µ ∂x
(2)
where p is the fluid pressure in the fracture. The fluid-flow (1) and (2) were supplemented with the boundary conditions of constant flow rate at the wellbore and zero fluid-pressure at the fluid front q(0) = q0 , p(l) = 0.
The effective fracture toughness in hydraulic fracturing 131 2.2. I NELASTIC
ROCK DEFORMATION
In soft rocks, such as clay or poorly consolidated sandstones, inelastic deformation, whose extent depends on rock properties and loading stresses, is expected to take place in the area near the crack tip due to excessive stress concentration. When the inelastic zone is small compared to the fracture length, solutions from linear elastic fracture mechanics can be used to analyze the fracture process. However, in weak rocks it is less likely that the inelastic zone will be small enough to justify the use of solutions based on linear elasticity. In such a case one should use plasticity theory which can describe properly the irreversible deformation due to excessive shear stresses around the fracture tip. In a stationary crack, plastic yielding takes place in the area near the crack tip, which is characterized by high stress concentration. When the fracture propagates the plastic zones unload elastically behind the advancing crack and the new area near the current tip deforms plastically Papanastasiou (1997a). In summary, the rock mass remote from the fracture may deform elastically, whereas the area near the body of the fracture is initially elastic then deforms plastically and finally unloads elastically after the fracture has advanced. Under such conditions, the plasticity model must be capable of dealing with nonproportional loading. Such capabilities are provided by an incremental flow theory of plasticity and finite element analysis. Among the different yield criteria, the Mohr–Coulomb model adequately describes the pressure-sensitive behaviour of rocks which exhibit dilatancy (increase in porosity) when sheared. The Mohr–Coulomb model is usually employed in cases of compressive stresses. Unlike most cases of classical fracture mechanics, the remote stress field in the hydraulic fracturing problem is compressive, due to the presence of the in-situ stresses. In such a case the use of Mohr–Coulomb yield criterion augmented with a tensile cut-off is justified. The tensile failure along the propagation line is modelled in the next section by a cohesive-type model. The Mohr–Coulomb criterion can be written as σe = σ1
1 + sin φ − σ3 1 − sin φ
(3)
σe is the equivalent stress, φ is the angle of internal friction and σ1 and σ3 are the maximum and minimum principal stresses. The value of equivalent stress under zero confining pressure is identical to the uniaxial compressive strength, σc , and it is related to the material cohesion c via σc = 2c
cos φ . 1 − sin φ
(4)
In the flow theory of plasticity the strain increment dij is decomposed into an elastic dije p and a plastic part dij p
dij = dije + dij .
(5)
The elastic strain increment dije can be obtained from Hooke’s law. The plastic strain p increments, dij , are generated when the yield surface σe in (3) is reached and can be generally expressed by a nonassociated flow rule in the form p
dij = dλ
∂Q , ∂σij
(6)
132 Panos Papanastasiou
Figure 3. Representation of the fracturing process of rock.
where Q is the plastic potential and the scalar function dλ is the so called plastic multiplier. The plastic potential Q can have a similar form to the yield surface σe if the dilation angle ψ replaces the friction angle φ. This allows nonassociated plasticity to be dealt with, and associated rules to be obtained as a special case, by making ψ = φ or Q = σe . The dilation angle is defined by the following equation sin ψ =
dvp , dγ p
(7)
where dvp is the increment of the plastic volumetric strain and dγ p is the increment of the plastic shear strain. In general, weak rocks obey a nonlinear yield criterion and exhibit nonassociative behaviour. Experimental results from triaxial compression tests show (a) that the dilation angle increases slightly with increasing plastic strain when low confining pressures are used but remains approximately constant in samples under high confining pressure and (b) the value of dilation angle is a strong function of the confining pressure. At tests of low confining pressure the measured dilation angle is greater than the friction angle but decreases rapidly with increasing confining pressure, eventually becoming negative. This indicates compaction at high confining pressure. In the problem of hydraulic fracturing, the initial in-situ mean pressure near the crack tip decreases during propagation. Under such conditions it is reasonable to assume an associative behaviour. Furthermore, earlier computations by Papanastasiou and Thiercelin (1993) showed that the nonassociative solution (e.g. zero dilation) was bounded by the associative solution and the elastic solution. 2.3. F RACTURE
PROPAGATION CRITERION
The most robust propagation criterion in nonlinear mechanics is based on constitutive modelling of the cohesive zone. Cohesive zone is the region ahead of crack tip which is characterized by micro-cracking around the crack tip and interlocking along a portion of the crack (Labuz et al., 1985; Bazant, 1986). For modelling the fracture process of rock, an effective crack length is assumed which consists of a length of true crack and a cohesive zone which is the nonlinear region where a closing stress acts (Figure 3) (Barenblatt, 1962). The closing stress is determined by the softening stress-displacement relation that various rocks exhibit in the post-peak region during a displacement control uniaxial tensile test. This model implies that normal stress continues to be transferred across a discontinuity which may or may not be visible. It is assumed that this stress transfer is due to grain bridging and aggregate interlocking of the three-dimensionality of the crack surfaces. In this separation process, the restraining
The effective fracture toughness in hydraulic fracturing 133
Figure 4. Constitutive model for the cohesive zone.
stress σ is a function of opening displacement δ, which falls off to zero at a critical opening displacement (δ = δc ) (Figure 4). It is also assumed that the cohesive zone localizes, due to the softening behaviour, into a narrow band ahead of the true crack tip. The last assumption is very convenient for finite element analysis, where the softening behaviour is modelled by one-dimensional interface elements lying in the direction of crack growth, ahead of the true crack tip (Hillerborg et al., 1976). The direction of crack growth is assumed to be known from the symmetry conditions. As already mentioned, the constitutive behaviour of interface elements can be defined by the stress-displacement curve which is usually derived from a displacement control uniaxial tensile test. Because of the confinement by the in-situ stresses in the problem of hydraulic fracturing, the stress path ahead of the fracture tip has more similarities with the triaxial extension test than with the direct tensile test. Furthermore, a direct tensile test is difficult to perform in weak or poorly consolidated sandstones. Due to the lack of available data, and in order to keep the number of input parameters minimum, a linear softening material was assumed (Figure 4). The area under the stress-displacement curve, σ (δ) equals the strain energy release rate GIC when the size of the cohesive zone is small compared to the crack length. For elastic solids the strain energy release rate, GIC is related to rock fracture toughness, KIC via: 2 KIC =
GIC E , 1 − ν2
(8)
where E is the Young’s modulus and ν is the Poisson ratio. The fracture toughness, KIC , is defined as a material parameter for elastic solids and can be calculated in a lab experiment as long as no plastic yielding takes place around the tip. If one assumes reasonable values for the rock fracture toughness KIC and the uniaxial tensile strength of the rock (approximately 8 to 10 times less than the uniaxial compressive strength, σc = 10σT ) the stress-displacement relation, σ (δ), is uniquely determined from equation: σ = σT (1 − δ/δc ),
(9)
where δc =
2 (1 − ν 2 ) 2 KIC . E σT
(10)
134 Panos Papanastasiou 2.4. I NITIAL
SOLUTION
The numerical model requires initial conditions to be specified at t = 0 for the initial fracture length `(0), width profile w(x, 0) and fluid pressure p(x, 0). In setting these initial conditions, the analytic solution derived by Desroches et al. (1994) for an elastic material with zero fracture toughness has been used. This solution is summarized as follows: # 0 1/3 " √ 6 2 1 E µv 1 p(x, 0) = σh min + , (11) − 3 π `(0)1/3 (`(0) − x)1/3 w(x, 0) = (`(0) − x)
2/3
µv 1/3
"
E0
x 7.21 − 3.17 1 − `(0)
5/6 # ,
(12)
where σh min is the in-situ stress perpendicular to the fracture axis, E 0 = E/(1 − v)2 is the plane strain modulus, µ is the fluid viscosity and v is the propagation velocity given in terms of flow rate q(0) by v = 0.835q(0)
3/4
`(0)
−1/2
E0 µ
1/4 .
(13)
The solution described by equations (11)–(13) has been derived under the assumption of constant velocity v = const and it pertains to the near tip region. This solution has recently been extended by Carbonell and Detournay (1998) and Detournay (1998) to apply to the whole of the hydraulic fracture. Finally, elastoplastic fracture model requires an initial fracture length. It has been shown by Papanastasiou (1997a), where propagations were carried out with different initial lengths, the initial fracture length had an effect on the fracture profiles in the initial fracture length interval, but, the fracture profiles were correctly calculated in the region where the fractures were propagated. 2.5. N UMERICAL
IMPLEMENTATION
The described time-dependent nonlinear hydraulic fracturing problem, which is characterized by fluid-flow, rock deformation and fracturing processes, is solved coupled by a combined finite difference-finite element scheme. The fluid-flow in the fracture is descritized by finite differences, the rock deformation by finite elements and the fracturing process by interface elements. A meshing/remeshing scheme is employed in order to permit efficient computations with fine mesh near the fracture tip during propagation. A special continuation method based on the volume of injected fluid in the fracture was used for direct coupling of the fuid-flow with rock deformation and for controlling the solution during propagation. 3. Effective fracture toughness A crack will propagate in an elastic solid when the stress intensity factor KI reaches a critical value of fracture toughness KIC : KI = KIC .
(14)
The effective fracture toughness in hydraulic fracturing 135 The stress intensity factor is related to the strain energy release rate, GI : KI2 =
GI E . 1 − ν2
(15)
For linear elastic situations the strain energy release rate equals the path-independent J integral introduced by Rice (1968): GI = J,
(16)
where
Z ∂ui J = U dy − ti ds , ∂x 0
(17)
in which U is the strain energy density, ti is the traction vector, ui is the displacement vector and ds is an element of the arc along the integration contour, 0. The direction of propagation is along the coordinate axis x. The J -integral is independent of the actual path chosen, provided that the initial and end points of the contour, 0, are on opposite faces of the crack and the contour contains the crack tip. The J -integral in elasto-plastic solids was originally employed as a measure of the intensity of stress and deformation fields near the tip of a stationary crack under monotonic loading. It can be used for characterizing crack initiation and limited amounts of crack extension (Kanninen and Popelar, 1985). In the problem of hydraulic fracturing due to significant elastic unloading which takes place behind the propagating crack, the so-called J -resistance curve is no longer a unique material property, but becomes a function of the flawed geometry. Under such conditions an alternate fracture criterion, independent of the geometry and flawed configuration must be utilized. Such a criterion is the cohesive model which was already described and used in this study. The J -integral is used here not as propagation criterion but to estimate the energy dissipated during fracture propagation and determine an EFT from (8) to (17) r JE KEFT = . (18) 1 − ν2 Using this value in an elastic analysis gives a good comparison with propagation pressures and fracture shapes obtained from the elastoplastic analysis. These comparisons will be shown for three examples in the next section. For elasto-plastic material the strain energy density U in (17) must be augmented by the plastic component: U = U e + U p,
(19)
where the elastic part U e is given by U e = 12 σij ije ,
(20)
and the value of U p is calculated by accumulating the incremental plastic work σe d p over the strain path Z p p U = σe d p , (21) 0
136 Panos Papanastasiou
Figure 5. Width profiles for elastic (solid-lines) and elasto-plastic (dashed-lines) fractures.
Figure 6. Effective fracture toughness vs fracture growth.
in which σe and p are the equivalent stress (3) and effective plastic strain, respectively. For elastic fractures the J -integral value is independent of crack length, whereas for the elasto-plastic fractures it increases initially and eventually reaches an asymptotic value.
4. Computational results Three examples with computational results of elasto-plastic fracture propagation and EFT determination, under different conditions will be presented in this section. The results of fracture propagation were analyzed extensively in Papanastasiou (1997a). In this study the limit value EFT will be calculated from the propagation of the elastoplastic fractures. This value will then be implemented in a pure elastic model and the results, in terms of fracture width and propagation pressure, will be compared with the results of the original elastoplastic model. The orientation of the crack tip with respect to the in-situ stresses is such that, in the first two examples plastic yielding takes place in the plane of deformation whereas in the third example plastic yielding takes place out of deformation plane. The reason being that in a deep rock formation where the vertical stress is usually the maximum, the conditions of the first two examples are met near the fracture front which propagates vertically, whereas the conditions of the third example are encountered near the fracture front which propagates horizontally (Figure 2).
The effective fracture toughness in hydraulic fracturing 137
Figure 7. Comparison of width profiles – example 1.
4.1. E XAMPLE 1 The parameters used in examples 1 and 2 are given in Table 1. Figure 5 shows the profiles of propagating elastic (solid-lines) and elasto-plastic (dashed-lines) fractures for the same fracture lengths. The fractures were propagated from an initial length of 0.2 m to a final length of 1.80 m with propagation steps of 0.01 m. In Figure 5 only four steps were drawn. The effective jump in the plastic solution of Figure 5 is due to the fact that the initial fracture length was 0.2 m. As mentioned above the fractures were propagated using a cohesive-type criterion, and the path independent J -integral was calculated to estimate the increase in the strain energy release rate due to plastic yielding. It was found that during fracture propagation the size of the plastic zones increases initially, but eventually reaches a constant size. The value of the EFT is directly related to the size of the plastic zones. After some propagation steps, the plastic zones fully develop and the EFT reaches an asymptotic value. This is shown in Figure 6 where the EFT is plotted as a function of the fracture growth. The EFT is calculated from the asymptotic value. In√the present example it was found √ that the fracture toughness increased from KIC = 2 MPa m to a value of KEFT = 7 MPa m. The initial value of EFT equals the rock fracture toughness KIC and dictates the energy required for propagating the elastic fractures. The energy required for propagating the long elasto-plastic fractures is given by the limit value. Higher energy is needed for propagating the elasto-plastic fractures because plastic yielding softens the material surrounding the tip. This reduces the level of stress in the direction of propagation, providing an effective shielding. √ The next natural step is to perform an elastic analysis with a KIC value of 7 MPa m to estimate whether plasticity effect in hydraulic fracturing can be embedded in the notion of EFT. The results of the EFT model are compared with the results of plasticity in Figure 7 in terms of fracture profiles, and in Figure 8 in terms of net-pressures. It is obvious that the EFT elastic model matches the results of the plasticity model quite closely. 4.2. E XAMPLE 2 The size of the plastic zone increases with increasing differences in the magnitude of the in-situ stresses and is strongly influenced by the strength of the rock, the elastic modulus and the pumping parameters (fluid viscosity and flow rate) (Papanastasiou and Thiercelin, 1993). For this reason in the next computation series, the fluid viscosity (µ) and flow rate at the wellbore (q0 ) were increased (Table 1), in order to trigger larger scale plastic yielding
138 Panos Papanastasiou Table 1. Input parameters of elastoplastic hydraulic fracturing model. Elastic constants Young’s modulus Poisson ratio
E = 25 GPa ν = 0.2
Plastic constants Friction angle Dilation angle Uniaxial compressive strength
φ = 30◦ ψ = 30◦ σc = 20 MPa
Fracturing parameters Fracture toughness Uniaxial tensile strength
KIC = 2 MPa σT = 3 MPa
In-situ stresses Maximum in-situ stress Minimum in-situ stress
σv max = 50 MPa σh min = 25 MPa
Example 1 Fluid viscosity Flow rate
µ = 5 10−8 MPa sec = 50 cp q0 = 0.0001 m3/sec m
Example 2 Fluid viscosity Flow rate
µ = 5 10−6 MPa sec = 5000 cp q0 = 0.001 m3/sec m
√ m
Figure 8. Comparison of net-pressure profiles – example 1.
and larger values of EFT. Figure 9 shows the profiles obtained for elastic (solid-lines) and elasto-plastic (dashed-lines) fractures, for the same fracture lengths. Figure 10 shows the EFT, calculated from the J -integral, as a function of the fracture growth. For the case √ of large scale yielding, the calculated EFT reaches an asymptotic value of KEFT = 34 MPa m which is more than an order of magnitude higher than the assumed fracture toughness. The elastic, elasto-plastic fracture profiles and the profile based on the EFT are compared in Figure 11, at the last propagation step. Better comparison between the plasticity model and the EFT model was obtained in term of the net-pressures (Figure 12). It is expected that the fracture profiles
The effective fracture toughness in hydraulic fracturing 139
Figure 9. Width profiles for elastic (solid-lines) and elasto-plastic (dashed-lines) fractures.
Figure 10. Effective fracture toughness vs fracture growth.
of the EFT model would more closely matching the plasticity model if the material was less dilatant. It is clear that the EFT model gives the best results in the case of a small scale yielding. As mentioned above, the EFT values are calculated from the J -integral which assumes the existence of K-dominance region. ‘K-dominant’ is the region around the crack-tip where the stress and displacement fields, which are function of the applied load σ and crack length α, can be expressed through the stress intensity factor K = K(σ, α). K-dominance exists only when the inelastic region can be contained within an annular region surrounding the crack tip (Kanninen and Popelar, 1985). Obviously, for a large scale yielding this conditions may not be satisfied. Nevertheless, the EFT elastic model matches the results of plasticity quite satisfactorily. 4.3. E XAMPLE 3 Example 3 considers the more realistic case of a vertical fracture which propagates horizontally in a deep rock formation. In such a case the maximum principal stress is usually the vertical in-situ stress. Plastic yielding near the crack-front is expected to take place out of the deformation plane (Figure 2). The parameters used in this example are real field data which were obtained from specialized coring, logging and stress measurements (Table 2). The fracturing zone was in a poorly consolidated sandstone located at 1500 m. Furthermore, a linear hardening plasticity model was assumed, in contrast to the first two examples where an elastic-perfectly plastic material was used. In this model, the equivalent stress σe is assumed to
140 Panos Papanastasiou
Figure 11. Comparison of width profiles – example 2.
Figure 12. Comparison of net-pressure profiles – example 2.
be a function of the accumulated effective plastic strain d p . The linear hardening modulus, h, was derived from the loading and unloading moduli of a uniaxial compression test according to the equation (Papanastasiou, 1997b): h=
Eloading dσe = . p d 1 − Eloading/Eunloading
(22)
Figure 13 shows the profiles of propagating fractures which were obtained using the elastoplastic linear hardening model. For clarity in Figure 13 the fracture profiles were plotted only once every eight propagation steps. The effect of the initial fracture length (0.5 m) on the propagating elasto-plastic fractures is obvious. Figure 14 shows the value of EFT as a function of the fracture growth. The energy required for propagating the long elasto-plastic fractures is given by the √ limit value. In this example one should use √ a value for the fracture toughness of 8.7 MPa m instead of the original value of 1 MPa m. Because of the wide fracture the calculated fluid pressure along the fracture is almost constant. The fluid front reaches the cohesive zone but does not penetrate it, since it is assumed to be impermeable. The net-pressure needed for propagating the elastoplastic fracture is shown in Figure 15. The net pressure increases during the first propagation steps, during which the plastic zones are growing but will start to decrease after the plastic zones have fully developed. In contrast, the net pressure of a propagating plane strain elastic fracture √ declines continuously. The predicted net-pressure by the elastic model with EFT of 8.7 MPa m is also plotted in Figure 15 as a function of the fracture length. This net-pressure is calculated using the well known equation of a uniformly loaded crack of length α
The effective fracture toughness in hydraulic fracturing 141 Table 2. Input parameters of example 3. Elastic constants Unloading modulus Poisson ratio
Eunloading = 16.2 GPa ν = 0.3
Plastic constants Friction angle Dilation angle Initial uniaxial compressive strength Loading modulus
φ = 28◦ ψ = 28◦ σc = 4 MPa Eloading = 1.8 GPa
Fracturing parameters Fracture toughness Uniaxial tensile strength
KIC = 1.0 MPa σT = 0.5 MPa
In-situ effective stresses Vertical stress Minimum horizontal stress Maximum horizontal stress
σv =14 MPa σh min = 3.7 MPa σh max = 9 MPa
Pumping parameters Fluid viscosity Flow rate
µ = 10−7 MPa sec = 100 cp q0 = 0.0005 m3/sec m
√
m
Figure 13. Width profiles of propagating elasto-plastic fractures.
√ Pnet = KEFT / π α.
(23)
The comparison of fracture profiles obtained by the plasticity model and the EFT model is shown in Figure 16 for a fracture length of 8 m. The fracture profile for a uniformly loaded fracture is given by w(x) = 4
p 1 − ν2 Pnet α 2 − x 2 . E
(24)
142 Panos Papanastasiou
Figure 14. Effective fracture toughness vs fracture growth.
Figure 15. Comparison of net-pressures vs fracture length – example 3.
It is obvious that the net-pressure obtained by the EFT elastic model is reasonably close to that of the plasticity model. The deviation observed can be attributed to the influence of the cohesive zone which enforces a fluid lag, which is not taken into account in (23) and (24). Furthermore, the plasticity model predicts a narrower opening than the elastic EFT model, especially near the crack-tip. This is due to the influence of rock dilation. It is expected that the difference will diminish with increasing fracture length. 5. Parametric studies It is clear that the increase in the EFT depends on the scale of plastic deformation taking place. The complexity of the problem, due to the material model and loading conditions, prevents the exact determination of the plastic zone boundaries and the resulting increase in the EFT. Instead, a scaling of the size of the plastic zones on the most important parameters will be attempted, assuming that the scaling of EFT is proportional. Furthermore, tables with values of EFT for a representative set of field data will be presented. The time dependence from the problem is eliminated here by considering a steady-state view. It is assumed that on a moving frame which takes the instantaneous velocity of the crack tip, the near-tip behaviour is slowly evolving (Desroches et al., 1994). The approximation of constant velocity is more appropriate near the fracture tip than the approximation of constant flow-rate in every cross-section. The latter may be more appropriate for the area near the
The effective fracture toughness in hydraulic fracturing 143
Figure 16. Comparison of width profiles – example 3.
wellbore (Geertsma and DeKlerk, 1969). However, the area of interest here is the region near the fracture tip. With the constant velocity and no fluid loss assumption, the continuity (1) can be written as q = v w,
(25)
where v is the propagation velocity. Equation (25) combined with the lubrication (2) provides the pressure gradient dp 12 µ v . =− dx w2
(26)
Finally, the steady state fluid-flow (26) was supplemented with the boundary conditions of constant velocity and zero fluid-pressure at the fluid front v = const.
(27)
p(l) = 0.
(28)
We emphasize here that the analysis in Section 4 was carried out with fluid-flow equations (1)–(3) and it pertains to the whole of the hydraulic fracture whereas the next computations are based on constant velocity assumption (25)–(28) and the analysis pertains to the near tip region. 5.1. S CALING
OF EFFECTIVE FRACTURE TOUGHNESS
Papanastasiou (1997a) has shown that the size of the plastic zone is a function of the deviator of the in-situ stresses. Small-scale yielding is expected if the stress field is nearly hydrostatic. A highly nonhydrostatic stress field is encountered in deep rock formations due to the significant difference between vertical and horizontal effective in-situ stresses. In addition, Anagnostou and Kovari (1993) suggested that in the structure of the solution of elastoplastic problems dealing with Mohr–Coulomb material the following term appears σ0 =
c +p tan φ
(29)
144 Panos Papanastasiou where p = −(σ1 + σ3 )/2. Hence the dependence of the shape of plastic zone on the deviation of remote stresses from uniformity will be through the ratio σ¯ =
τ , σ0
(30)
where τ = (σ1 − σ3 )/2. Accordingly, a typical length, `σ , for the size of plastic zones depends on upon the remote stress field and rock parameters in the ratio σ¯ such that `σ increases when σ¯ increases. The higher the fracture toughness, the more difficult it will be for the fracture to propagate and the greater the size of the plastic zone that develops around the tip. For example, for an elastic-perfectly plastic material with a Mises yield surface and a tensile yield stress σT , the plastic zone size is proportional to K 2 . For a plane strain problem the approximate distance to the boundary ahead of the crack tip, `p , is 1 `p = 3π
K σT
2 .
(31)
Equation (31) is valid for the classical problem of linear fracture mechanics for an unloaded or uniformly loaded fracture. It has been shown, in an intermediate asymptotic analysis of the brittle hydraulic fracturing problem with zero fracture toughness (Desroches et al., 1994), that the tip region stress is dominated by a singularity which is weaker than the inverse square root singularity of linear elastic fracture mechanics. In the case of a Newtonian fracture fluid with viscosity µ, the stress σyy can be expressed as σyy = C +
2E 02/3 (µv)1/3 −1/3 x , 31/6
(32)
where the constant C is related to the wellbore pressure. E 0 is the plane strain modulus E 0 = E/(1 − ν 2 ). Along this line, a characteristic length of the elastoplastic problem is defined as `p ∼
E 02 µv . σc3
(33)
This is analogous to the argument leading to expression (31) from which an elastic analysis with an elastic stress intensity fator K leads to the result (31) for the elastoplastic zone length for an elastic-plastic (small scale) yielding solution. According to (32) and (33), high values of effective elastic modulus E 0 , viscosity µ and propagation velocity v; and low values of rock strength σe , lead to higher stress concentrations near the fracture tip. Consequently, larger plastic zones are expected to develop. The strong dependence of EFT on the pumping parameters was demonstrated in the presented examples 1 and 2 where only the viscosity µ and flow rate q0 had been varied. It is noted that (32) and (33) are valid only in the region near the fracture-tip, where the singular behaviour dominates. Far away from the fracture-tip, extra terms must be added to the leading terms. Strong correlations between elevations in EFT and fluid rheology have been documented experimentally by Holder et al. (1993). Although, the plane strain modulus E 0 normally determines the deformation rather than the stresses, its strong appearance in (32) and (33) can be justified in the following simple way; the elastic
The effective fracture toughness in hydraulic fracturing 145 Table 3. Variation of EFT with stress field and rock strength. plane strain modulus
E 0 = 31.25 GP a
friction and dilation angles rock fracture toughness
φ = ψ = 30o √ KIC = 2MP a m
pumping parameters
µv = 10−8 ∼ 10−7 MP a m rock strength σσTc
stress field σσ3
60 6
20 6
20 2
30 = 1.0 30
2.0
2.0
2.0
45 = 1.5 30
2.0
4.60
7.31
60 = 2.0 30
2.0
7.03
15.48
1
modulus E 0 influences the fracture width which in turn determines the load distribution on the fracture from the profile of the fluid pressure. 5.2. TABLES
OF EFFECTIVE FRACTURE TOUGHNESS
As mentioned above, the EFT is estimated once the plastic zones are fully developed. Therefore, the EFT is independent of the fracture length. Tables with values of EFT were derived for a representative set of field data. The input parameters on which the EFT depends are: the Young’s modulus E; Poisson ratio ν; friction angle φ; dilation angle ψ; rock fracture toughness KIC (laboratory measured value); the ratio of uniaxial compressive strength over the uniaxial tensile strength σc /σT and the ratio of maximum to minimum in-situ stress σ1 /σ3 . It was assumed that the cohesive zone, which is related to the rock fracture toughness and uniaxial tensile strength σT , is fully impermeable. It is also more convenient to derive the EFT, as a function of the product of fluid viscosity times propagation velocity µv, using the fluid-flow equation (26). The main conclusion from Table 3 is that in a strong rock formation (σc /σT = 60/6), plastic yielding does not take place even in a highly nonhydrostatic stress-field. In weak formations (σc /σT = 20/2; Table 4) large scale plastic yielding occurs under a highly nonhydrostatic field and high values of pumping parameters resulting in significant increase in the EFT. On comparing the last column of Table 3 with the 2nd and 3rd columns of Table 4, where the only parameter changed is the rock fracture toughness KIC , it is obvious that, when the rock fracture toughness increases the fracture propagates with more difficulty resulting in larger plastic zones and significant increase in EFT value. Finally, it is shown that plastic yielding does not take place under a hydrostatic field, even in weak formations. However, it is expected that the near crack-tip region is under a nonhydrostatic stress field since in most of cases (deep rock formations) the vertical geostatic stress is higher than the horizontal stresses. The last remark is in favour of fracture containment especially if the barriers are clay layers which plastified more readily than sandstones.
146 Panos Papanastasiou Table 4. Variation of EFT with stress field and pumping parameters. plane strain modulus
E 0 = 31.25 GP a
friction and dilation angles rock fracture toughness
φ = ψ = 30o √ KIC = 1MP a m
rock strength
σc 20 σT = 2
pumping parameters µ v MP a m stress field σσ3
10−8
10−7
10−6
30 = 1.0 30
1.0
1.0
1.0
45 = 1.5 30
1.0
1.88
5.66
60 = 2.0 30
3.87
5.07
13.04
1
6. Conclusions In the present study, a fully deterministic analysis of effective fracture toughness approach was presented. The analysis was based on a coupled elastoplastic hydraulic fracturing model and finite element analysis. The value of effective fracture toughness was calculated using the path independent J -integral. It has been shown that plastic yielding near the tip of a propagating fracture provides an effective shielding, resulting in an increase of the effective fracture toughness (EFT) by more than an order of magnitude. The value of the EFT is directly related to the size of the plastic zones. After some propagation steps the plastic zones fully develop and the EFT reaches an asymptotic value. The size of the plastic zones and the resulting EFT increase with the contrast of the magnitude of the in-situ stresses. They are also strongly influenced by the strength of the rock, the elastic modulus and the pumping parameters (fluid viscosity and flow rate). It has been demonstrated that an elastic model using the concept of EFT, matches the results of plasticity quite well. Finally, scaling and tables with values of EFT for a representative set of physical parameters were presented. Acknowledgement The author is grateful to Dr. M. Thiercelin and Prof. E. Detournay for useful discussions. References Anagnostou, G. and Kovari, K. (1993). Significant parameters in elastoplastic analysis of underground openings. Journal of Geotechnology Engineering ASCE 119(3), 401–419. Barenblatt, G.I. (1962). Mathematical theory of equilibrium cracks in brittle fracture. Advances in Applied Mechanics 7, 55–129. Barr, D.T. (1991). Leading-Edge Analysis for Correct Simulation of Interface Separation and Hydraulic Fracturing. Ph.D. thesis, MIT. Bazant, Z.P. (1986). Mechanics of distributed cracking. Applied Mechanics Review 39(5), 675–705.
The effective fracture toughness in hydraulic fracturing 147 Carbonell, R. and Detournay, E. (1998). Self-similar solution of a fluid-driven fracture. Proceedings of the Royal Society of London A (submitted). Cleary, M.P., Wright, C.A. and Wright, T.B. (1991). Experimental and modeling evidence for major changes in hydraulic fracturing design and field procedures. Proceedings SPE Gas Technology Symposium, Houston, 131–146. Desroches, J., Detournay, E., Lenoach, B., Papanastasiou, P., Pearson, J.R.A, Thiercelin, M. and Cheng, A.H.-D. (1994). The crack tip region in hydraulic fracturing. Proceedings of the Royal Society of London A 447, 39–48. Detournay, E. (1998). Fluid and solid singularities at the tip of a fluid-driven fracture. Proceedings IUTAM Symposium on Non-Linear Singularities in Deformation and Flow (Edited by D. Durban and J.R.A. Pearson), Dordrecht, Kluwer Academic Publishers. Detournay, E. and Cheng, A.H.-D. (1991). Plane strain analysis of a stationary hydraulic fracture in a poroelastic medium. International Journal of Solids and Structures 37, 1645–1662. Gardner, D.C. (1992). High fracturing pressures for shales and which tip effects may be responsible. Proceedings of the 67th Annual Technology Conference and Exhibition of the SPE, Washington DC, 879–893. Geertsma, J. and DeKlerk, F. (1969). A rapit method of predicting width and extend of hydraulic induced fractures. JPT 21(12), 1571–1581. Hillerborg, A., Modeer, M. and Petersson, P.E. (1976). Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research 6, 773–782. Holder, J., Morita, N., Kenrick, A.J., Thallak, S.G. and Gray, K.E. (1993). Measurements of effective fracture toughness values for hydraulic fracture. SPE Production Operations Symposium, Oklahoma C., SPE 25491. Kanninen, M.F. and Popelar, C.H. (1985). Advanced fracture mechanics, Oxford University Press, 386. Kristianovitch, S.A. and Zheltov, Y.P. (1955). Formation of vertical fractures by means of highly viscous fluid. Proceedings of the 4th World Petroleum Congress 2, 579. Labuz, J.F., Shah, S.P. and Dowding, C.H. (1985). Experimental analysis of crack propagation in granite. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 22(2), 85–98. Papanastasiou, P. and Thiercelin, M. (1993). Influence of inelastic rock behaviour in hydraulic fracturing. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 30(7), 1241–1247. Papanastasiou, P. (1997a). The influence of plasticity in hydraulic fracturing. International Journal of Fracture 84, 61–97. Papanastasiou, P. (1997b). A coupled elastoplastic hydraulic fracturing model. International Journal of Rock Mechanics and Mining Sciences 34(3–4), paper No. 240. Rice, J.R. (1968). A path-independent integral and the approximate analysis of strain concentration by notches and cracks. Journal of Applied Mechanics Transition ASME 35, 379–386. Shlyapobersky, J. (1985). Energy analysis of hydraulic fracturing. Proceedings of the 26th US Symposium on Rock Mechanics, Rapid City, 539–546. Spence, D.A. and Turcotte, D.L. (1985). Magma-driven propagation of cracks. Journal of Geophysics Research 90, 575–580. Valco, P. and Economides, M.J. (1994). Propagation of hydraulically induced fractures – a continuum damage mechanics approach. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 31(3), 221–229. Van Dam, D.B., de Pater, C.J. and Romijn, R. (1997). Experimental study of the impact of plastic rock deformation on hydraulic fracture geometry. International Journal of Rock Mechanics and Mining Sciences 34(3–4), paper No. 318.