c Allerton Press, Inc., 2008. ISSN 0027-1330, Moscow University Mechanics Bulletin, 2008, Vol. 63, No. 1, pp. 6–12. c A.V. Akulich, A.V. Zvyagin, 2008, published in Vestnik Moskovskogo Universiteta, Matematika. Original Russian Text Mekhanika, 2008, Vol. 63, No. 1, pp. 43–49.
Numerical Simulation of Hydraulic Fracture Crack Propagation A. V. Akulich and A. V. Zvyagin Moscow State University, Faculty of Mechanics and Mathematics, Leninskie Gory, Moscow, 119899 Russia Received July 10, 2006; in final form, October 15, 2007 Abstract—The plane problem of crack motion in an elastic medium under the pressure of a viscous fluid is considered. Under the condition of a constant fluid flow rate, the fluid is injected at the center of the crack. Contrary to other formulations of the problem, this paper attempts to take into account a possible fluid lag behind the crack tip. The resulting numerical solution is compared with a semianalytic one. It is found that the proposed numerical model can be used to predict the characteristics of a hydraulic fracture crack formed in a medium of a prescribed strength.
DOI: 10.3103/S0027133008010020 One of the most efficient ways to enhance the production rate of oil wells is the artificial formation of highpermeability cracks in an oil-bearing stratum; these cracks expand the surface of oil seepage and, as a result, increase the rate of oil extraction from wells. The problem of hydraulic fracturing of oil-bearing formations was first solved in the works of Yu.P. Zheltov and S.A. Khristianovich [1, 2]. This problem remains important and continues to be studied in various formulations. For example, some classes of self-similar solutions are discussed in [3–9]. In this paper we consider the plane problem of crack motion in an elastic medium under the pressure of a viscous fluid. Under the condition of a constant fluid flow rate, the fluid is injected at the center of the crack. Contrary to other formulations of the problem, this paper attempts to take into account a possible fluid lag behind the crack tip. The resulting numerical solution is compared with the semi-analytic one discussed in [9]. These two approaches differ in the way they model the processes near the crack tip. The analytic solution is based on the assumption that the entire crack is filled with a fluid and that the strength of an elastic medium can be ignored, whereas our numerical model assumes that the strength is present and there exists an unknown empty space between the fluid front and the crack tip (this space should be determined in the course of solving the problem). In this paper we show that our numerical results are close to the self-similar ones and that the analytic and numerical solutions are in good agreement when the fracture strength is small. FORMULATION OF THE PROBLEM Let us assume that a viscous incompressible fluid is injected at the crack center x = 0 and Q0 (t) is its prescribed volume flow rate (Fig. 1).
Fig. 1. Geometric interpretation of hydraulic fracture crack propagation.
6
NUMERICAL SIMULATION OF HYDRAULIC FRACTURE
7
It is assumed that the current crack width w(x, t) is small compared to its current half-length l(t); then the problem of fluid motion can be considered in a one-dimensional formulation. In addition to pressure, a viscous force F dx is exerted on a small fluid element of length dx from the crack wall direction. Assuming that the fluid flow is sufficiently slow, we ignore the inertial forces of the fluid. In this case the volume balance equation and the momentum balance equation can be represented in the form ∂w ∂(wV ) + = 0, ∂t ∂x ∂Pf = −F, ∂x
(1)
where t is time, V (x, t) is the velocity of the fluid, Pf is the pressure in the fluid, and F (x, t) is the density of the viscous forces exerted on the fluid length unit from the direction of the crack walls. This density can be taken as that in the problem of flow of a viscous fluid in a bearing: F (x, t) =
12μ V (x, t). w2 (x, t)
(2)
Here μ is the fluid viscosity. It is assumed that pressure is given at infinity, i.e., σxx = σyy = −σ0 . The system of Eqs. (1) and (2) is not closed. If the medium under consideration is elastic, then, in the framework of linear mechanics of cracks and a quasistatic approximation, these equations are solved together with the equilibrium equations ∂σxx ∂σxy + = 0, ∂x ∂y
∂σyx ∂σyy + =0 ∂x ∂y
(3)
and Hooke’s law εxx =
1+ν 1+ν ∂ux ∂uy = (1 − ν)σxx − νσyy , εyy = = (1 − ν)σyy − νσxx , ∂x E ∂y E 1 ∂ux ∂uy 1+ν εxy = εyx = + = σxy . 2 ∂y ∂x E
(4)
Here x and y are the coordinates, σxx , σxy , and σyy are the stress-tensor components, εxx , εxy , and εyy are the strain-tensor components, ux and uy are the displacement vector components, E is Young’s modulus, and ν is Poisson’s ratio. Equations (1), (2) and (3), (4) are related by boundary conditions on the crack faces (y = 0± , |x| < l). The following parameter notation is introduced on the crack faces: ± σij = σij (x, 0± ),
± u± i = ui (x, 0 ),
w(x, t) = u+ (x, 0, t) − u− (x, 0, t).
Then, these boundary conditions take the form y = 0± ,
|x| < l(t),
± σxy = 0,
± σyy = −Pf (x, t).
(5)
± ± Note that the exact boundary condition is σxy = −F rather than σxy = 0. Following [3–6, 8, 9], however, we assume that the force F is much less in absolute value than the pressure Pf (x, t); hence, this force can be ignored. Moreover, it is necessary to take into account the conditions of fluid inflow and a criterion of medium fracture. Let V0 (t) be the velocity of the fluid, w0 (t) be the opening of the crack at the point x = 0, and Q = Q0 /2 be half the fluid flow rate. Since the problem is symmetric, we have
V0 (t) =
Q . w0 (t)
(6)
As a fracture criterion, we adopt the following condition: the stress intensity factor is equal to its critical value (the material strength of the elastic medium). Thus, it is assumed that the crack moves if l l+ξ 1 dξ = KIc , Pf (ξ, t) − σ0 KI = √ (7) l−ξ πl −l
where KI is the stress intensity factor and KIc is its critical value. MOSCOW UNIVERSITY MECHANICS BULLETIN
Vol. 63
No. 1
2008
8
AKULICH, ZVYAGIN
In the framework of our formulation of the problem, the position of the fluid front is assumed to be unknown and is specified by the coordinate x = L(t) (Fig. 1). From the volume conservation law for the injected fluid, we obtain the relation L
t w(x, t) dx =
0
Q(t) dt,
(8)
0
where L(t) is the position of the fluid front in the crack. On this front the pressure should be equal to zero: Pf (L, t) = 0.
x = L(t),
(9)
At infinity, moreover, in the elastic medium the following boundary conditions should be fulfilled: x2 + y 2 → ∞,
σxx = σyy = −σ0 ,
σxy = 0.
(10)
Thus, our problem can be reduced to solving Eqs. (1)–(4) with the boundary conditions (5)–(10). A METHOD OF SOLUTION Our analysis of the problem shows that, at the initial instant of time when the crack length is equal to zero, the pressure in the crack should be infinite; when the problem is solved numerically, therefore, it is necessary to choose an initial approximation at a nonzero instant t = t0 in order to start computations. As such initial conditions, we take the exact solution obtained in [8]: Q3/4 √ , 4.39 l0 (μ/E )1/4 1/3 1/3 μV0 (E )2/3 V0 μ1/3 w0 (x) = 2(l0 − x)2/3 37/6 , p (x) = p − . 0 E (l0 − x)1/3 31/3 √ 1/3 6 2 (E )2 μV0 E , E = , V0 is the initial velocity of the fluid, w0 (x) is the initial Here p = σ0 + π 3l0 1 − ν2 opening profile of the crack, and p0 (x) is the initial pressure profile. Our computations were performed for l0 = 1 m (this length corresponded to t0 = 0.63 s). In this paper we use the displacement discontinuity method [10, 11] to find the displacements and stresses of an elastic medium. This method is very efficient for solving the problems of crack mechanics. The massconservation equation (1) is linear in x and can be easily integrated. At each instant of time, then, the velocity V (x, t) can be expressed by the formula
x ˙ t) dx c − 0 w(x, . (11) V (x, t) = w(x, t) V0 =
Taking into account the boundary condition (6) in the fluid injection zone, we get x = 0,
V (0, t) =
Q . w0 (t)
From here we conclude that in (11) the integration constant c is equal to Q. In order to determine the fluid velocity inside the crack, we assume that the crack opening can be represented in the form (see [7]) x w(x, t) = w0 (t)f (12) = w0 f (ξ), l(t) where ξ =
x is the dimensionless coordinate. In (11) the integral takes the form l(t) x
ξ f (ξ) dξ − w0 l˙
w(x, ˙ t) dx = w˙ 0 l 0
ξ
0
ξ f (ξ) dξ − w0 ξf (ξ)l.˙
˙ ξf (ξ) dξ = (w˙ 0 l + w0 l) 0
(13)
0
MOSCOW UNIVERSITY MECHANICS BULLETIN
Vol. 63
No. 1
2008
NUMERICAL SIMULATION OF HYDRAULIC FRACTURE
9
From the volume conservation law it follows that for the constant fluid flow rate Q we have 1 w0 l
f (ξ) dξ = Qt,
w˙ 0 l + w0 l˙ = 1 0
0
Q
.
(14)
f (ξ) dξ
Substituting (14) into (13) and (12), we can find the velocity distribution as a function of Q and l˙ and the dependence of the fluid velocity profile on the crack opening shape:
ξ Q 0 f (ξ) dξ ˙ V (ξ, t) = 1− 1 + lξ. (15) w f (ξ) dξ 0
Here l˙ is the velocity of the crack tip. A NUMERICAL METHOD Below we use the fact that the elasticity theory problem is elliptic; therefore, its solution depends continuously on the parameters in use. In our formulation of the problem, time is a parameter, since all equations are assumed to be static (the only exception is the continuity equation). It should be taken into account that Eqs. (1) and (2) describing the fluid motion and the static elasticity theory equations (3) and (4) should be simultaneously solved at each instant of time. In our computations we used a method of successive approximations. As an initial approximation, a selfsimilar solution [7] was adopted at each instant of time. In order to enhance the accuracy of computations in the case when the general number of boundary elements remains the same, we used a nonuniform partition (the smallest boundary elements were situated near the crack tip). To accomplish this, in the change ξ = sin ϕ we used a variable step in ϕ, where 0 ≤ ϕ ≤ π/2. Equation (15) can be represented in the following discrete form for the chosen partition along the crack length: Q I k (ξ N ) k k ˙ (16) V =ξ L+ k 1− k k . w I1 (ξ ) Here I k (ξ N ) =
N k
wk−1 + wk dξ k k k wk−1 + wk dξ k , I1 (ξ ) = , and k is the node number. In order to 2 w0 2 w0 m=1 m=1
find wk , the value of wk determined at the previous time step is used as a first approximation. To integrate the momentum conservation equation (1), we rewrite it in the following discrete form: k−1 k−1 2 V w +Vk + wk pk = pk−1 − 12μ L dξ k−1 . (17) 2 2 The displacement discontinuity method is used to solve the elasticity theory problem (3), (4). This method is based on the representation of a solution to the above problem in the form of a sum of basis analytic solutions with unknown coefficients. These coefficients are determined from the boundary conditions on the discrete set of boundary points. As basis solutions, we take the solutions to the constant displacement discontinuity problems on a finite interval in an elastic plane (in a boundary element). The crack is divided into N intervals; each of these intervals is used as a boundary element. If an analytic solution is known for the unit discontinuity of displacements, then we find a numerical solution by summing the effects of all N elements. A detailed description of the above method can be found in [7]. For each time step, the problem is solved by successive refinements of the chosen initial approximation to minimize the error in the numerical solution. The following iteration parameters are chosen for a fixed value ˙ and the fluid front velocity L(t). ˙ of t: the crack velocity l(t) If the solution at the previous time step is known ˙ and L(t) ˙ and the velocities l(t) are given, then Eqs. (16) and (17) and Eqs. (3) and (4) can be integrated with consideration of the boundary conditions (5)–(10) (note that this solution is used to construct an initial approximation). The above two iteration parameters are found from the fracture criterion at the crack tip and from the condition that the fluid volume in the crack is equal to the volume of the injected fluid: L(t) t w(x, t) dξ − (t) dt < ε2 . |KI − KIc | < ε1 , 0 0 MOSCOW UNIVERSITY MECHANICS BULLETIN
Vol. 63
No. 1
2008
10
AKULICH, ZVYAGIN
Here ε1 and ε2 are the corresponding errors. The above computational scheme was implemented as a computer program. In order to verify its correct operation, a number of test problems were solved and the thus-obtained numerical results were compared with the results of other authors. NUMERICAL RESULTS AND THEIR ANALYSIS A self-similar solution to the hydraulic fracture crack problem is solved in [4, 9] for the case of a constant fluid flow rate under the following assumptions: (i) a crack is entirely filled with an injected fluid (in other words, the fluid front coincides with the crack tip); (ii) the fracture strength of the elastic medium under study is absent (in other words, the critical stress intensity factor is equal to zero). Note that the first assumption leads to an infinitely large negative pressure in the fluid at the crack tip. The second assumption is not always justified physically in the case of high-strength rocks. A detailed description of the self-similar solution is given in [4]. Since the results obtained for the self-similar solution are represented in terms of dimensionless variables, these results are interpreted below to compare with our numerical results. The following quantities were introduced to solve the self-similar problem: T =
12μ , E
L =
Q0 T .
(18)
Here T is the time scale and L is the length scale. The dimensionless coordinate ξ and the dimensionless time τ were defined as t x (19) τ= . ξ = , L T The dimensionless crack length χ(τ ), the dimensionless pressure Π(ξ, τ ), and the dimensionless crack opening Ω(ξ, t) were introduced in a similar way: χ=
L , L
Π=
Pf − σ0 , E
Ω=
w . L
(20)
The self-similar solution for opening and for pressure is Ω(ξ, τ ) = λτ 1/3 Ωξ (ξ),
Π(ξ, τ ) = τ −1/3 Πξ (ξ),
(21)
where the function Πξ (ξ) is represented by a double series and is of a rather complicated form [4]. According to [4], we have λ ∼ = 0.629. For comparison, the fluid flow rate and the parameters for the elastic medium and for the fluid were chosen as in [4]; the curves of various quantities were constructed depending on time or on a coordinate for a fixed value of time. Our numerical results were represented in the form (18)–(22) and were compared with the analytic ones. The following input data were used in our numerical model: the parameters of the elastic medium were E = 25 GPa, ν = 0.15, σ0 = 50 MPa, and Q0 = 4 × 10−3 m2 /s; the fluid viscosity μ was taken equal to 8.53 × 10−7 MPa s; and the medium strength KIc was taken equal to 0.001, 0.1, 1, 5, and 10 (these values were expressed in MPa m1/2 ). In its dimensionless form, the strength was found by the formula K=
KIc . E (L )1/2
(22)
In other words, in our computations the following values of the dimensionless strength were used: K = 3.5 × 10−5 , 3.5 × 10−3 , 3.5 × 10−2 , 1.7 × 10−1 , and 3.5 × 10−1 . The time dependencies of the pressure and the crack opening at the center of the crack and the time dependence of the crack length are illustrated in Fig. 2. Our comparison of the numerical and self-similar solutions for K up to 1.7 × 10−1 shows that these solutions are in good agreement.
MOSCOW UNIVERSITY MECHANICS BULLETIN
Vol. 63
No. 1
2008
NUMERICAL SIMULATION OF HYDRAULIC FRACTURE
11
Fig. 2. The solution at the fluid injection point. The dependence of (I ) the dimensionless ˆ = 104 Π = 104 Pf − σ0 , (II ) the dimensionless crack length χ ˆ = 10−7 χ = pressure Π E l ˆ = 10−3 Ω = 10−3 w on the dimen10−7 , and (III ) the dimensionless crack opening Ω L L −12 −12 t . The solid line corresponds to the analytic solution. sionless time τˆ = 10 τ = 10 T
The curves of Πξ (ξ) and Ωξ (ξ) are represented in Fig. 3a for several time values when K = 3.5 × 10−5 . A significant distinction between the analytic solution and the numerical solution can be observed for small values of time when the medium strength and the fluid lag behind the crack tip exert a noticeable influence on the process. It is interesting to note that these solutions coincide again when the crack length is large. If strength is very small, then these solutions are almost the same; this fact is a good test for our computer program. The results obtained for K = 3.5 × 10−3 , 3.5 × 10−2 , and 1.7 × 10−1 are almost the same as those obtained for K = 3.5 × 10−5 .
Fig. 3. The solutions for the dimensionless strength (a) K = 3.5 × 10−5 and (b) K = 3.5 × 10−1 . The ˆ ξ (ξ) = 10Πξ (ξ) = 10τ 1/3 Π(ξ, τ ) and (II ) the dependence of (I ) the truncated dimensionless pressure Π 1 x truncated dimensionless crack opening Ωξ (ξ) = τ −1/3 Ω(ξ, τ ) on the dimensionless coordinate ξ = . λ L The solid line corresponds to the analytic solution.
The above curves are represented in Fig. 3b for K = 3.5 × 10−1 . In this case the numerical results are in wide disagreement with the semi-analytic solution. This can be explained by the fact that the analytic solution obtained without consideration of strength is unacceptable. The results of this paper show that the above numerical model can be used to predict the hydraulic fracture crack characteristics in a medium of a given strength. Its advantage is that it can also be used in the case when the hydraulic fracture crack is curvilinear or changes its direction under the effect of medium inhomogeneities (for example, when there exists a natural fracture of rock). MOSCOW UNIVERSITY MECHANICS BULLETIN
Vol. 63
No. 1
2008
12
AKULICH, ZVYAGIN
REFERENCES 1. S. A. Khristianovich and Yu. P. Zheltov, “Formation of Vertical Fractures by Means of Highly Viscous Liquid,” in Proc. 4th World Petrol. Congr., June 6–16, 1995 (Carlo Colombo Publishers, Rome, 1955), Vol. 2, pp. 579–586. 2. Yu. P. Zheltov and S. A. Khristianovich, “On the Hydraulic Fracturing of an Oil-Bearing Stratum,” Izv. Akad. Nauk SSSR, Otd. Tekh. Nauk, No. 5, 3–41 (1955). 3. T. K. Perkins and L. R. Kern, “Widths of Hydraulic Fractures,” J. Petrol. Technol. 13 (9), 937–949 (1961). 4. J. I. Adachi and E. Detournay, “Self-Similar Solution of a Plane-Strain Fracture Driven by a Power-Law Fluid,” Int. J. Numer. Anal. Meth. Geomech. 26, 579–604 (2002). 5. D. A. Spence and P. Sharp, “Self-Similar Solution for Elastohydrodynamic Cavity Flow,” Proc. Roy. Soc. London. Ser. A 400, 289–313 (1985). 6. O. E. Ivashnev and N. N. Smirnov, “Formation of Hydraulic Fracture Cracks in a Porous Medium,” Vestn. Mosk. Univ. Ser. 1: Mat. Mekh., No. 6, 28–36 (2003) [Moscow Univ. Mech. Bull. 58 (6), 12–20 (2003)]. 7. A. V. Zvyagin, “Motion of a Viscous Fluid in a Channel with Elastic Boundaries,” Vestn. Mosk. Univ. Ser. 1: Mat. Mekh., No. 1, 50–54 (2005) [Moscow Univ. Mech. Bull. 60 (1), 9–13 (2005)]. 8. J. Desroches, E. Detournay, B. Lenoach, et al., “The Crack Tip Region in Hydraulic Fracturing,” Proc. Roy. Soc. London. Ser. A 447, 39–48 (1994). 9. R. Carbonell, J. Desroches, and E. Detournay, “A Comparison Between a Semi-Analytical and a Numerical Solution of a Two-Dimensional Hydraulic Fracture,” Int. J. Solids Struct. 36, 4869–4888 (1999). 10. S. L. Crouch and A. M. Starfield, Boundary Element Methods in Solid Mechanics (Allen & Unwin, London, 1983). 11. A. V. Bogdanov, A. V. Zvyagin, and M. Thiercelin, “Mutual Influence of a Crack System and the Stress Intensity Factor,” Vestn. Mosk. Univ. Ser. 1: Mat. Mekh., No. 6, 44–49 (2004) [Moscow Univ. Mech. Bull. 59 (6), 13–18 (2004)].
Translated by O. Arushanyan
MOSCOW UNIVERSITY MECHANICS BULLETIN
Vol. 63
No. 1
2008
