**1. Introduction**

12 Effective and Sustainable Hydraulic Fracturing

640 Effective and Sustainable Hydraulic Fracturing

Cambridge UK, 1967.

edition, 1959.

[16] J. R. Lister. Buoyancy-driven fluid fracture: The effects of material toughness and of

[17] G. K. Batchelor. *An Introduction to Fluid Dynamics*. Cambridge University Press,

[19] H. Carslaw and J. C. Jaeger. *Conduction of Heat in Solids*. Oxford University Press, 2nd

[20] H. Stehfest. Numerical inversion of Laplace transform. *Communications of the Association*

[21] L.X. Nghiem, P.A. Forsyth Jr., and A. Behie. A fully implicit hydraulic fracture model.

[22] D. Garagash and E. Detournay. An analysis of the influence of the pressurization rate on the borehole breakdown pressure. *Int. J. Solids Structures*, 34(24):3099–3118, 1997.

low-viscosity precursors. *J. Fluid Mech.*, 210:263–280, 1990.

[18] J. Bear. *Dynamics of Fluids in Porous Media*. Dover Publications, 1988.

*for Computing Machinery (ACM)*, 13:47–49, January 1970.

*Journal of Petroleum Technology*, 36(7):1191–1198, July 1984.

Analytical models and numerical simulation are important means to increase understanding and to enhance efficiency of hydraulic fracturing. The reasons for developing and using them are clearly explained, for instance, by Mack and Warpinski [1]. Thus, there is no need to dwell on them. Rather we focus on improving analytical and numeric methods used to the date. Our objective is to suggest new approaches for developing accurate, robust and stable simulators on the basis of recent analytical and computational findings [2-6].

The approaches discussed in the paper stem from the fact [2,3] that the conventional formu‐ lation of the hydraulic fracture problem (for example, see [7]), when neglecting the lag and fixing the position of the fracture front at a time step, is ill-posed. This feature has not been reported for more than three decades of studying hydraulic fractures because of two rea‐ sons. Numerical simulators, based on the conventional formulation (for example, see [7,8]), employ quite rough meshes, which themselves serve as specific 'regularizators'. On the oth‐ er hand, rare solutions of model problems either also employed rough meshes [9], or they were obtained by solving the initial value (Cauchy) problem [10-12] rather than the boun‐ dary value problem (a discussion of the difference may be found in references [3,4]). The disclosure of the mentioned fact has led to (i) explicit formulation of the speed equation (SE) in its general form1 , (ii) comprehension of its significance for proper numerical simulation of hydraulic fractures and (iii) distinguishing the particle velocity as a preferable variable2 . It

<sup>1</sup> To the authors'knowledge, Kemp [11] was the first who clearly distinguished the speed equation when revisiting the Nordgren's problem. When introducing the SE, numbered (5) in his paper, Kemp wrote (p. 289): "Nowhere is (5) mentioned. (5) is called the Stefan condition and is always present in moving boundary problems". Kemp used the SE to solve the Nordgren's problem as an initial value rather than boundary value problem. This excluded solving ill-posed boundary value problem.

<sup>© 2013</sup> Linkov and Mishuris; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

has also led to the efficient means, called ε–regularization [2,3], to overcome the analytical and computational difficulties caused by ill-posedness. Finally, the entire conventional for‐ mulation of the problem has changed to the modified formulation, which opens new analyt‐ ical and computational options for solving hydraulic fracture problems.

general case, keeping the velocity, which is a smooth function near the fluid front, as un‐ known is more convenient than to replace it by grad*p*, which is commonly strongly singu‐ lar. Note only that the operator resulting from substitution of (2) into (1) is of first order in time and of second order and elliptic in spatial coordinates. Consequently it requires one *initial* condition defining the initial distribution of the opening. In terms of the modified

Modified Formulation, ε-Regularization and the Efficient Solution of Hydraulic Fracture Problems

The elliptic operator requires only one *boundary* condition (BC) on the fluid front *Le*. For in‐ stance, when neglecting the lag between the fluid front and the fracture contour, the BC may

where q₀(x) is a known function at *L <sup>e</sup>*; specifically, at the injection points, it is defined by the injection regime. At a point *x\** of the propagating fluid front, coinciding with the fracture contour *L <sup>f</sup>* , we have *w* = 0 and equation (2) implies *q <sup>f</sup>* (*x*, *t*)=*wv* =0. Then (4) becomes:

Of special importance is that the limit value of the particle velocity *vn\** at the fluid front rep‐

Herein, *xn\** is the normal component of a point *x\** at the front. It is assumed that sucking or evaporation through the front is negligible. (6) is the *speed equation* (SE), which is fundamen‐ tal for moving boundaries. Using (2) in (6) specifies the SE for a flow of incompressible fluid

) *<sup>D</sup>*(*w*, *<sup>p</sup>*) <sup>∂</sup> *<sup>p</sup>*

Thus we have the *local* condition (7) at points of the propagating fluid front. This allows one to trace the propagation by well-established methods of the theory of propagating surfaces [17]. In contrast, the conventional formulation employs the *global* mass balance (for example, see [7,8,10,12]), which is a single equation. The latter is sufficient when considering 1D prob‐ lems with one point of the front to be traced. However, in the general case of a 2D fracture, it is preferable to employ the SE, which is formulated at each of the points of the fluid front.

∂ *n x*=*x*\*

*v\** <sup>=</sup> *<sup>d</sup> xn\**

be the condition of the prescribed normal component *qn*(*x*, *<sup>t</sup>*)= *<sup>y</sup> <sup>α</sup>vn* of the flux:

*<sup>y</sup> <sup>α</sup>vn*(*x*, *<sup>t</sup>*)=*q*<sup>0</sup>

resents the speed of the front propagation *v\** [2-4]:

*<sup>v</sup>*\*<sup>=</sup> *<sup>d</sup> xn*\*

*dt* <sup>=</sup> - <sup>1</sup> *w*(*x*\*

in a narrow channel [2,3]:

*y*(*x*, *t*0) = *y*0(*x*) (3)

http://dx.doi.org/10.5772/56218

643

(*x*, *t*), *x* ∈ *L <sup>e</sup>* (4)

*dt* =*vn\** (6)

(7)

*y*(*x\**, *t*) =0, *x\** ∈ *L <sup>f</sup>* (5)

opening it is:

Below we employ these options. Section 2 contains a concise review of the modified formu‐ lation. In Section 3, we illustrate its analytical advantages by simple solutions for the Per‐ kins-Kern-Nordgren (PKN) [9,13] and Khristianovich-Geertsma-de Klerk (KGD) [14,15] models. In Section 4, we turn to computational advantages of the modified formulation, present new computational schemes and illustrate their efficiency by numerical results for the PKN model. Section 5 contains the extension of these schemes to the pseudo-three-di‐ mensional (P3D) models [1], the importance of which grows nowadays because of their em‐ ployment in simulators for hydraulic fractures in low permeable shales [16]. A brief summary concludes the paper (Section 6).

#### **2. Modified formulation of hydraulic fracture problem**

The modified equations [3,4] use as variables the velocity *v* of fluid particles averaged through the channel width (fracture opening) and the modified opening *y* =*w* 1/*α*, where *α* is the exponent defined by the asymptotic behavior of the opening *w* at the fluid front (*w* =*C*(*t*)*r <sup>α</sup>*, *r* is the distance from the front, *t* is the time). Then the continuity equation for a flow in a narrow channel reads:

$$\frac{\partial y}{\partial t} + \mathbf{v} \bullet \mathbf{grad} y + \frac{1}{\alpha} y \text{div} \mathbf{v} + \frac{y^{\text{l-a}}}{\alpha} q\_l = 0 \tag{1}$$

where *ql* is the intensity of distributed sinks or sources (below this term is assumed positive to account for leak off); the divergence and gradient are defined in the tangent plane to the channel (hydraulic fracture). The Poiseuille type relation in terms of the particle velocity is:

$$\mathbf{v} = -\frac{1}{w} D(w, p) \mathbf{grad} p \tag{2}$$

It is obtained by integration of the Navier-Stokes equations for a flow of a viscous fluid in a narrow channel. Herein, *p*(*x,t*) is the fluid pressure, *x* is the position vector on the surface of the flow in immovable coordinates, *D* is a function or operator, such that *D*(0, *p*)grad*p* =0.

The equations (1) and (2) represent the modified system of PDE. Substitution of (2) into (1) gives the modified lubrication equation. We shall not write it down explicitly because in the

<sup>2</sup> Mack and Warpinski [1] have noticed the beneficial property of the velocity. On p. 6-21 of their fundamental work they wrote: "Detailed numerical simulations have shown that the velocity varies much more slowly than the flow rate". These authors made the best of this property by using the velocity as the unknown variable; they actually employed the speed equations, as well, although not writing it explicitly.

general case, keeping the velocity, which is a smooth function near the fluid front, as un‐ known is more convenient than to replace it by grad*p*, which is commonly strongly singu‐ lar. Note only that the operator resulting from substitution of (2) into (1) is of first order in time and of second order and elliptic in spatial coordinates. Consequently it requires one *initial* condition defining the initial distribution of the opening. In terms of the modified opening it is:

has also led to the efficient means, called ε–regularization [2,3], to overcome the analytical and computational difficulties caused by ill-posedness. Finally, the entire conventional for‐ mulation of the problem has changed to the modified formulation, which opens new analyt‐

Below we employ these options. Section 2 contains a concise review of the modified formu‐ lation. In Section 3, we illustrate its analytical advantages by simple solutions for the Per‐ kins-Kern-Nordgren (PKN) [9,13] and Khristianovich-Geertsma-de Klerk (KGD) [14,15] models. In Section 4, we turn to computational advantages of the modified formulation, present new computational schemes and illustrate their efficiency by numerical results for the PKN model. Section 5 contains the extension of these schemes to the pseudo-three-di‐ mensional (P3D) models [1], the importance of which grows nowadays because of their em‐ ployment in simulators for hydraulic fractures in low permeable shales [16]. A brief

The modified equations [3,4] use as variables the velocity *v* of fluid particles averaged through the channel width (fracture opening) and the modified opening *y* =*w* 1/*α*, where *α* is the exponent defined by the asymptotic behavior of the opening *w* at the fluid front (*w* =*C*(*t*)*r <sup>α</sup>*, *r* is the distance from the front, *t* is the time). Then the continuity equation for a

1

*<sup>α</sup> <sup>y</sup>*div**<sup>v</sup>** <sup>+</sup> *<sup>y</sup>* 1-*<sup>α</sup>*

to account for leak off); the divergence and gradient are defined in the tangent plane to the channel (hydraulic fracture). The Poiseuille type relation in terms of the particle velocity is:

It is obtained by integration of the Navier-Stokes equations for a flow of a viscous fluid in a narrow channel. Herein, *p*(*x,t*) is the fluid pressure, *x* is the position vector on the surface of the flow in immovable coordinates, *D* is a function or operator, such that *D*(0, *p*)grad*p* =0.

The equations (1) and (2) represent the modified system of PDE. Substitution of (2) into (1) gives the modified lubrication equation. We shall not write it down explicitly because in the

2 Mack and Warpinski [1] have noticed the beneficial property of the velocity. On p. 6-21 of their fundamental work they wrote: "Detailed numerical simulations have shown that the velocity varies much more slowly than the flow rate". These authors made the best of this property by using the velocity as the unknown variable; they actually employed the speed

is the intensity of distributed sinks or sources (below this term is assumed positive

*<sup>α</sup> ql* =0 (1)

*<sup>w</sup> D*(*w*, *p*)grad*p* (2)

ical and computational options for solving hydraulic fracture problems.

**2. Modified formulation of hydraulic fracture problem**

<sup>∂</sup> *<sup>t</sup>* + **v**∙grad*y* +

**<sup>v</sup>**<sup>=</sup> - <sup>1</sup>

summary concludes the paper (Section 6).

642 Effective and Sustainable Hydraulic Fracturing

flow in a narrow channel reads:

equations, as well, although not writing it explicitly.

where *ql*

∂ *y*

$$y(x, t\_0) = y\_0(x) \tag{3}$$

The elliptic operator requires only one *boundary* condition (BC) on the fluid front *Le*. For in‐ stance, when neglecting the lag between the fluid front and the fracture contour, the BC may be the condition of the prescribed normal component *qn*(*x*, *<sup>t</sup>*)= *<sup>y</sup> <sup>α</sup>vn* of the flux:

$$y^{\alpha}v\_n(\mathbf{x},t) = q\_0(\mathbf{x},t), \qquad \qquad \mathbf{x} \in L\_{\epsilon} \tag{4}$$

where q₀(x) is a known function at *L <sup>e</sup>*; specifically, at the injection points, it is defined by the injection regime. At a point *x\** of the propagating fluid front, coinciding with the fracture contour *L <sup>f</sup>* , we have *w* = 0 and equation (2) implies *q <sup>f</sup>* (*x*, *t*)=*wv* =0. Then (4) becomes:

$$\mathbf{y}(\mathbf{x}\_{\ast},t) = \mathbf{0}, \qquad \qquad \qquad \mathbf{x}\_{\ast} \in L\_{\frac{\pi}{f}} \tag{5}$$

Of special importance is that the limit value of the particle velocity *vn\** at the fluid front rep‐ resents the speed of the front propagation *v\** [2-4]:

$$
\omega\_{\mathcal{U}\*} = \frac{d \,\, \mathbf{x}\_{n^\*}}{dt} = \mathcal{U}\_{n^\*} \tag{6}
$$

Herein, *xn\** is the normal component of a point *x\** at the front. It is assumed that sucking or evaporation through the front is negligible. (6) is the *speed equation* (SE), which is fundamen‐ tal for moving boundaries. Using (2) in (6) specifies the SE for a flow of incompressible fluid in a narrow channel [2,3]:

$$\mathbf{w}\_{\ast} = \frac{d\,\mathbf{x}\_{\ast}\prime}{dt} = \mathbf{-}\left[\frac{1}{w\{\mathbf{x}\_{\ast}\}}D\begin{pmatrix}\mathbf{w}\_{\prime} & p\end{pmatrix}\frac{\partial p}{\partial \mathbf{n}}\right]\_{\mathbf{x}=\mathbf{x}\_{\ast}}\tag{7}$$

Thus we have the *local* condition (7) at points of the propagating fluid front. This allows one to trace the propagation by well-established methods of the theory of propagating surfaces [17]. In contrast, the conventional formulation employs the *global* mass balance (for example, see [7,8,10,12]), which is a single equation. The latter is sufficient when considering 1D prob‐ lems with one point of the front to be traced. However, in the general case of a 2D fracture, it is preferable to employ the SE, which is formulated at each of the points of the fluid front. This gives another evidence that using the particle velocity is beneficial from the computa‐ tional point of view.

The speed equation (6) also yields important implications for numerical simulation of hy‐ draulic fractures by finite differences. Indeed, when at a time step we have known both *x\** and *v\** (*x\**), the SE (7) becomes a boundary condition additional to the boundary condition (5) on the front. Thus, as noted in [2], a boundary value problem may appear overdeter‐ mined and ill-posed in the Hadamard sense [18]. To avoid difficulties, it is reasonable to use *ε*-regularization, suggested and employed in [2,3].

The *ε*-regularization is performed as follows. At each point of the fluid front, an exact BC is changed to an approximate equality at a small distance *rε* behind the front:

$$\int\_{p\_0}^{p\_\star} \frac{1}{w} D\left(w, p\right) dp = v\_\* r\_\varepsilon \tag{8}$$

*KI* = *KIc*, *KII* =0 (11)

http://dx.doi.org/10.5772/56218

645

Modified Formulation, ε-Regularization and the Efficient Solution of Hydraulic Fracture Problems

<sup>∂</sup> *<sup>x</sup>* )1/*<sup>n</sup>* (12)

*p* =*krw* (13)

<sup>∂</sup> *<sup>x</sup>* )1/*<sup>n</sup>* (14)

=*E* /(1 - *ν* 2), *E* is the elasticity modulus, *ν* is the Poisson's ratio, *h* is

(for example, see [12]), the *x-*axis is in

where *KI* is the tensile stress intensity factor (SIF), *KIc* is its critical value, *KII* is the shear SIF. The modified problem, when neglecting the lag (*xe* = *x\**), consists of solving the PDE (1) and (2) together with the solid mechanics equation (10) under the initial condition (3), BC (4) at the contour part with prescribed influx, and the *ε*-regularized BC (8) at the fluid (and frac‐ ture) front. The conditions (11) define the possibility and the direction of the fracture propa‐ gation. The *ε*-regularized SE (9) serves to find the front propagation on the time steps.

In particular cases of 1D problems for the PKN [9,13] and KGD [14,15] models, one can solve an initial value rather than boundary value problem (for example, see [10-12]). Then the problem is well-posed and it does not require regularization. Still in these cases, the modi‐ fied variables are of use to obtain simple analytical solutions of problems, which otherwise requite involved calculations. For a Newtonian fluid, the analytical solutions are given in [4] both for the PKN and KGD models. In a similar way, by employing the modified variables,

Consider a fluid with the viscosity law of power-type *τ* =*M γ*˙ *<sup>n</sup>*, where *τ* is the shear stress, *γ*˙ is the shear strain rate, *M* is the consistency index, *n* is the behaviour index. The common derivations for a flow in a narrow channel yield the Poiseuille type dependence between the

the direction of the fracture propagation, *v* is the component of the particle velocity in the

The geometrical scheme of the PKN model is given in Figure 1. In this model the planestrain conditions occur in cross-sections parallel to the fracture front. Then the elasticity

∂ *w <sup>n</sup>*+2

particle velocity *v*, averaged through the channel width *w*, and the pressure gradient:

*<sup>v</sup>* =(-*<sup>k</sup> <sup>f</sup> <sup>w</sup> <sup>n</sup>*+1 <sup>∂</sup> *<sup>p</sup>*

we may solve these 1D problems for a non-Newtonian fluid.

where *<sup>k</sup> <sup>f</sup>* =1 /(*θ<sup>M</sup>* ), for a plane channel, *<sup>θ</sup>* =2 2(*<sup>n</sup>* <sup>+</sup> 1) *<sup>n</sup>*

the fracture height. Substitution of (13) into (12) yields:

*v* =( *k <sup>f</sup> kr n* + 2

**3. Analytical solutions**

direction of *x*.

equation (10) takes the form [9]:

where *kr* = 2 /(*πh* ) *E'*, *E '*

where *p*<sup>0</sup> is the pressure at the front, *p<sup>ε</sup>* is the pressure at a distance *rε* from the front. Equa‐ tion (8) is obtained by combining the boundary condition at the fluid front, particular for a problem, with the SE, which is general. The distance (absolute *rε* or relative *ε*) is taken small enough to use the equality sign in the derived approximate condition. The SE is also as‐ sumed to be met at the distance *rε* to an accepted accuracy:

$$\text{Cov}\_\*(t) = \frac{d \,\, \text{x}\_{w^\*}}{dt} = -\left[\frac{1}{w} \, D\begin{pmatrix} w & p \end{pmatrix} \frac{\partial \, p}{\partial w}\right]\_{t = \pi\_x} \tag{9}$$

The *ε*-regularized BC (8) allows one to avoid unfavorable computational effects. The *ε*-regu‐ larized SE (9) serves to find the front propagation. In particular, it is basic when applying the level set and fast marching methods [17].

The described modification concerns mostly with the *fluid* equations, which are sufficient when the opening *w* is known. However, it is not known in advance in problems of hydraul‐ ic fracturing. To find it, the fluid equations are complimented with solid mechanics and frac‐ ture equations.

In terms of the modified opening, the *solid mechanics* equation with an integral operator *A*

$$A\,y^{\alpha} = p\tag{10}$$

is solved under the BC of zero opening at points of the fracture contour. When neglecting the lag, this condition coincides with the condition (5) of zero flux at the fluid front.

The *fracture mechanics* equations define the critical state and the perspective direction of the fracture propagation. In the case of tensile mode of fracture, these are:

Modified Formulation, ε-Regularization and the Efficient Solution of Hydraulic Fracture Problems http://dx.doi.org/10.5772/56218 645

$$\mathbf{K}\_{I} = \mathbf{K}\_{Ic'} \qquad \quad \quad \mathbf{K}\_{II} = \mathbf{0} \tag{11}$$

where *KI* is the tensile stress intensity factor (SIF), *KIc* is its critical value, *KII* is the shear SIF. The modified problem, when neglecting the lag (*xe* = *x\**), consists of solving the PDE (1) and (2) together with the solid mechanics equation (10) under the initial condition (3), BC (4) at the contour part with prescribed influx, and the *ε*-regularized BC (8) at the fluid (and frac‐ ture) front. The conditions (11) define the possibility and the direction of the fracture propa‐ gation. The *ε*-regularized SE (9) serves to find the front propagation on the time steps.

#### **3. Analytical solutions**

This gives another evidence that using the particle velocity is beneficial from the computa‐

The speed equation (6) also yields important implications for numerical simulation of hy‐ draulic fractures by finite differences. Indeed, when at a time step we have known both *x\**

(5) on the front. Thus, as noted in [2], a boundary value problem may appear overdeter‐ mined and ill-posed in the Hadamard sense [18]. To avoid difficulties, it is reasonable to use

The *ε*-regularization is performed as follows. At each point of the fluid front, an exact BC is

where *p*<sup>0</sup> is the pressure at the front, *p<sup>ε</sup>* is the pressure at a distance *rε* from the front. Equa‐ tion (8) is obtained by combining the boundary condition at the fluid front, particular for a problem, with the SE, which is general. The distance (absolute *rε* or relative *ε*) is taken small enough to use the equality sign in the derived approximate condition. The SE is also as‐

*<sup>w</sup> <sup>D</sup>*(*w*, *<sup>p</sup>*) <sup>∂</sup> *<sup>p</sup>*

The *ε*-regularized BC (8) allows one to avoid unfavorable computational effects. The *ε*-regu‐ larized SE (9) serves to find the front propagation. In particular, it is basic when applying

The described modification concerns mostly with the *fluid* equations, which are sufficient when the opening *w* is known. However, it is not known in advance in problems of hydraul‐ ic fracturing. To find it, the fluid equations are complimented with solid mechanics and frac‐

In terms of the modified opening, the *solid mechanics* equation with an integral operator *A*

is solved under the BC of zero opening at points of the fracture contour. When neglecting

The *fracture mechanics* equations define the critical state and the perspective direction of the

the lag, this condition coincides with the condition (5) of zero flux at the fluid front.

fracture propagation. In the case of tensile mode of fracture, these are:

∂ *n r*=*r<sup>ε</sup>*

( ) <sup>0</sup> \* <sup>1</sup> , *<sup>e</sup> <sup>p</sup> <sup>p</sup> D w p dp v r <sup>w</sup>*

changed to an approximate equality at a small distance *rε* behind the front:

(*x\**), the SE (7) becomes a boundary condition additional to the boundary condition

e

<sup>=</sup> ò (8)

*Ay <sup>α</sup>* = *p* (10)

(9)

tional point of view.

644 Effective and Sustainable Hydraulic Fracturing

*ε*-regularization, suggested and employed in [2,3].

sumed to be met at the distance *rε* to an accepted accuracy:

(*t*)= *<sup>d</sup> xn\**

*dt* <sup>=</sup> - <sup>1</sup>

*v\**

the level set and fast marching methods [17].

ture equations.

and *v\**

In particular cases of 1D problems for the PKN [9,13] and KGD [14,15] models, one can solve an initial value rather than boundary value problem (for example, see [10-12]). Then the problem is well-posed and it does not require regularization. Still in these cases, the modi‐ fied variables are of use to obtain simple analytical solutions of problems, which otherwise requite involved calculations. For a Newtonian fluid, the analytical solutions are given in [4] both for the PKN and KGD models. In a similar way, by employing the modified variables, we may solve these 1D problems for a non-Newtonian fluid.

Consider a fluid with the viscosity law of power-type *τ* =*M γ*˙ *<sup>n</sup>*, where *τ* is the shear stress, *γ*˙ is the shear strain rate, *M* is the consistency index, *n* is the behaviour index. The common derivations for a flow in a narrow channel yield the Poiseuille type dependence between the particle velocity *v*, averaged through the channel width *w*, and the pressure gradient:

$$w = \left( \*k\_f \, w^{n+1} \frac{\partial \, p}{\partial x} \right)^{1/n} \tag{12}$$

where *<sup>k</sup> <sup>f</sup>* =1 /(*θ<sup>M</sup>* ), for a plane channel, *<sup>θ</sup>* =2 2(*<sup>n</sup>* <sup>+</sup> 1) *<sup>n</sup>* (for example, see [12]), the *x-*axis is in the direction of the fracture propagation, *v* is the component of the particle velocity in the direction of *x*.

The geometrical scheme of the PKN model is given in Figure 1. In this model the planestrain conditions occur in cross-sections parallel to the fracture front. Then the elasticity equation (10) takes the form [9]:

$$p = k\_r w \tag{13}$$

where *kr* = 2 /(*πh* ) *E'*, *E '* =*E* /(1 - *ν* 2), *E* is the elasticity modulus, *ν* is the Poisson's ratio, *h* is the fracture height. Substitution of (13) into (12) yields:

$$\mathbf{U}\_{\mathbf{U}} = \left( \mathbf{-} \frac{k\_f \cdot k\_r}{n + 2} \frac{\partial \mathbf{u}^{\prime \mathbf{w}^{\prime \mathbf{w}}}}{\partial \mathbf{x}} \right)^{1/n} \tag{14}$$

**Figure 1.** Geometrical scheme of the PKN model

Physically significant speed of the front is neither zero, nor infinite. According to (14) and the SE (7) it is possible only when the function *y* =*w <sup>n</sup>*+2 is linear in *x*. Hence, in the consid‐ ered problem, the exponent *α* in the modified opening *y* =*w* 1/*α* is *α* =1 / (*n* + 2). Then the con‐ tinuity equation (1), the velocity equation (2), the initial condition (3), the BC (4), (5) and the SE (6) become, respectively,

$$
\frac{\partial \frac{\partial y}{\partial t}}{\partial t} + \upsilon \frac{\partial \frac{\partial y}{\partial x}}{\partial x} + \frac{y}{\alpha} \frac{\partial \upsilon}{\partial x} + \frac{y^{1-\alpha}}{\alpha} q\_l = 0 \tag{15}
$$

$$w = \left( -k\_f \, k\_r \alpha \frac{\partial \, y}{\partial x} \right)^{1/n} \tag{16}$$

Following [4], introduce the normalized variables *xd* = *x* / *xn*, *x*\**<sup>d</sup>* = *x*\* / *xn*, *td* =*t* / *tn*, *vd* =*v* / *v*, *v*\**<sup>d</sup>* =*v*\* / *vn*, *yd* = *y* / *yn*, *qd* =*q* / *qn*, *qld* =*ql* / *qln*, where the normalizing quantities are:

values of the time and influx per unit height, respectively. In terms of the normalized varia‐ bles, the problem (15)-(20) has the same form except that the multiplier *k <sup>f</sup> kr* is changed to the unit. This excludes the consistency factor, the elasticity modulus and the height. Here‐

where *β<sup>q</sup>* is a prescribed constant; for a constant influx, *β<sup>q</sup>* =0. In the case of zero leak-off, the

with *β<sup>w</sup>* = 1 + (*n* + 1)*β<sup>q</sup>* /(2*n* + 3), *β\** = 2(*n* + 1) + (*n* + 2)*β<sup>q</sup>* /(2*n* + 3). In (22), *ξ\** and *V\** are con‐

*y* =*Y* (*ξ*)*t*

*βw*/*α*

(*ξ*)=*o*((*ξ\** - *ξ*)*<sup>α</sup>*-1). For the exponent *β<sup>l</sup>* , it follows that *β<sup>l</sup>* =*β<sup>w</sup>* - 1.

Consider the case when the influx is prescribed by the power dependence in time:

*q*0 (*t*)=*t*

solution of (15)-(20) with the influx (21) may be found in self-similar variables:

stants, expressing respectively the self-similar fracture length and the front speed.

We prescribe the functions *Y* (*ξ*) and *V* (*ξ*) by power series in the variable *τ* =1 - *ξ* / *ξ\**:

(for zero leak-off, *q <sup>j</sup>* =0, *j* = 0,1,...). Then the coefficients *b <sup>j</sup>*

for *j* = 2, 3,… are:

( *j* - *k* + 1 + *αk*)*akb <sup>j</sup>*-*<sup>k</sup>* +1 + (*αj* -

recurrently from the equations (15), (16) re-written in self-similar variables. Omitting techni‐

, *V* (*ξ*)=*V*\*∑

*j*=0 ∞ *b j τ j*

(*ξ*) is represented as *Ql*

*βw β*\*

)*a <sup>j</sup>* - *Cl* ∑ *k*=1 *j*

*β*\* -1 *v*\*=*V*\* *t β*\* -1

, *vn* = *xn* / *tn*, *yn* =(*qntn* / *xn*)1/*α* with *tn* and *qn* being arbitrary typical

Modified Formulation, ε-Regularization and the Efficient Solution of Hydraulic Fracture Problems

*<sup>β</sup><sup>q</sup>* (21)

*q* =*Y* (*ξ*)*αV* (*ξ*)*t*

, where *Ql*(*ξ*) is a prescribed function, which may be

*<sup>β</sup><sup>q</sup>* ()

http://dx.doi.org/10.5772/56218

647

is also represented in the form

, (22)

and *a <sup>j</sup>*

*ckq <sup>j</sup>*-*<sup>k</sup>* } (23)

with known

are found

(*ξ*)=*τ <sup>α</sup>* ∑ *j*=0 *∞ q j τ j*

*xn* =(*k <sup>f</sup> krqn*

, (21)

*n*+2 *tn* 2*n*+2

*x* =*ξt β*\* *x*\*=*ξ*\* *t*

with separated variables: *ql* =*Ql*

singular at the front as *Ql*

coefficients *q <sup>j</sup>*

)1/(2*n*+3)

after in this section we use (15)-(20) assuming *k <sup>f</sup>* =1 and *kr* =1.

*<sup>β</sup>*\* *v* =*V* (*ξ*)*t*

We may account for leak-off by assuming that the term *ql*

*<sup>Y</sup>* (*ξ*)= *<sup>ξ</sup>*\*

where *V\** =*ξ\*β\**. The leak-off function *Ql*

cal details, the coefficients *b <sup>j</sup>*

*<sup>b</sup> <sup>j</sup>* <sup>=</sup> - <sup>1</sup>

*<sup>j</sup>* <sup>+</sup> *<sup>α</sup>* {∑ *k*=2 *j*

(*ξ*)*t βl*

*n*+1 *β*\* *n <sup>α</sup>* ∑ *j*=1 ∞ *a j τ j*

$$y(\mathbf{x}, \ t\_0) = y\_0(\mathbf{x}) \tag{17}$$

$$\left[y^{\alpha}\upsilon\right]\_{x=0} = q\_0(t),\tag{18}$$

$$y(\mathbf{x}\_\*) = 0\tag{19}$$

$$\upsilon\_{\ast}(t) = \frac{d\,\,\,\omega\_{\ast}}{dt} = \left[ \left( \cdot k\_{\,\,f} k\_{r} \alpha \frac{\partial \,\,y}{\partial \,\,x} \right)^{1/n} \right]\_{\,\,\,\pi\,\,x\_{\ast}} \tag{20}$$

Following [4], introduce the normalized variables *xd* = *x* / *xn*, *x*\**<sup>d</sup>* = *x*\* / *xn*, *td* =*t* / *tn*, *vd* =*v* / *v*, *v*\**<sup>d</sup>* =*v*\* / *vn*, *yd* = *y* / *yn*, *qd* =*q* / *qn*, *qld* =*ql* / *qln*, where the normalizing quantities are: *xn* =(*k <sup>f</sup> krqn n*+2 *tn* 2*n*+2 )1/(2*n*+3) , *vn* = *xn* / *tn*, *yn* =(*qntn* / *xn*)1/*α* with *tn* and *qn* being arbitrary typical values of the time and influx per unit height, respectively. In terms of the normalized varia‐ bles, the problem (15)-(20) has the same form except that the multiplier *k <sup>f</sup> kr* is changed to the unit. This excludes the consistency factor, the elasticity modulus and the height. Here‐ after in this section we use (15)-(20) assuming *k <sup>f</sup>* =1 and *kr* =1.

Consider the case when the influx is prescribed by the power dependence in time:

$$q\_0(t) = t^{\beta\_q} \tag{21}$$

, (21)

Physically significant speed of the front is neither zero, nor infinite. According to (14) and

ered problem, the exponent *α* in the modified opening *y* =*w* 1/*α* is *α* =1 / (*n* + 2). Then the con‐ tinuity equation (1), the velocity equation (2), the initial condition (3), the BC (4), (5) and the

∂ *y*

is linear in *x*. Hence, in the consid‐

*<sup>α</sup> ql* =0 (15)

<sup>∂</sup> *<sup>x</sup>* )1/*<sup>n</sup>* (16)

(*t*), (18)

(20)

*y*(*x*, *t*0) = *y*0(*x*) (17)

*y*(*x\**) =0 (19)

the SE (7) it is possible only when the function *y* =*w <sup>n</sup>*+2

∂ *y* <sup>∂</sup> *<sup>t</sup>* + *v* ∂ *y* <sup>∂</sup> *<sup>x</sup>* <sup>+</sup> *<sup>y</sup> α* ∂ *v* <sup>∂</sup> *<sup>x</sup>* <sup>+</sup> *<sup>y</sup>* 1-*<sup>α</sup>*

*v*\* (*t*)= *<sup>d</sup> <sup>x</sup>*\*

*v* =(-*k <sup>f</sup> krα*

*<sup>y</sup> <sup>α</sup><sup>v</sup> <sup>x</sup>*=0 <sup>=</sup>*q*<sup>0</sup>

*dt* = (-*k <sup>f</sup> krα*

∂ *y* <sup>∂</sup> *<sup>x</sup>* )1/*<sup>n</sup>*

*x*=*x*\*

SE (6) become, respectively,

**Figure 1.** Geometrical scheme of the PKN model

646 Effective and Sustainable Hydraulic Fracturing

where *β<sup>q</sup>* is a prescribed constant; for a constant influx, *β<sup>q</sup>* =0. In the case of zero leak-off, the solution of (15)-(20) with the influx (21) may be found in self-similar variables:

$$\mathbf{x} = \boldsymbol{\xi} \, t^{\beta\_{\ast}} \mathbf{x}\_{\ast} = \boldsymbol{\xi}\_{\ast} t^{\beta\_{\ast}} \boldsymbol{\upsilon} = V(\boldsymbol{\xi}) t^{\beta\_{\ast} \cdot 1} \boldsymbol{\upsilon}\_{\ast} = V\_{\ast} t^{\beta\_{\ast} \cdot 1} \boldsymbol{\ y} = Y(\boldsymbol{\xi}) t^{\beta\_{\ast}/\alpha} \boldsymbol{\eta} = Y(\boldsymbol{\xi})^{\alpha} V(\boldsymbol{\xi}) t^{\beta\_{\ast}} \tag{6}$$

with *β<sup>w</sup>* = 1 + (*n* + 1)*β<sup>q</sup>* /(2*n* + 3), *β\** = 2(*n* + 1) + (*n* + 2)*β<sup>q</sup>* /(2*n* + 3). In (22), *ξ\** and *V\** are con‐ stants, expressing respectively the self-similar fracture length and the front speed.

We may account for leak-off by assuming that the term *ql* is also represented in the form with separated variables: *ql* =*Ql* (*ξ*)*t βl* , where *Ql*(*ξ*) is a prescribed function, which may be singular at the front as *Ql* (*ξ*)=*o*((*ξ\** - *ξ*)*<sup>α</sup>*-1). For the exponent *β<sup>l</sup>* , it follows that *β<sup>l</sup>* =*β<sup>w</sup>* - 1.

We prescribe the functions *Y* (*ξ*) and *V* (*ξ*) by power series in the variable *τ* =1 - *ξ* / *ξ\**:

$$Y(\xi) = \frac{\xi^{n+1}\beta^{n}}{\alpha} \sum\_{j=1}^{\alpha} a\_{j}\pi^{-j}\prime\prime\prime\prime\prime\prime = V\*\sum\_{j=0}^{\alpha} b\_{j}\pi^{-j}\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\prime\tag{22}$$

where *V\** =*ξ\*β\**. The leak-off function *Ql* (*ξ*) is represented as *Ql* (*ξ*)=*τ <sup>α</sup>* ∑ *j*=0 *∞ q j τ j* with known coefficients *q <sup>j</sup>* (for zero leak-off, *q <sup>j</sup>* =0, *j* = 0,1,...). Then the coefficients *b <sup>j</sup>* and *a <sup>j</sup>* are found recurrently from the equations (15), (16) re-written in self-similar variables. Omitting techni‐ cal details, the coefficients *b <sup>j</sup>* for *j* = 2, 3,… are:

$$b\_{j}b\_{j} = -\frac{1}{j+\alpha} \left| \sum\_{k=2}^{j} \left[ (j-k+1+\alpha k)a\_{k}b\_{j+k+1} + \left(\alpha j - \frac{\beta\_{w}}{\beta\_{w}}\right)a\_{j} \right] - \mathbf{C}\_{l} \sum\_{k=1}^{j} c\_{k}q\_{j+k} \right| \tag{23}$$

with the starting values *a*<sup>1</sup> <sup>=</sup>*b*<sup>0</sup> =1, *b*<sup>1</sup> <sup>=</sup> <sup>1</sup> <sup>1</sup> <sup>+</sup> *<sup>α</sup>* ( - *α* + *βw <sup>β</sup>*\* + *Cl q*0 ), *c*<sup>1</sup> =1. In (24), *Cl* =( *<sup>α</sup> ξ*\* *n*+1 *β*\* *<sup>n</sup>*+1/*<sup>α</sup>* ) *<sup>α</sup>* and the coefficients *a <sup>j</sup>* and *c <sup>j</sup>* are found recurrently from equations: ∑ *k*=0 ∞ (*<sup>k</sup>* <sup>+</sup> 1)*ak* +1*<sup>τ</sup> <sup>k</sup>* =(<sup>∑</sup> *j*=0 ∞ *b j τ j* ) *n* , ∑ *k*=1 ∞ *ck<sup>τ</sup> <sup>k</sup>* <sup>=</sup>*τ*(<sup>∑</sup> *j*=0 ∞ *<sup>a</sup> <sup>j</sup>*+1*<sup>τ</sup> <sup>j</sup>* ) *α* . In particular, for the first five coefficients *a <sup>j</sup>* and *c <sup>j</sup>* we have:

$$\begin{aligned} a\_1 &= b\_0 = 1, \ a\_2 = \frac{1}{2} n b\_1, \ a\_3 = \frac{1}{6} n \left[ (n - 1) b\_1^2 + 2 b\_2 \right] \\ a\_4 &= \frac{1}{24} n \left[ (n - 1)(n - 2) b\_1^3 + 6(n - 1) b\_1 b\_2 + 6 b\_3 \right] \end{aligned} \tag{24}$$

faces [17]. As to the authors' knowledge, to this date, the only paper, in which the level set method has been applied to the hydraulic fractures, is that by Peirce and Detournay [20]. Since these authors used the conventional formulation not including the SE, special techni‐ que, called by them "implicit level set method", was suggested and successfully used. Still, direct employing of the SE looks superior in simplicity and possibility to use the standard

Modified Formulation, ε-Regularization and the Efficient Solution of Hydraulic Fracture Problems

From now on, we focus on other computational advantages of the modified formulation which evidently appear when considering the PKN model. Recall that this model is basic for a wide class of simulators employing the P3D models (for example, see [1]). Since the Nordgren's problem is 1D, it is convenient to use the spatial coordinate *x* normalized by the

1

 Va

a

,*l*


<sup>∂</sup> *<sup>ς</sup>* )1/*<sup>n</sup>* (28)

http://dx.doi.org/10.5772/56218

649

(*t*), (30)

*<sup>ς</sup>*=1 (32)

*y*(*ς*, *t*0) = *y*0(*ς*) (29)

*y*(1, *t*)=0 (31)

*x\**(*t*0) = *x\**<sup>0</sup> (33)

fracture length *x\**: *ς* = *x* / *x\**. In this spatial variable, the problem (15)-(20) reads:

aV

*<sup>v</sup>* <sup>=</sup> *<sup>k</sup> x\** 1/*<sup>n</sup>* (- ∂ *y*

*<sup>y</sup> <sup>α</sup><sup>v</sup> <sup>ς</sup>*=0 <sup>=</sup>*q*<sup>0</sup>

opening is zero. Thus we actually have the initial condition for the length:

where *k* =(*k <sup>f</sup> krα*)1/*n*. When writing (27)-(32), we have used the symbols *y* and *v* for the func‐ tions *y*˜(*ς*, *t*)= *y*(*x*(*ς*), *t*) and *v*˜(*ς*, *t*)=*v*(*x*(*ς*), *t*) omitting the tilde. Note that the initial condi‐ tion (28) defines also the initial value *x\**0 of the fracture length as the end point where the

In the considered problem the dependence (28) between the unknowns is explicit. Therefore, we may substitute (28) into (27). This yields the PDE of the first order in time and of the sec‐

*v*\* (*t*)= *<sup>d</sup> <sup>x</sup>*\* *dt* <sup>=</sup> *<sup>k</sup> x*\* 1/*<sup>n</sup>* (- ∂ *y* <sup>∂</sup> *<sup>ς</sup>* )1/*<sup>n</sup>*

(*ς*)=0, we have *x\**<sup>0</sup> =0.

\* \* \*

*y y <sup>v</sup> v v y y <sup>q</sup> tx x*

V

well-established technique.

In the case when *y*<sup>0</sup>

ond order in ς:

$$a\_5 = \frac{1}{120}n\left[ (n-1)(n-2)(n-3)b\_1^4 + 12(n-1)(n-2)b\_1^2b\_2 + 24(n-1)b\_1b\_3 + 24b\_4 \right] \tag{25}$$

$$c\_1 = a\_1 = 1, \quad c\_2 = \begin{bmatrix} 1 \ -\alpha \end{bmatrix} a\_2, \quad c\_3 = \frac{1}{2} (1 \ -\alpha) \begin{bmatrix} -\alpha a\_2^2 + 2a\_3 \end{bmatrix}, \quad c\_4 = \frac{1}{6} (1 \ -\alpha) \begin{bmatrix} \alpha(\alpha + 1) a\_2^3 \ -6\alpha a\_2 a\_3 + 6a\_4 \end{bmatrix} \tag{26}$$
 
$$c\_5 = \frac{1}{24} (1 \ -\alpha) \begin{bmatrix} -\alpha(\alpha + 1)(\alpha + 2) a\_2^4 + 12\alpha(\alpha + 1) a\_2^2 a\_3 \ -24\alpha a\_2 a\_4 + 24a\_5 \end{bmatrix} . \tag{27}$$

Starting from *a*<sup>1</sup> =*b*<sup>0</sup> =*c*<sup>1</sup> =1, *b*<sup>1</sup> =( - *α* + *β<sup>w</sup>* / *β\** + *Cl q*0) /(1 + *α*), we find *a*2 from the second of (25) and *c*<sup>2</sup> from the second of (26). Then (24) gives *b*2, the third of (25) defines *a*3, the third of (26) defines *c*3 and so on.

The value *ξ\** of the self-similar fracture length is found from the self-similar BC at the inlet: *<sup>Y</sup>* (0) *<sup>α</sup><sup>V</sup>* (0)=1. By using the solution (23)-(26) for various *ξ\**, we find the one, which meets the BC with a prescribed tolerance.

In the case of a *Newtonian fluid* (*n* = 1, *α* =1 / 3) equations (23)-(26) extend the analytical solu‐ tion, obtained in [4], to the case of non-zero leak-off. In the case of a *perfectly plastic fluid* (*n* = 0, *α* =1 / 2), we have *ak* =0, *ck* =0 for *k* > 1. Then for constant influx (*β<sup>q</sup>* =0), the solution is: *<sup>ξ</sup>\** =(9 / 8)1/3 , *Y* (*ξ*) =2(*ξ\** - *ξ*), *V* (*ξ*) =*V\** 1 + 1 *β\** 2*ξ\** ∑ *j*=1 *∞* 2 <sup>2</sup> *<sup>j</sup>* <sup>+</sup> <sup>1</sup> *<sup>q</sup> <sup>j</sup>*-1*<sup>τ</sup> <sup>j</sup>* , *V\** =2 / <sup>3</sup>*ξ\**. Calculations show that the difference between the self-similar solutions for these two limiting cases is quite

small. For instance, for zero leak-off, the self-similar fracture length is *ξ\*p* =1.04004 for a per‐ fectly plastic fluid; it is *ξ\*N* =1.00101 for a Newtonian fluid. In both cases, the particle veloci‐ ty is almost constant along the fracture. Thinning fluids (0 < *n* < 1), being intermediate between Newtonian and perfectly plastic fluids, the conclusions hold for an arbitrary thin‐ ning fluid. A detailed discussion of the solution is given in [19].

#### **4. Increasing efficiency of numerical simulations**

As mentioned, in the general case, using the SE (7) opens options for tracing the fracture in the 3D space by level set, fast marching and other methods of the theory of propagating sur‐ faces [17]. As to the authors' knowledge, to this date, the only paper, in which the level set method has been applied to the hydraulic fractures, is that by Peirce and Detournay [20]. Since these authors used the conventional formulation not including the SE, special techni‐ que, called by them "implicit level set method", was suggested and successfully used. Still, direct employing of the SE looks superior in simplicity and possibility to use the standard well-established technique.

From now on, we focus on other computational advantages of the modified formulation which evidently appear when considering the PKN model. Recall that this model is basic for a wide class of simulators employing the P3D models (for example, see [1]). Since the Nordgren's problem is 1D, it is convenient to use the spatial coordinate *x* normalized by the fracture length *x\**: *ς* = *x* / *x\**. In this spatial variable, the problem (15)-(20) reads:

$$\frac{\partial \hat{y}}{\partial t} = -\frac{y}{\alpha \infty} \frac{\partial \upsilon}{\partial \underline{\varsigma}} + \frac{\underline{\varsigma} \upsilon\_\* - \upsilon}{\upsilon\_\*} \frac{\partial \underline{y}}{\partial \underline{\varsigma}} - \frac{y^{1-\alpha}}{\alpha} q\_{1\prime} \tag{27}$$

$$
\omega \upsilon = \frac{k}{\chi\_\*^{1/n}} \left( -\frac{\partial \cdot y}{\partial \cdot \varsigma} \right)^{1/n} \tag{28}
$$

$$y\_0(\varepsilon\_\prime \ t\_0) = y\_0(\varepsilon) \tag{29}$$

$$\left[\underline{\mathfrak{z}}\,^{\alpha}\underline{v}\,\right]\_{\prec 0} = q\_0(t),\tag{30}$$

$$y(1, t) = 0\tag{31}$$

$$\text{Cov}\_\*(t) = \frac{d \,\, \text{x}\_\*}{dt} = \frac{k}{\chi\_\*^{1/n}} \mathbb{E}\left[ \left( -\frac{\partial \,\, \text{y}}{\partial \,\, \zeta} \right)^{1/n} \right]\_{\zeta = \mathbf{1}} \tag{32}$$

where *k* =(*k <sup>f</sup> krα*)1/*n*. When writing (27)-(32), we have used the symbols *y* and *v* for the func‐ tions *y*˜(*ς*, *t*)= *y*(*x*(*ς*), *t*) and *v*˜(*ς*, *t*)=*v*(*x*(*ς*), *t*) omitting the tilde. Note that the initial condi‐ tion (28) defines also the initial value *x\**0 of the fracture length as the end point where the opening is zero. Thus we actually have the initial condition for the length:

$$\mathbf{x}\_{\ast}(t\_0) = \mathbf{x}\_{\ast 0} \tag{33}$$

In the case when *y*<sup>0</sup> (*ς*)=0, we have *x\**<sup>0</sup> =0.

with the starting values *a*<sup>1</sup> <sup>=</sup>*b*<sup>0</sup> =1, *b*<sup>1</sup> <sup>=</sup> <sup>1</sup>

648 Effective and Sustainable Hydraulic Fracturing

and *c <sup>j</sup>*

*<sup>a</sup>*<sup>1</sup> <sup>=</sup>*b*<sup>0</sup> =1, *<sup>a</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup>

<sup>24</sup> *n* (*n* - 1)(*n* - 2)*b*<sup>1</sup>

<sup>2</sup> (1 - *α*) -α*a*<sup>2</sup>

<sup>24</sup> (1 - *α*) -α(α + 1)(α + 2)a2

*<sup>a</sup>*<sup>4</sup> <sup>=</sup> <sup>1</sup>

<sup>120</sup> *n* (*n* - 1)(*n* - 2)(*n* - 3)*b*<sup>1</sup>

Starting from *a*<sup>1</sup> =*b*<sup>0</sup> =*c*<sup>1</sup> =1, *b*<sup>1</sup> =( - *α* + *β<sup>w</sup>* / *β\** + *Cl*

, *Y* (*ξ*) =2(*ξ\** - *ξ*), *V* (*ξ*) =*V\** 1 +

ning fluid. A detailed discussion of the solution is given in [19].

**4. Increasing efficiency of numerical simulations**

the coefficients *a <sup>j</sup>*

*j*=0 ∞ *<sup>a</sup> <sup>j</sup>*+1*<sup>τ</sup> <sup>j</sup>* ) *α*

*<sup>a</sup>*<sup>5</sup> <sup>=</sup> <sup>1</sup>

defines *c*3 and so on.

*<sup>ξ</sup>\** =(9 / 8)1/3

*<sup>c</sup>*<sup>1</sup> <sup>=</sup>*a*<sup>1</sup> =1, *<sup>c</sup>*<sup>2</sup> =(1 - <sup>α</sup>)*a*2, *<sup>c</sup>*<sup>3</sup> <sup>=</sup> <sup>1</sup>

*<sup>c</sup>*<sup>5</sup> <sup>=</sup> <sup>1</sup>

the BC with a prescribed tolerance.

*ck<sup>τ</sup> <sup>k</sup>* <sup>=</sup>*τ*(<sup>∑</sup>

∑ *k*=1 ∞

<sup>1</sup> <sup>+</sup> *<sup>α</sup>* ( - *α* +

are found recurrently from equations: ∑

. In particular, for the first five coefficients *a <sup>j</sup>*

<sup>4</sup> <sup>+</sup> 12(*<sup>n</sup>* - 1)(*<sup>n</sup>* - 2)*b*<sup>1</sup>

<sup>2</sup> <sup>+</sup> <sup>2</sup>*a*<sup>3</sup> , *<sup>c</sup>*<sup>4</sup> <sup>=</sup> <sup>1</sup>

and *c*<sup>2</sup> from the second of (26). Then (24) gives *b*2, the third of (25) defines *a*3, the third of (26)

The value *ξ\** of the self-similar fracture length is found from the self-similar BC at the inlet: *<sup>Y</sup>* (0) *<sup>α</sup><sup>V</sup>* (0)=1. By using the solution (23)-(26) for various *ξ\**, we find the one, which meets

In the case of a *Newtonian fluid* (*n* = 1, *α* =1 / 3) equations (23)-(26) extend the analytical solu‐ tion, obtained in [4], to the case of non-zero leak-off. In the case of a *perfectly plastic fluid* (*n* = 0, *α* =1 / 2), we have *ak* =0, *ck* =0 for *k* > 1. Then for constant influx (*β<sup>q</sup>* =0), the solution is:

that the difference between the self-similar solutions for these two limiting cases is quite small. For instance, for zero leak-off, the self-similar fracture length is *ξ\*p* =1.04004 for a per‐ fectly plastic fluid; it is *ξ\*N* =1.00101 for a Newtonian fluid. In both cases, the particle veloci‐ ty is almost constant along the fracture. Thinning fluids (0 < *n* < 1), being intermediate between Newtonian and perfectly plastic fluids, the conclusions hold for an arbitrary thin‐

As mentioned, in the general case, using the SE (7) opens options for tracing the fracture in the 3D space by level set, fast marching and other methods of the theory of propagating sur‐

1 *β\** 2*ξ\** ∑ *j*=1 *∞* 2

<sup>4</sup> <sup>+</sup> <sup>12</sup>α(<sup>α</sup> <sup>+</sup> 1)*a*<sup>2</sup>

<sup>2</sup> *nb*1, *<sup>a</sup>*<sup>3</sup> <sup>=</sup> <sup>1</sup>

*βw <sup>β</sup>*\* + *Cl q*0

<sup>6</sup> *n* (*n* - 1)*b*<sup>1</sup>

<sup>3</sup> <sup>+</sup> 6(*<sup>n</sup>* - 1)*b*1*b*<sup>2</sup> <sup>+</sup> <sup>6</sup>*b*<sup>3</sup>

2

2

<sup>2</sup> <sup>+</sup> <sup>2</sup>*b*<sup>2</sup>

<sup>6</sup> (1 - *α*) α(*α* + 1)*a*<sup>2</sup>

*a*<sup>3</sup> - 24α*a*2*a*<sup>4</sup> + 24*a*<sup>5</sup> .

*q*0) /(1 + *α*), we find *a*2 from the second of (25)

<sup>2</sup> *<sup>j</sup>* <sup>+</sup> <sup>1</sup> *<sup>q</sup> <sup>j</sup>*-1*<sup>τ</sup> <sup>j</sup>* , *V\** =2 / <sup>3</sup>*ξ\**. Calculations show

), *c*<sup>1</sup> =1. In (24), *Cl* =( *<sup>α</sup>*

and *c <sup>j</sup>*

*b*<sup>2</sup> + 24(*n* - 1)*b*1*b*<sup>3</sup> + 24*b*<sup>4</sup> (25)

<sup>3</sup> - <sup>6</sup>α*a*2*a*<sup>3</sup> <sup>+</sup> <sup>6</sup>*a*<sup>4</sup> ,

*k*=0 ∞

*ξ*\* *n*+1 *β*\* *<sup>n</sup>*+1/*<sup>α</sup>* )

we have:

*j*=0 ∞ *b j τ j* ) *n* ,

(*<sup>k</sup>* <sup>+</sup> 1)*ak* +1*<sup>τ</sup> <sup>k</sup>* =(<sup>∑</sup>

*<sup>α</sup>* and

(24)

(26)

In the considered problem the dependence (28) between the unknowns is explicit. Therefore, we may substitute (28) into (27). This yields the PDE of the first order in time and of the sec‐ ond order in ς:

$$\frac{\partial \, y}{\partial t} = \frac{k}{z\_{\cdot}} \left\{ \frac{y}{\alpha n} \left( - \frac{\partial \, y}{\partial \cdot \zeta} \right)^{\left(1 \cdot n\right)/n} \frac{\partial^2 \, y}{\partial \cdot \zeta^2} + \left[ \zeta \left( - \frac{\partial \, y}{\partial \cdot \zeta} \right)\_{\zeta = 1}^{1/n} - \left( - \frac{\partial \, y}{\partial \cdot \zeta} \right)^{\left(1 \cdot n\right)} \right] \frac{\partial \, y}{\partial \cdot \zeta} \right\} - \frac{y^{1 \cdot n}}{\alpha} q\_{\vert} \tag{34}$$

where *z\** = *x\** (1+*n*)/*n* is the modified fracture length. For it, the SE (32) becomes:

$$\frac{d^2z\_\*}{dt^2} = \frac{n+1}{n}k\left[\left(-\frac{\partial \, \, \underline{y}}{\partial \, \underline{z}}\right)^{1/n}\right]\_{\zeta=1} \tag{35}$$

The results are summarized in Figure 2. The thick solid line presents the Nordgren's curve (Fig. 2 of reference [9]), obtained in the time range [0.01, 5.0]. The author did not comment on the accuracy of his calculations. The accuracy, as stated in [3], is certainly below 1% which is also evident from the fact that the Nordgren's curve intersects the asymptotic dot‐ ted line, corresponding to small time, and the asymptotic dashed line, corresponding to

Modified Formulation, ε-Regularization and the Efficient Solution of Hydraulic Fracture Problems

The results obtained by employing the computational scheme described are presented in Figure 2 by the thin solid line. It can be seen that the solution starting from the small time asymptotics tends to that corresponding to the large time one. The calculations are per‐

instability are noted. Even though the calculations were performed in the MATLAB environ‐ ment with use of the standard ODE solver *ode15s*, the run time to cover the time range of 12 orders is 10 seconds. The said confirms high efficiency of the computational scheme suggest‐

3 Previously, Kovalyshen and Detournay [21] also have solved the same problem in the range of time wider than that covered by Nordgren. The authors performed calculations starting from the small-time asymptote as an initial condition at τ =10-8 and presented numerical results in the range 10-5 < τ < 500. Their results are shown by circles in Figure 2. They are indistinguishable from our numerical solution. The authors do not discuss the accuracy, stability and robustness of their calculations. From our results, to which the relative error does not exceed 10-3, we may conclude that the relative error of the results, given for the fracture length in Table 1 of the paper [21], does not exceed 1% in the whole range of

. Specially designed numeri‐

http://dx.doi.org/10.5772/56218

651

. No signs of

formed by using *ε*-regularized BC (36) and SE (37) with *ε* =10-4

**Figure 2.** Dimensionless fracture length vs dimensionless time

cal experiments have shown that the relative error does not exceed 10−<sup>3</sup>

large time.

ed. 3

the calculations.

From (34) and (35) it follows that after spatial discretization of these equations and the BC (30), (31), we obtain a well-posed system of ordinary differential equations (ODE) of the first order in time to be solved under the initial conditions (29) and (33). Actually, the problem does not require regularization. Still, as noted in [6], not to have too stiff system of ODE, it is reasonable to employ the *ε*-regularized forms of the BC (31) and the SE (34):

$$\int y\left(1-\varepsilon, t\right) = \left[\frac{n}{k(n+1)}\frac{d\mathbf{z}\_\*}{dt}\right]^n \,\mathrm{s}\tag{36}$$

$$\frac{\partial^d z\_\*}{\partial t} = \frac{n+1}{n} k \left[ \left( -\frac{\partial \cdot y}{\partial \cdot \zeta} \right)^{1/n} \right]\_{\zeta = 1 \cdot \varepsilon} \tag{37}$$

where *ε* is a small positive value in the range from 10-3 to 10-6 (it may yet be less depending on round-off errors of computer calculations).

Solving the system of ODE, resulting from spatial discretization of (34), (35), under the ini‐ tial conditions (30), (33) may be performed by well-developed methods, like the Runge-Kut‐ ta method. Standard solvers are of immediate use. Emphasize that this option has appeared only due to employing the *local* SE rather than the global mass balance for tracing the front propagation. As show numerous numerical experiments, summarized in [6] for a Newtoni‐ an fluid (*n* = 1), this computational scheme provides highly accurate and stable results with small time expense.

As an example, we use this scheme to extend the Nordgren's numerical results [9] on the dependence of the fracture length on time for a Newtonian fluid and Carter's leak-off. To this end, the normalizing length *xn* and normalizing time *tn* are taken as those in the paper by Nordgren:

$$\chi\_n = \pi \left(\frac{\mu \, q\_0^{5} \hbar}{128 E \, ^\circ \text{C}^{-8}}\right)^{1/3}, \; t\_n = \pi \, ^2 \left(\frac{\mu \, q\_0^{2} \hbar}{16 E \, ^\circ \text{C}^{-5}}\right)^{2/3}.$$

Herein, *μ* is the dynamic viscosity, *C* is the fluid-loss coefficient in the Carter's leak-off term *ql* =2*C* / *t* - *τ*, *τ* is the time at which the fracture front reaches a point *x*. In accordance with (18), we have used the influx *q*0 per unit height, while Nordgren wrote his normalizing val‐ ues in terms of the total influx *q*1 through the entire height *h* (*q*<sup>1</sup> =*q*0*h* ).

The results are summarized in Figure 2. The thick solid line presents the Nordgren's curve (Fig. 2 of reference [9]), obtained in the time range [0.01, 5.0]. The author did not comment on the accuracy of his calculations. The accuracy, as stated in [3], is certainly below 1% which is also evident from the fact that the Nordgren's curve intersects the asymptotic dot‐ ted line, corresponding to small time, and the asymptotic dashed line, corresponding to large time.

The results obtained by employing the computational scheme described are presented in Figure 2 by the thin solid line. It can be seen that the solution starting from the small time asymptotics tends to that corresponding to the large time one. The calculations are per‐ formed by using *ε*-regularized BC (36) and SE (37) with *ε* =10-4 . Specially designed numeri‐ cal experiments have shown that the relative error does not exceed 10−<sup>3</sup> . No signs of instability are noted. Even though the calculations were performed in the MATLAB environ‐ ment with use of the standard ODE solver *ode15s*, the run time to cover the time range of 12 orders is 10 seconds. The said confirms high efficiency of the computational scheme suggest‐ ed. 3

**Figure 2.** Dimensionless fracture length vs dimensionless time

∂ *y* <sup>∂</sup> *<sup>t</sup>* <sup>=</sup> *<sup>k</sup> z*\* { *y <sup>α</sup><sup>n</sup>* (- ∂ *y* <sup>∂</sup> *<sup>ς</sup>* )

650 Effective and Sustainable Hydraulic Fracturing

where *z\** = *x\**

(1-*n*)/*n* ∂<sup>2</sup> *y*

*d z\* dt* <sup>=</sup> *<sup>n</sup>* <sup>+</sup> <sup>1</sup>

*d z*\* *dt* <sup>=</sup> *<sup>n</sup>* <sup>+</sup> <sup>1</sup>

on round-off errors of computer calculations).

, *tn* <sup>=</sup>*<sup>π</sup>* 2( *μq*<sup>0</sup>

2 *h* 16*E* ' *C* 5 ) 2/3 .

ues in terms of the total influx *q*1 through the entire height *h* (*q*<sup>1</sup> =*q*0*h* ).

small time expense.

by Nordgren:

5 *h* 128*E* ' *C* 8 ) 1/3

*xn* <sup>=</sup>*π*( *μq*<sup>0</sup>

<sup>∂</sup> *<sup>ς</sup>* <sup>2</sup> <sup>+</sup> *<sup>ς</sup>*(-

∂ *y* <sup>∂</sup> *<sup>ς</sup>* )*ς*=1 1/*<sup>n</sup>* - (- ∂ *y* <sup>∂</sup> *<sup>ς</sup>* )1/*<sup>n</sup>* <sup>∂</sup> *<sup>y</sup>*

(1+*n*)/*n* is the modified fracture length. For it, the SE (32) becomes:

∂ *y* <sup>∂</sup> *<sup>ς</sup>* )1/*<sup>n</sup>*

From (34) and (35) it follows that after spatial discretization of these equations and the BC (30), (31), we obtain a well-posed system of ordinary differential equations (ODE) of the first order in time to be solved under the initial conditions (29) and (33). Actually, the problem does not require regularization. Still, as noted in [6], not to have too stiff system of ODE, it is

*n dz*

*n*

 e

*<sup>n</sup> k* (-

reasonable to employ the *ε*-regularized forms of the BC (31) and the SE (34):

( ) \* 1 , ( 1)

é ù - = ê ú

∂ *y* <sup>∂</sup> *<sup>ς</sup>* )1/*<sup>n</sup>*

where *ε* is a small positive value in the range from 10-3 to 10-6 (it may yet be less depending

Solving the system of ODE, resulting from spatial discretization of (34), (35), under the ini‐ tial conditions (30), (33) may be performed by well-developed methods, like the Runge-Kut‐ ta method. Standard solvers are of immediate use. Emphasize that this option has appeared only due to employing the *local* SE rather than the global mass balance for tracing the front propagation. As show numerous numerical experiments, summarized in [6] for a Newtoni‐ an fluid (*n* = 1), this computational scheme provides highly accurate and stable results with

As an example, we use this scheme to extend the Nordgren's numerical results [9] on the dependence of the fracture length on time for a Newtonian fluid and Carter's leak-off. To this end, the normalizing length *xn* and normalizing time *tn* are taken as those in the paper

Herein, *μ* is the dynamic viscosity, *C* is the fluid-loss coefficient in the Carter's leak-off term *ql* =2*C* / *t* - *τ*, *τ* is the time at which the fracture front reaches a point *x*. In accordance with (18), we have used the influx *q*0 per unit height, while Nordgren wrote his normalizing val‐

*y t k n dt* e

*<sup>n</sup> k* (-

<sup>∂</sup> *<sup>ς</sup>* } - *<sup>y</sup>* 1-*<sup>α</sup>*

ë û <sup>+</sup> (36)

*<sup>ς</sup>*=1-*<sup>ε</sup>* (37)

*<sup>ς</sup>*=1 (35)

*<sup>α</sup> ql* (34)

<sup>3</sup> Previously, Kovalyshen and Detournay [21] also have solved the same problem in the range of time wider than that covered by Nordgren. The authors performed calculations starting from the small-time asymptote as an initial condition at τ =10-8 and presented numerical results in the range 10-5 < τ < 500. Their results are shown by circles in Figure 2. They are indistinguishable from our numerical solution. The authors do not discuss the accuracy, stability and robustness of their calculations. From our results, to which the relative error does not exceed 10-3, we may conclude that the relative error of the results, given for the fracture length in Table 1 of the paper [21], does not exceed 1% in the whole range of the calculations.

A similarly efficient computational scheme consists of solving the PDE (27), (28) with *ε*regularized BC and SE by the Crank-Nicolson method [3,6].

(iii) The dependence (12) between the particle velocity and the pressure, after averaging the

Modified Formulation, ε-Regularization and the Efficient Solution of Hydraulic Fracture Problems

where for sufficiently small *wav*, we have *Fv*(*wav*) =1. Since *wav* →0 when approaching the

(iv) The continuity equation is integrated over the cross-sectional height. As a result, the to‐ tal flux through a cross-section and the total leak-off replace, respectively, the flux and leakoff per unit height. The average velocity *vav*, being defined as the ratio of the total flux to the area *A* of the cross-section, the area *A* replaces the opening present in the PKN model. To preserve connection with the PKN model, we may use the flux, leak-off losses and area div‐

*<sup>h</sup> <sup>r</sup>* <sup>=</sup>*wvav*, *ql* <sup>=</sup> *Ql*

we have the modified continuity equation (15) in the unchanged form. The BC (18) at the inlet is also the same when denoting *q*0(*t*) =*Q*0(*t*) / *hr*, where *Q*0(*t*) is the prescribed total in‐

Take into account that by the definition of the average opening, we have *A*=*wavh <sup>f</sup>* , where *h <sup>f</sup>* is the fracture height in a considered cross-section. Thus, by the first of (42), *w* = *Ah <sup>f</sup>* / *hr*. Then, in view of (39) and (38), we may use the argument *w* instead of *wav* in the equations

*Gv*(*w*), and *v*\*=(-*<sup>k</sup> <sup>f</sup> <sup>w</sup> <sup>n</sup>*+1 <sup>∂</sup> *<sup>p</sup>*

where to simplify notation, we have omitted the subscript in the averaged velocity *vav*. The functions *Gp*(*w*) and *Gv*(*w*) are found in advance through the functions *F <sup>p</sup>*(*wav*) and *Fv*(*wav*).

Finally, by introducing the variables *<sup>y</sup>* <sup>=</sup>*<sup>w</sup>* 1/*α* and *<sup>ς</sup>* <sup>=</sup> *<sup>x</sup>* / *x\**, we arrive at the system similar to

 Va

\* \* \*

*y y <sup>v</sup> v v y y <sup>q</sup> tx x*

V

aV

1

a

,*l*


*h r*

<sup>∂</sup> *<sup>x</sup>* )1/*<sup>n</sup> x*=*x*\* ,

*n*+1 ∂ *p* <sup>∂</sup> *<sup>x</sup>* )1/*<sup>n</sup>*

*Fv*(*wav*) (40)

http://dx.doi.org/10.5772/56218

653

*<sup>x</sup>*=*x*\* (41)

, (42)

*n*+1 ∂ *p* <sup>∂</sup> *<sup>x</sup>* )1/*<sup>n</sup>*

velocity over the height of a cross-section, obtains a factor *Fv*(*wav*):

fracture front, the SE becomes:

flux.

(27)-(33):

*vav* =(-*k <sup>f</sup> wav*

*v*\*=(-*k <sup>f</sup> wav*

ided by a fixed reference height *hr*, say the height of the pay-layer. Then denoting

, *<sup>q</sup>* <sup>=</sup> *Avav*

*<sup>w</sup>* <sup>=</sup> *<sup>A</sup> h r*

(38), (40) and(41) writing them, respectively, as

<sup>∂</sup> *<sup>x</sup>* )1/*<sup>n</sup>*

They are such that *Gp*(*w*)=1 and *Gv*(*w*)=1 for sufficiently small *w*.

*<sup>p</sup>* <sup>=</sup>*krwGp*(*w*), *<sup>v</sup>* =(-*<sup>k</sup> <sup>f</sup> <sup>w</sup> <sup>n</sup>*+1 <sup>∂</sup> *<sup>p</sup>*

*Comment.* In view of the accepted linear dependence (13) between the net-pressure and the opening, the pressure may replace the opening as an unknown variable. Actually, under the normalization, yielding *kr* =1, there is no difference between *p* and *w*.
