Open access

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

Written By

Submitted: 19 July 2012 Published: 17 May 2013

DOI: 10.5772/56218

From the Proceeding

## Effective and Sustainable Hydraulic Fracturing

Edited by Andrew P. Bunger, John McLennan and Rob Jeffrey

Chapter metrics overview

View Full Metrics

## 1. Introduction

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 formulation 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 reasons. Numerical simulators, based on the conventional formulation (for example, see [7,8]), employ quite rough meshes, which themselves serve as specific ‘regularizators’. On the other 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 boundary 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 form

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.

, (ii) comprehension of its significance for proper numerical simulation of hydraulic fractures and (iii) distinguishing the particle velocity as a preferable variable

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.

. It 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 formulation of the problem has changed to the modified formulation, which opens new analytical and computational options for solving hydraulic fracture problems.

Below we employ these options. Section 2 contains a concise review of the modified formulation. In Section 3, we illustrate its analytical advantages by simple solutions for the Perkins-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-dimensional (P3D) models [1], the importance of which grows nowadays because of their employment 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=w1/α, where α is the exponent defined by the asymptotic behavior of the opening w at the fluid front (w=Ctrα, r is the distance from the front, t is the time). Then the continuity equation for a flow in a narrow channel reads:

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:

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)gradp=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 general case, keeping the velocity, which is a smooth function near the fluid front, as unknown is more convenient than to replace it by gradp, which is commonly strongly singular. 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:

yx,t0=y0(x)E3

The elliptic operator requires only one boundary condition (BC) on the fluid front Le. For instance, when neglecting the lag between the fluid front and the fracture contour, the BC may be the condition of the prescribed normal component qnx,t=yαvn of the flux:

yαvnx,t=q0x,t,                    xLeE4

where q₀(x) is a known function at Le; 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 Lf, we have w = 0 and equation (2) implies qfx,t=wv=0. Then (4) becomes:

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

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

v*=dxn*dt=vn*E6

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 fundamental for moving boundaries. Using (2) in (6) specifies the SE for a flow of incompressible fluid in a narrow channel [2,3]:

v*=dxn*dt=-1w(x*)D(w,p)pnx=x*E7

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 problems 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 computational point of view.

The speed equation (6) also yields important implications for numerical simulation of hydraulic 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 overdetermined 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:

p0pe1wD(w,p)dp=v*rεE8

where p0 is the pressure at the front, pε  is the pressure at a distance rε  from the front. Equation (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 assumed to be met at the distance rε to an accepted accuracy:

v*t=dxn*dt=-1wD(w,p)pnr=rεE9

The ε-regularized BC (8) allows one to avoid unfavorable computational effects. The ε-regularized 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 hydraulic fracturing. To find it, the fluid equations are complimented with solid mechanics and fracture equations.

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

Ayα=pE10

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:

KI=KIc,         KII=0E11

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 fracture) front. The conditions (11) define the possibility and the direction of the fracture propagation. The ε-regularized SE (9) serves to find the front propagation on the time steps.

## 3. Analytical solutions

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 modified 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γ˙n, 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:

v=-kfwn+1px1/nE12

where kf=1/(θM), for a plane channel, θ=2[2n+1]n (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 plane-strain conditions occur in cross-sections parallel to the fracture front. Then the elasticity equation (10) takes the form [9]:

p=krwE13

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:

v=-kfkrn+2wn+2x1/nE14

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=wn+2 is linear in x. Hence, in the considered problem, the exponent α in the modified opening y=w1/α is α=1/(n+2). Then the continuity equation (1), the velocity equation (2), the initial condition (3), the BC (4), (5) and the SE (6) become, respectively,

yt+vyx+yαvx+y1-ααql=0E15
v=-kfkrαyx1/nE16
yx,t0=y0(x)E17
[yαv]x=0=q0t, E18
y(x*)=0E19
v*t=dx*dt=-kfkrαyx1/nx=x*E20

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

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

q0t=tβqE21
, (21)

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

x=ξtβ*x*=ξ*tβ* v=Vξtβ*-1v*=V*tβ*-1 y=Y(ξ)tβw/αq=Y(ξ)αV(ξ)tβq

with βw=1+n+1βq/(2n+3), β*=2(n+1)+n+2βq/(2n+3). In (22), ξ* and V* are constants, 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(ξ*-ξ)α-1. For the exponent βl ,  it follows that βl=βw-1.

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

Yξ=ξ*n+1β*nαj=1ajτj,Vξ=V*j=0bjτj,E22

where V*=ξ*β*. The leak-off function Qlξ is represented as Qlξ=ταj=0qjτj with known coefficients qj (for zero leak-off, qj=0, j = 0,1,...). Then the coefficients bj and aj are found recurrently from the equations (15), (16) re-written in self-similar variables. Omitting technical details, the coefficients bj for j = 2, 3,… are:

bj=-1j+αk=2jj-k+1+αkakbj-k+1+(αj-βwβ*)aj-Clk=1jckqj-kE23

with the starting values a1=b0=1, b1=11+α(-α+βwβ*+Clq0), c1=1. In (24), Cl=αξ*n+1β*n+1/αα and the coefficients aj and cj are found recurrently from equations: k=0k+1ak+1τk =j=0bjτj n, k=1ckτk =τj=0aj+1τj α. In particular, for the first five coefficients aj and cj we have:

a1=b0=1,a2=12nb1,a3=16n[n-1b12+2b2]a4=124n[n-1n-2b13+6n-1b1b2+6b3]E24
a5=1120n[n-1n-2n-3b14+12n-1n-2b12b2+24n-1b1b3+24b4]E25
c1=a1=1, c2=1-αa2, c3=121-α-αa22+2a3,   c4=161-ααα+1a23-6αa2a3+6a4, c5=1241-α-αα+1α+2a24+12αα+1a22a3-24αa2a4+24a5.E26

Starting from a1=b0=c1=1, b1=(-α+βw/β*+Clq0)/(1+α), we find a2 from the second of (25) and c2 from the second of (26). Then (24) gives b2, the third of (25) defines a3, the third of (26) defines c3 and so on.

The value ξ* of the self-similar fracture length is found from the self-similar BC at the inlet: [Y(0)]αV0=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 solution, 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 (βq=0), the solution is: ξ*=(9/8)1/3, Y(ξ)=2(ξ*-ξ),  V(ξ)=V*1+1β*2ξ*j=122j+1qj-1τj, V*=2/3ξ*. 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 perfectly plastic fluid; it is ξ*N=1.00101 for a Newtonian fluid. In both cases, the particle velocity 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 thinning 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 surfaces [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 technique, 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:

yt=yαx*vς+ςv*vx*yςy1ααql,E27
v=kx*1/n-yς1/nE28
yς,t0=y0(ς)E29
[yαv]ς=0=q0t,   E30
y(1,t)=0E31
v*t=dx*dt=kx*1/n-yς1/nς=1E32

where k=kfkrα1/n. When writing (27)-(32), we have used the symbols y and v for the functions y˜(ς,t)=y(x(ς),t) and v˜(ς,t)=v(x(ς),t) omitting the tilde. Note that the initial condition (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:

x*t0=x*0E33

In the case when y0ς=0, we have x*0=0.

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 second order in ς:

yt=kz*yαn-yς(1-n)/n2yς2+ς-yςς=11/n--yς1/nyς-y1-ααqlE34

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

dz*dt=n+1nk-yς1/nς=1E35

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):

y(1ε,t)=[nk(n+1)dz*dt]nεE36
dz*dt=n+1nk-yς1/nς=1-εE37

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 initial conditions (30), (33) may be performed by well-developed methods, like the Runge-Kutta 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 Newtonian 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:

xn=πμq05h128E'C81/3,tn=π2μq02h16E'C52/3.

Herein, μ is the dynamic viscosity, C is the fluid-loss coefficient in the Carter’s leak-off term ql=2C/t-τ, τ is the time at which the fracture front reaches a point x. In accordance with (18), we have used the influx q0 per unit height, while Nordgren wrote his normalizing values in terms of the total influx q1 through the entire height h (q1=q0h).

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 dotted 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 performed by using ε-regularized BC (36) and SE (37) with ε=10-4. Specially designed numerical experiments have shown that the relative error does not exceed 103. No signs of instability are noted. Even though the calculations were performed in the MATLAB environment 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 suggested.

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].

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.

## 5. Extension to P3D models

In the previous section we have stated that the modified formulation provides prerequisites for efficient solving the Nordgren’s problem. The latter, being the basis for the P3D models, we may extend the efficient schemes to these models.

The detailed description of the P3D models is given by Mack and Warpinski [1]. Thus it is sufficient to list only those their features, which distinguish them from the PKN model (Fig. 1) and which are significant for the extension.

(i) In P3D models, the in-situ rock stresses are different in various layers while the fluid pressure is assumed constant in any vertical cross-section. This implies that, in contrast with the PKN model, the net-pressure is now not constant in a vertical cross-section. Thus the P3D models employ the pressure itself rather than the net-pressure. Alternatively, one may employ the difference of the fluid pressure with a fixed value of in-situ rock pressure, say the pressure in the pay-layer. Below to keep clear connection with the PKN model, we shall use this option and conditionally call the difference the ‘net-pressure’.

(ii) The linear dependence (13) between the net-pressure p and the average opening wav of a cross-section is changed to a non–linear dependence p=p(wav). To have clear connection with the initial PKN model, we re-write it as

p=krwavFp(wav)E38

where for sufficiently small wav, in particular near the fracture front, Fpwav=1. A specific form of the dependence (38) is found from solving plane-strain elasticity equation for a straight vertical crack under the fracture conditions of the form KIKIc at the upper and bottom tips. Although looking for this dependence may be involved, it is found in advance for a prescribed geometry of layered stratum, in-situ stresses in layers and critical SIFs. These preliminary calculations also give the positions zu(p) and zl(p) of the upper and bottom tips, respectively. Consequently, the height hfp of the fracture in a cross section is a known function of the pressure:

hfp=zu(p)-zl(p)E39

(iii) The dependence (12) between the particle velocity and the pressure, after averaging the velocity over the height of a cross-section, obtains a factor Fvwav:

vav=-kfwavn+1px1/nFvwavE40

where for sufficiently small wav, we have Fvwav=1. Since wav0 when approaching the fracture front, the SE becomes:

v*=-kfwavn+1px1/nx=x*E41

(iv) The continuity equation is integrated over the cross-sectional height. As a result, the total flux through a cross-section and the total leak-off replace, respectively, the flux and leak-off 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 divided by a fixed reference height hr, say the height of the pay-layer. Then denoting

w=Ahr,q=Avavhr=wvav,ql=Qlhr, E42

we have the modified continuity equation (15) in the unchanged form. The BC (18) at the inlet is also the same when denoting q0(t)=Q0(t)/hr, where Q0(t) is the prescribed total influx.

Take into account that by the definition of the average opening, we have A=wavhf, where hf is the fracture height in a considered cross-section. Thus, by the first of (42), w=Ahf/hr. Then, in view of (39) and (38), we may use the argument w instead of wav in the equations (38), (40) and(41) writing them, respectively, as

p=krwGp(w), v=-kfwn+1px1/nGvw, and v*=-kfwn+1px1/nx=x*,

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 Fpwav and Fvwav. They are such that Gpw=1 and Gvw=1 for sufficiently small w.

Finally, by introducing the variables y=w1/α and ς=x/x*, we arrive at the system similar to (27)-(33):

yt=yαx*vς+ςv*vx*yςy1ααql,E43
v=kx*1/nH(y)-yς1/nE44
yς,t0=y0ςE45
[yαv]ς=0=q0t,   E46
y(1,t)=0E47
v*t=dx*dt=kx*1/n-yς1/nς=1E48
x*t0=x*0E49

where

H(y)=[d(wGp(w))dw]1/nGv(w)E50

is a function to be found in advance through the functions Gpw  and Gvw  of the argument w=yα. By the properties of Gpw  and Gvw, we have H0 =1 which explains its absence in the SE (48).

The problem (43)-(49) differs from the problem (27)-(33) in the only detail: equation (44) for the velocity contains the function Hy defined by (50). The latter function, being smooth and tending to the unity at the front, the efficient numerical schemes, discussed in the previous section, are of use for the P3D models.

Comment. In some cases, it may be convenient to use the net-pressure rather than the opening. Since dp=dpdwdw, reformulation of the equations and computational schemes in terms of the modified pressure P=p1/α is obvious.

## 6. Conclusions

The discussion above demonstrates the analytical and computational advantages of using the modified formulation. The analytical advantages are evident from the obtained simple analytical solutions for the PKN and KGD models, which otherwise require involved calculations. The computational advantages include: (i) the possibility to use the well-established theory of propagating surfaces, (ii) avoiding deterioration of numerical solution caused by ill-posedness of the problem when neglecting the lag and fixing the fracture contour at a time step, (iii) avoiding singularities on the fluid front, (iv) the possibility to use highly efficient numerical schemes for the PKN and P3D models. These beneficial features are of significance when developing simulators able to efficiently solve truly 3D and pseudo-three-dimensional problems in real time. The work on simulators of both types is in progress.

## Acknowledgments

The authors appreciate support of the EU Marie Curie IAPP program (Grant # 251475). The first author (AL) also thankfully acknowledges support of the Russian Fund of Fundamental Investigations as concerns with the results on non-Newtonian fluids (Grant # 12-05-00140). Both authors are grateful to Piotr Kusmierczyk and Michal Wrobel for the help in performing many numerical experiments for the Nordgren’s problem with Carter’s leak-off.

## References

1. 1. Mack MG., Warpinski NR. Mechanics of hydraulic fracturing. In: Economides MJ., Nolte KG. (eds) Reservoir simulation, 3-rd edn. John Willey & Sons; 2000. p6.1-6.49.
2. 2. Linkov AM. Speed equation and its application for solving ill-posed problem of hydraulic fracturing. Doklady Physics 2011; 56 (8) 436-438.
3. 3. Linkov AM. Use of speed equation for numerical simulation of hydraulic fractures. Available at: http://arxiv.org/abs/1108.6146.
4. 4. Linkov AM. On efficient simulation of hydraulic fracturing in terms of particle velocity. Int. J. Engineering Sci. 2012; 52, 77-88.
5. 5. Linkov AM. Numerical modeling of hydraulic fractures: State of art and new results. In: Advanced Problems in Mechanics, APM 2012: proceedings of International XL Summer School-Conference, 1-7 July 2012, Institute for Problems of Mechanical Engineering, RAS; 2012. CD-ROM. p225-236.
6. 6. Mishuris G., Wrobel M., Linkov A. On modeling hydraulic fracture in proper variables: stiffness, accuracy, sensitivity. Int. J. Engineering Sci. 2012; 61, 10-23.
7. 7. Adachi J., Siebrits E., Peirce A., Desroches J. Computer simulation of hydraulic fractures. Int. J. Rock Mech. Mining Sci. 2007; 44, 739-757.
8. 8. Jamamoto K., Shimamoto T., Sukemura S. Multi fracture propagation model for a three-dimensional hydraulic fracture simulator. Int. J. Geomech. ASCE 2004; 1, 46-57.
9. 9. Nordgren RP. Propagation of a vertical hydraulic fracture. SPEJ Aug., 1972; 10(4), 306-314.
10. 10. Spence DA., Sharp PW. Self-similar solutions for elastohydrodynamic cavity flow. Proc. Roy Soc. London, Ser. A 1985; 400, 289-313.
11. 11. Kemp LF. Study of Nordgren’s equation of hydraulic fracturing. SPE 18959, Society of Petroleum Engineers 1989; 287-298.
12. 12. Adachi J., Detournay E. Self-similar solution of a plane-strain fracture driven by a power law fluid. Int. J. Numer. Anal. Methods Geomechanics 2002; 26, 579-604.
13. 13. Perkins TK., Kern LF. Widths of hydraulic fractures. J. Pet. Tech. September 1961; 937-949.
14. 14. Khristianovich SA., Zheltov VP. Formation of vertical fractures by means of highly viscous liquid. In: Proc. 4-th World Petroleum Congress, Rome, 1955. p. 579-586
15. 15. Geertsma J., de Klerk F. A rapid method of predicting width and extent of hydraulically induced fractures. J. Pet. Tech. December 1969; 1571-1581.
16. 16. Kresse O., Cohen C., Weng X. et al. Numerical modeling of hydraulic fracturing in naturally fractured formations. In: Proc. 5-th US Rock Mechanics Symposium, San Francisco, CA, June 26-29, 2011, American Rock Mechanics Association. Paper ARMA 11-363. 11p.
17. 17. Sethian JA. Level set methods and fast marching methods. Cambridge: Cambridge University Press; 1998.
18. 18. Hadamard J. Sur les problemes aux derivees partielles et leur signification physique. Princeton University Bulletin 1902; 49-52.
19. 19. Linkov AM. Analytical solution of PKN problem for non-Newtonian fluid. Russian Mining Science 2013; 1 (in press).
20. 20. Peirce A., Detournay E. An implicit level set method for modeling hydraulically driven fractures. Comput. Methods Appl. Mech. Engng. 2008; 197, 2858-2885.
21. 21. Kovalyshen Y., Detournay E. A reexamination of the classical PKN model of hydraulic fracture. Transp Porous Med. 2010; 81, 317–339.

## Notes

• 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.
• 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.
• 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.

Written By