An extended J-integral for evaluating fluid-driven cracks in hydraulic fracturing

1. Introduction

Parameters of stress intensity factor (SIF), crack tip opening displacement(CTOD), energy release rate and J-integral are most commonly used to evaluate crack geometry, stress field, stability state of fracture tip, etc. The usage and calculation of each of these parameters are often limited to certain conditions such as elastic materials, small-displacement assumption or small-scale plastic region. Among these parameters, J-integral, however, has been found to be the most effective method in characterizing fracture tip due to its wider application range.

J-integral is a path integral along the contour starting from any point on bottom surface of the crack and ending in top surface. It was initially proposed independently by Cherepanov (1967) and Rice (1968a) as an alternative approach to determine the energy release rate. The original form of J-integral is expressed as(1) $J = \int_{Γ} (w d y - T \frac{\partial u}{\partial x} d s)$ where w is the strain energy density, Γ is the integral path from the bottom to the top surface around the crack, T is the traction along the path, u is the displacement vector along the path, and ds is the line element along the path.

Rice (1968a) showed that the value of J-integral is equal to the energy release rate during crack extension, regardless of its integration path. The path-independence made the J-integral a powerful tool to investigate a wide range of studies, such as different types of loading, material laws and field problems, in both linear elastic and elastoplastic conditions, especially in connection with numerical analysis. One of the main advantages of J-integral is that it can be calculated in a region far enough from crack tip so that the numerical accuracy is not compromised. In case of small-scale yielding, in which stress field beyond yield zone is still governed by stress singularity, Rice (1968a) showed that the path-independent J-integral can be correlated to the SIF (KI) in mode-I fracture loading as follows:(2) $J = \frac{1 - ν^{2}}{E} K_{I}^{2}$ where ν is the Poisson's ratio, and E is the Young's modulus. In case of plane-stress condition, Eq. (2) applies by replacing 1−ν2 with unity. Similarly, Eq. (2)can be generalized for mode-II and mode-III loading cases.

In case of elastoplastic conditions, J-integral is also capable of characterizing the intensity of the crack tip stress field. Therefore, the critical value of J-integral (Jc) plays a comparable role in determining fracture toughness in elastoplastic condition as critical SIF does in elastic condition (Kuna, 2013). Landes and Begley (1977) proposed experimental procedures using a series of specimens to determine Jc. Subsequently, standard experimental procedure was developed (ASTM E1737-96; 1996). Either a critical value or a resistance curve can be obtained by using particular standardized specimen (Tada et al., 2000) for a test material. J measurement is also a substitute where the value of SIF (K) is not available, especially when the plastic zone around the crack tip is too large.

Numerous attempts have been made in the literature to extend the applications of the path-independent J-integral. Begley and Landes (1972) applied the J-integral to determining the fracture toughness in elastoplastic materials. Their results were validated against experimental determination of JIC (for mode-I) through different measuring points (Landes and Begley, 1974, Landes and Begley, 1977). In addition, Landes and Begley (1976) proposed a modification of J-integral under nonlinear viscous flow rule, denoted as C*, in which strain and displacement were replaced by their time derivatives, respectively. Schapery (1984) developed a generalized J-integral which is applicable to a wide range of viscoelastic materials. By such modification, the application of J-integral was extended to high-temperature creep cracking phenomena (Taira et al., 1979). For bi-material interface crack problems, M1 integral, another extended version of J-integral, was developed (Miyazaki et al., 1993). More relevant applications of J-integral have been continuously developed by different researchers (Roberti et al., 1984, Landis, 2004, Prawoto and Onn, 2012, Zimmerman and Jones, 2013, Ochensberger and Kolednik, 2014).

Apart from the application of the J-integral as fracture criterion, it also offers a potential characterization tool to study stress field in yield region. Stress field around crack tip under elastic assumption can be simply acquired from the SIFs KI, KII and KIII, which were given by Irwin (1957) using Westergaard (1939)'s method. In reality, however, yield or plastic zone exists around crack tip and therefore, elastoplastic fracture mechanics needs to be considered to take into account the yield region. In order to address such effects, Hutchinson (1968)and Rice and Rosengren (1968) applied power-law deformation theory of plasticity and demonstrated that the near-tip stress field can be expressed as a function of J-integral. This near-tip area is called HRR field (HRR for initials of the authors), or J-dominance zone, as shown in Fig. 1. It is noteworthy that the HRR field solution is applied under the assumptions of small strain, infinite plate, and monotonical as well as proportional loading since the deformation theory cannot be used when unloading occurs.

The numerical computation of J-integral, in particular with finite element method (FEM), has also been widely studied. J-integral can be calculated directly by interpolating the variables on predefined arbitrary path using numerical computation. Such an interpolative approach is, however, associated with uncertainty as stress magnitudes are calculated for Gauss points, whereas displacements are solved over nodal points. In order to address such an uncertainty, equivalent domain integral (EDI) (Hutchinson, 1968, Li et al., 1985, Nikishkov and Atluri, 1987, Raju and Shivakumar, 1990) has been proposed to compute the J-integral value. EDI is a generalization of virtual crack extension method (Parks, 1977, Parks, 1978, deLorenzi, 1982, Li et al., 1985), which essentially evolves from the calculation of energy release rate directly. The domain integration alleviates the local solution error.

The application of the J-integral to hydraulic fracturing needs further improvements as there are issues such as internal traction, porous media and, in some cases, thermal strain. It is noteworthy that hydraulic fracturing is a technique used to facilitate oil and gas extraction by creating cracks under induced fluid pressure and it is now an indispensable process in the development of tight gas and shale gas reservoirs. Karlsson and Bäcklund (1978)developed a revised J-integral for internally loaded cracks in two-dimensional (2D) cases. However, a more detailed extension for three-dimensional (3D) internally pressurized cracks and the cracks under thermal conditions are not available currently.

In this study, J-integral is further extended to make it suitable for investigating 3D hydraulically pressurized fractures and its path-independence is successfully examined. Conditions of porous media and thermal strains are discussed. The application of the extended J-integral is assessed and the results are compared against those obtained by the currently available methodologies. The main purpose of the proposed methodology is to apply the concept of J-integral to a much wider range of fracture treatments.

2. Formulation of 3D J-integral and its application to hydraulic fracturing

J-integral can be represented by the subtraction of its two components as follows:(3) $\begin{array}{l} J = J_{1} - J_{2} \\ J_{1} = \int_{Γ} w d y \\ J_{2} = \int_{Γ} T \frac{\partial u}{\partial x} d s \end{array}}$

Rice (1968a) indicated that J-integral values are equal along all paths around crack tip. This conclusion was derived based on conservative enclosed path, that is, for an enclosed path without singular point inside, J is proved to be zero. This statement was made based on the following assumptions: (1) strain-displacement relation is linear (small strain); (2) no body force exists; and (3) no singularity point exists along/within the integration line. A detailed derivation is shown in Appendix A.

Rice (1968a) suggested that J-integral is a measure of average strain energy of crack tip in the case of a blunt crack. In case of a sharp tip, the contour can shrink to as near as crack tip and yet not reaching the sharp tip:(4) $J = \lim_{r \to 0} \int_{Γ_{r}} (w d y - T \frac{\partial u}{\partial x} d s)$ where r is the distance to the crack tip.

Therefore, we could consider that J-integral is a characterizing parameter of the crack tip field. From an energy point of view, J-integral represents the potential energy release rate.

Under linear elastic condition, J is related to KI, as shown in Eq. (2). By working out the integration of J, it is found that if the integration route is circular and within K-dominance zone, J1 and J2 are also distance-independent of crack tip, which can be used as an alternative technique to calculate KI:(5) $J_{1} = \int_{c} w d y = r \int_{- π}^{π} w \cos θ d θ = \frac{K_{I}^{2} (1 + ν) (1 - 2 ν)}{4 E}$ (6) $J_{2} = \int_{c} T \frac{\partial u}{\partial x} d s = r \int_{- π}^{π} T (\frac{\partial u}{\partial r} \frac{\partial r}{\partial x} + \frac{\partial u}{\partial θ} \frac{\partial θ}{\partial x}) d θ = \frac{K_{I}^{2} (1 + ν) (- 3 + 2 ν)}{4 E}$ where c is the circular integration route.

The derivations of Eqs. (5), (6) are presented in Appendix B.

2.1. Crack pressurized with fluid

For internally pressurized crack (as shown in Fig. 2), the path-independence of the J-integral around the crack tip is not applicable as the traction-free condition at the crack faces is not met. Therefore, an addition term is required to maintain this path-independence.

Due to the traction force (hydraulic pressure) over the crack surfaces (see Fig. 2a), the corresponding components of the J-integral on routes Γ2 and Γ4 are no longer zero and they can be expressed as follows:(7) $\begin{array}{l} J_{Γ_{2}} = - \int_{Γ_{2}} T \frac{\partial u}{\partial x} d s = - \int_{Γ_{2}} p \frac{\partial u_{y, Γ_{2}}}{\partial x} d s \\ J_{Γ_{4}} = - \int_{Γ_{4}} T \frac{\partial u}{\partial x} d s = - \int_{Γ_{4}} (- p) \frac{\partial u_{y, Γ_{4}}}{\partial x} d s \end{array}}$ where p is the traction force (hydraulic pressure) over the crack surfaces.

Assuming that the crack is symmetric, and the end points of the route are also symmetric, it can be obtained that(8) $J_{Γ_{2}} + J_{Γ_{4}} = - 2 \int_{Γ_{2}} p \frac{\partial u_{y}}{\partial x} d s$

Assuming Γ3 to shrink to an infinitely small contour around the crack tip and yet excluding crack tip singularity (as shown in Fig. 2b), or become a rounded notch (as shown in Fig. 2c), the J-integral over Γ3 can be expressed as follows (for enclosed path integral conservation, refer to Appendix A):(9) $J_{Γ_{3}} = - (J_{Γ_{1}} + J_{Γ_{2}} + J_{Γ_{4}}) = - \int_{Γ_{1}} (w d y - T \frac{\partial u}{\partial x} d s) + 2 \int_{Γ_{2}} p \frac{\partial u_{y}}{\partial x} d s$

Since the path Γ3 is in clockwise direction (as shown in Fig. 2b and c), J value in the opposite direction Γ3∗ is reversed, i.e.(10) $J_{Γ_{3}^{*}} = \int_{Γ_{1}} (w d y - T \frac{\partial u}{\partial x} d s) - 2 \int_{Γ_{2}} p \frac{\partial u_{y}}{\partial x} d s$

The length of path Γ2 in Eq. (10) can be approximated as the distance from the J-contour/crack face intersection to the crack tip. The extended J-integral is then presented as(11) $J_{r} = J - \int_{0}^{r} p \frac{\partial δ}{\partial x} d r$ where J is the original value along the path, and δ represents the crack opening width (aperture). In fluid-driven cracks, before the cracks reach propagation criteria, fluid pressure builds up inside the cracks without flowing, therefore, the hydraulic traction force p is the same along the crack surfaces. With application of integration by parts, Eq. (11) can be written as follows:(12) $J_{r} = J - (p δ |_{- r}^{0} - \int_{- r}^{0} δ \frac{\partial p}{\partial r} d r) = J + p δ_{- r}$

This extended J-integral is used to characterize local strain energy in case of hydro-traction crack faces. Eq. (12) is applicable when the start and end points are symmetric. For deviating end points, however, one readily obtains(13) $J_{r} = J + p u_{y, 1} - p u_{y, 2}$ where uy,1 and uy,2 denote the y-displacements of upper and lower end points, respectively. In this study, for the sake of simplicity, it is assumed that all contour edge points are symmetrically distributed on both crack faces, unless otherwise mentioned.

The extended J value can also be derived from an energy point of view. In case of the blunt fracture tip, J-integral can be shrunk to the flat surfaced notch and therefore, it can be used as a parameter to characterize the local strain energy. From the energy point of view, it can be deemed as the potential energy release rate. The potential energy of the area around the crack tip is(14) $π = \int_{A} w d x d y - \int_{Γ} T u d l$

The first term at the right-hand side of Eq. (14) is the free energy of the body, while the second term is the energy coming from boundary loading. The body and the load together are treated as a composite system. If the crack propagates by a small amount da, the potential energy release rate is expressed as follows in terms of field variables:(15) $- \frac{\partial π}{\partial a} = \frac{\partial}{\partial a} \int_{A} (σ : \nabla u - w) d A$ where σ is the stress, and a is the distance that the crack propagates.

According to Rice (1968b), J-integral is equal to the energy release rate and can be expressed as(16) $J = - \frac{\partial π}{\partial a} = \frac{\partial}{\partial a} \int_{A} (σ : \nabla u - w) d A$

In case of internally pressurized cracks, the term $\int_{Γ}$ Tu dl in Eq. (14) includes the pressure work 2 $\int_{Γ}$ puydl. Therefore, the extended J-integral is also equal to the potential energy release rate of pressurized crack:(17) $J_{r} = - \frac{\partial π}{\partial a} = J + p δ_{- r}$

Since the energy interpretation of extended J-integral remains the same, and the domain integral method (Nikishkov and Atluri, 1987), commonly used for calculation of J, is essentially based on the energy release rate, the domain integral method can, therefore, still be used in hydraulically pressurized fractures. In the following section, the extended J-integral for hydraulically pressurized crack surfaces is generalized to 3D conditions.

2.2. Point-wise 3D extended J-integral

J-integral was initially proposed for 2D problems. It becomes surface integral if a crack is in 3D, which can be derived in the same manner because of volume integral conservatory as(18) $J^{3 D} = \int_{s} (w n_{x} d A - T \frac{\partial u}{\partial x} d A)$ where s represents the wrapping surface around the crack front segment, and nxis the x-component of the unit vector normal to the wrapping surface. A schematic representation of cylindrical wrapping surfaces around cracks with constant and changing heights is shown in Fig. 3a and b, respectively. Wrapping surfaces around the crack can be arbitrarily shaped, yet the surface geometry shown in Fig. 3 is simplified as a cylinder for better understanding. In case of traction-free crack surfaces, J value is obtained by integration over side surface (B) as well as over top and bottom surfaces (A and C), respectively. In case of internally pressurized cracks, another term, which represents the effect of the loaded crack surfaces M and N (see Fig. 3a and b), needs to be added to Eq. (18)as follows:(19) $J_{r}^{3 D} = \int_{s} (w n_{x} d A - T \frac{\partial u}{\partial x} d A) - 2 \int_{s} p \frac{\partial u_{y}}{\partial x} d A = J^{3 D} - \int_{Γ (M)} p \frac{\partial δ}{\partial x} d Γ$

The surface integral in Eq. (19), J3D, a global parameter for 3D cracks, represents the accumulation of strain energy along the crack front, and can be expressed as follows:(20) $J^{3 D} = \lim_{r \to 0} \int_{s} (w n_{x} d A - T \frac{\partial u}{\partial x} d A)$

Numerical implementation of the 3D J-integral is, however, associated with the difficulties due to the presence of singularities. For instance, integration over the top and bottom parts of the wrapping surface in Fig. 3a and the side cylindrical surface in Fig. 3b needs displacement gradient to be calculated near the vertex singularity. Therefore, the numerical solution of 3D J-integral is relatively inaccurate and unstable. However, this is not the case with point-wise 3D J-integral.

For a specific point on the crack front, point-wise J-integral is defined to characterize the state of a point, as shown in Fig. 4a. The point-wise J-integral can be calculated as follows (Nikishkov and Atluri, 1987):(21) $J_{point}^{3 D} Δ = \lim_{\begin{array}{l} r / Δ \to 0 \\ Δ \to 0 \end{array}} \int_{s} (w n_{x} d A - T \frac{\partial u}{\partial x} d A)$ where Δ is the disk thickness, and r is the radius of the disk. Due to the conservation of the volume integral, and the application of the divergence theorem, the point-wise J-integral can be represented as(22) $J_{point}^{3 D} Δ = \int_{A + B + C} (w n_{x} d S - T \frac{\partial u}{\partial x} d S)$

Since the normal vectors of faces A and C (top and bottom parts of the wrapping surface, as shown in Fig. 4a) have opposite directions, the point-wise J-integral can be expressed as(23) $J_{point}^{3 D} Δ = Δ \int_{A} \frac{\partial}{\partial z} (w n_{x} - T \frac{\partial u}{\partial x}) d A + Δ \int_{Γ_{B}} (w n_{x} - T \frac{\partial u}{\partial x}) d Γ_{B}$ where nx is zero in areas A and C.

Eq. (23) can be further simplified as(24) $J_{point}^{3 D} = \int_{Γ} (w n_{x} - T \frac{\partial u}{\partial x}) d Γ - \int_{A} (T \frac{\partial u}{\partial x}) d A$

The path Γ resides on a plane, of which the normal vector is tangent to the crack front, as shown in Fig. 4b.

Taking into account the hydraulic pressure, Eq. (22) is written as follows:(25) $J_{point}^{3 D} Δ = \int_{A + B + C + M + N} (w n_{x} d S - T \frac{\partial u}{\partial x} d S)$ which can further be deduced to(26) $J_{point}^{3 D} Δ = Δ \int_{A} \frac{\partial}{\partial z} (- T \frac{\partial u}{\partial x}) d A + Δ \int_{Γ_{B}} (w n_{x} - T \frac{\partial u}{\partial x}) d Γ_{B} + Δ \int_{Γ_{r}} (- p \frac{\partial δ}{\partial x}) d A = Δ [\int_{Γ} (w n_{x} - T \frac{\partial u}{\partial x}) d Γ - \int_{A} \frac{\partial}{\partial z} (T \frac{\partial u}{\partial x}) d A + p δ]$

Therefore, the extended point-wise J-integral for 3D hydraulically pressurized cracks can be expressed as(27) $J_{point}^{3 D} = \int_{Γ} (w n_{x} - T \frac{\partial u}{\partial x}) d Γ - \int_{A} \frac{\partial}{\partial z} (T \frac{\partial u}{\partial x}) d A + p δ$

The line integral term and pressure term can be calculated away from the crack front, yet the area integral term is still related to area near the crack front. From the second term on the right side of Eq. (27), it can be noted that taking the gradient over z-direction eliminates strain singularity, which makes the numerical calculation of the point-wise J-integral more stable. This local value of point-wise J-integral indicates the stress state of a point at the crack front, and it changes along the crack front.

In the following part, we will explain that the point-wise 3D extended J-integral is also the energy release rate at that point along a fluid-driven crack front. We take out a volume around a crack line with constant thickness Δ, as shown in Fig. 5a. The potential energy of the crack tip front volume is(28) $π = \int_{v} w d v - \int_{A} p δ d A - \int_{B} T u d B - \int_{c} T u d C$ where A, B and C are shown in Fig. 5. If the crack advances by an infinite small area δA, the change of potential energy is π(A+δA) − π(A). If we transform Fig. 5a so that the crack front coincides with that in Fig. 5b, the difference volume represents the difference in elastic energy, which is shown in Fig. 5c. The difference volume can be acquired by considering a sweeping motion between Fig. 5a and b. Therefore, we obtain the difference in elastic energy, i.e. $\int_{Γ} [- w n_{x} (δ A)] d Γ$ . Work done by the external traction also changes during sweeping. As to the work done by the traction on the circumferential area, the contributing difference is $\int_{Γ} T (\partial u / \partial x) (δ A) d Γ$ . For the top surface, consequently, the difference in the work done by the traction is $\int_{C} T (\partial u / \partial x) (δ a) d C$ , where $δ a$ is the increased crack length, $δ a = δ A / Δ$ . Similar expression can be derived for the bottom surface. Finally, the contribution of the internal pressure to the potential energy difference is $\int_{δ A} (- p) δ d (δ A)$ .

With all the contributions described above, the energy release rate is(29) $G = - \frac{π (A + δ A) - π (A)}{δ A} = - [\int_{Γ} (- w) n_{x} (δ A) d Γ + \int_{Γ} T \frac{\partial u}{\partial x} (δ A) d Γ \sqrt{b^{2} - 4 a c} + \int_{C} T \frac{\partial u}{\partial x} (δ a) d C + \int_{C^{'}} T \frac{\partial u}{\partial x} (δ a) d C^{'} + \int_{δ A} (- p) δ d (δ A)] / (δ A)$ where C ′ is the bottom surface.

If the thickness is infinite small, the summation of integrals for the top and bottom surfaces can be expressed as $\int_{C} \frac{\partial}{\partial z} (T \frac{\partial u}{\partial x}) (δ A) d C$ since the tractions are in opposite directions. Therefore, the energy release rate can be further simplified as(30) $G = - [\int_{Γ} (- w) n_{x} (δ A) d Γ + \int_{Γ} T \frac{\partial u}{\partial x} (δ A) d Γ + \int_{C} \frac{\partial}{\partial z} (T \frac{\partial u}{\partial x}) (δ A) d C + \int_{δ A} (- p) δ d (δ A)] / (δ A) = \int_{Γ} w n_{x} d Γ - \int_{Γ} T \frac{\partial u}{\partial x} d Γ - \int_{C} \frac{\partial}{\partial z} (T \frac{\partial u}{\partial x}) d C + \frac{1}{δ A} \int_{δ A} p δ d (δ A)$

Eq. (30) is in accordance with the extended 3D point-wise J-integral. If we assume constant pressure and infinitely small $δ A,$ the formula would be exactly the same with the extended 3D point-wise J-integral. It indicates that the energy release rate is naturally a point-wise characterizing parameter, and the extended 3D point-wise J-integral can be used to calculate the energy release rate in fluid-driven cracks.