An approximate solution for a penny-shaped hydraulic fracture that accounts for fracture toughness, fluid viscosity and leak-off.

This paper develops a closed-form approximate solution for a penny-shaped hydraulic fracture whose behaviour is determined by an interplay of three competing physical processes that are associated with fluid viscosity, fracture toughness and fluid leak-off. The primary assumption that permits one to construct the solution is that the fracture behaviour is mainly determined by the three-process multiscale tip asymptotics and the global fluid volume balance. First, the developed approximation is compared with the existing solutions for all limiting regimes of propagation. Then, a solution map, which indicates applicability regions of the limiting solutions, is constructed. It is also shown that the constructed approximation accurately captures the scaling that is associated with the transition from any one limiting solution to another. The developed approximation is tested against a reference numerical solution, showing that accuracy of the fracture width and radius predictions lie within a fraction of a per cent for a wide range of parameters. As a result, the constructed approximation provides a rapid solution for a penny-shaped hydraulic fracture, which can be used for quick fracture design calculations or as a reference solution to evaluate accuracy of various hydraulic fracture simulators.

Introduction

Hydraulic fractures are the fluid-filled cracks that propagate under the influence of a fluid pressure acting along the crack’s surface. The most common and well-known application of hydraulic fracturing is the stimulation of oil and gas wells to enhance production of hydrocarbons [1]. Other industrial applications include waste remediation process [2], waste disposal [3] and preconditioning in rock mining [4]. Hydraulic fractures also occur in nature in the process of magma ascent through the lithosphere owing to a buoyancy force [5-8,9] or as fluid-filled cracks in glacier beds [10].

Different hydraulic fracture geometries have been considered over the years. Research efforts shifted from the development of simple models, such as the Khristianovich–Zheltov–Geertsma–De Klerk (KGD) model for a plane-strain crack [11], towards more complex models that consider a planar hydraulic fracture [12-14], multiple hydraulic fractures [15-17] or fracture networks [18]. Detailed reviews of various hydraulic fracturing models can be found in [19-21]. In addition, there are techniques in which cracks are not modelled explicitly, such as the phase field [22,23], distinct element method [24,25] and peridynamics [26]. The primary advantage of such methodologies is the ability to tackle complex crack geometries easier than the conventional methods. At the same time, predictions of such approaches should be thoroughly tested against reference solutions to ensure that the new techniques are able to capture all the features of the conventional approaches.

Even for the simplest geometries, hydraulic fractures are known to obey a complex multiscale behaviour, see e.g. a thorough review paper [27]. This multiscale nature arises both in time, where multiple timescales determine fracture evolution, and space, where the solution may experience variations at different length scales in the tip region. As pointed out in [27], the time and length scales are related; that is, a particular timescale in the fracture evolution corresponds to dominance of one length scale in the tip region. Recognizing the significance and multiscale nature of the tip region, many studies were devoted specifically to quantify hydraulic fracture behaviour in the tip region [28-35]. On the other hand, time evolution and regimes of propagation were studied for a plane-strain KGD fracture in [36-40], and for a penny-shaped hydraulic fracture in [40-44]. A recent review paper [27] provides a more detailed summary of the findings and shows the complexity of the behaviour that occurs even for simple fracture geometries.

In view of the multiscale behaviour of hydraulic fractures, the primary goal of this paper is to quantify such a behaviour for the case of a penny-shaped hydraulic fracture, where the latter is driven by a Newtonian fluid and propagates in a permeable medium assuming no fluid lag. Most of the previous studies that addressed the problem of a penny-shaped fracture focused on the limiting regimes of propagation and asymptotic closed-form solutions (or accurate approximations) for the problem [40-43]. One exception is the study [44], in which a numerical solution for the full problem is obtained. The numerical solution is, however, relatively difficult to obtain owing to the temporal and spatial multiscale nature of the solution. In contrast, this study develops a closed-form approximate solution that provides results virtually instantly and accurately captures the complex behaviour of a radial hydraulic fracture at all length scales and timescales. In particular, the developed solution is able to describe all existing limiting solutions and all possible transitions between them, so that it covers all the trajectories in the parametric space for the problem under consideration. This development is made possible by using a closed-form approximation for the tip asymptotic solution obtained in [34], which is used to approximate fracture width profile. Once the fracture geometry is known, the global fluid volume balance is used to determine behaviour of the solution.

The importance of the obtained solution can be summarized as follows. It first shows that the tip region plays a crucial role in hydraulic fracture modelling. It also allows one to obtain a solution rapidly, which can be useful for quick estimations of fracture geometry for any values of fracture toughness, fluid viscosity and leak-off. Owing to a relatively simple implementation of the solution, it can be used as a reference solution to evaluate accuracy of other hydraulic fracturing simulators and, at the same time, as an initial condition to improve stability of the numerical schemes at early times. Finally, the developed approximation leads to construction of the solution map, which indicates validity regions of the limiting solutions and permits one to easily determine whether the solution for given problem parameters corresponds to one of the limiting cases.

This paper is organized as follows. Section 2 describes governing equations for a radial hydraulic fracture with leak-off. Section 3 outlines procedure for obtaining the approximate solution. Comparison of the developed approximation with the existing limiting solutions is presented in §4. Section 5 contains description of the solution structure, where applicability zones of the limiting solutions are indicated. Finally, §6 evaluates accuracy of the approximation by comparing its predictions to the reference numerical solution, which is followed by the summary.

Governing equations for a penny-shaped hydraulic fracture

This study considers propagation of an axisymmetric (‘radial’ or ‘penny-shaped’ are also used throughout the paper) hydraulic fracture in a permeable rock [27,44]. There are four primary material parameters that appear in the model, which can be introduced in the scaled form, for convenience, as where μ is the fluid viscosity, E is the Young modulus, ν is the Poisson’s ratio, KI is the mode I fracture toughness of the rock and CL is the Carter’s leak-off parameter.

Volume balance for an incompressible Newtonian fluid inside the crack can be written as where w(r,t) denotes the fracture width, q is the flux in the radial direction, the term that contains C′ accounts for leak-off via Carter’s model, t0(r) is the time instant at which the fracture front was located at point r, p is the fluid pressure and Q0 is the fluid injection rate (assumed to be constant in time).

The elasticity equation, which characterizes elastic equilibrium of the rock, relates the fluid pressure inside the crack to the fracture width as [41,44,45] where R is the fracture radius, and the kernel is The functions K(⋅) and E(⋅) denote the complete elliptic integrals of the first and the second kind, respectively.

Fracture propagation is modelled by a classical result of the linear elastic fracture mechanics (LEFM), in which the fracture opening in the tip region follows the square root solution [46] which implies that the stress intensity factor is equal to the fracture toughness for a propagating fracture. The propagation condition (2.5) should also be complemented by a zero flux condition at the fracture tip, i.e. q(R,t)=0.

For future reference, it is useful to consider the global fluid balance, which can be obtained by integrating (2.2) with respect to time and space as where the relations q(R,t)=0 and w(R,t)=0 were used to obtain the result.

Approximate solution for a penny-shaped hydraulic fracture

Outline of the methodology

To construct the approximate solution, it is assumed that the fracture evolution is primarily determined by the near-tip behaviour and the global fluid balance (2.6). In this situation, the approximation for the fracture width can be taken in the form where wa is the tip asymptotic solution that, in addition to distance to the tip R−r, depends on the material parameters (2.1) and time through R(t) and , whereas λ is the parameter that is determined later. Here, the tip asymptotic solution, wa, is the solution of (2.2), (2.3) and (2.5) in the tip region defined as (R−r)/R≪1, which, as shown in [13], is equivalent to the problem of a semi-infinite hydraulic fracture that propagates steadily under plane-strain elastic conditions [32,34]. As a result, because w(R−s,t)=wa(s) for s≪R, the approximation (3.1) precisely solves the governing equations (2.2), (2.3) and (2.5) in the tip region and approximates the solution inside the fracture away from the tip.

In order to construct the approximation, the closed-form approximation for wa, obtained in [34], is employed. As shown in [34], the three-process tip asymptotic solution, which captures the effects of fracture toughness, fluid viscosity and leak-off, has the maximum error of 0.14% for all regimes of propagation and all possible transitions between them. One of the properties of the tip solution, wa, is that , where s=R−r is the distance from a point inside the fracture to the tip and is a slowly varying function, which is also part of the solution (i.e. known). This fact allows to reduce (3.1) to Motivated by the results of the limiting asymptotic solutions, for which R(t)∝t and α is a number that is equal to either 1/4, 2/5, or 4/9 [44,27], it is further assumed that R(t)∝t, where α is a slowly varying function. With this approximation, the function t0(r) is determined from the relation r/R=(t0/t), which together with (3.2) reduces the global fluid balance (2.6) to where ρ=r/R is the scaled spatial coordinate. The integrals in equation (3.3) can be evaluated, in which case (3.3) reduces to where B(a,b) denotes the beta function, and where B(x;a,b) is the incomplete beta function.

Procedure for obtaining the approximate solution can be outlined as follows. First, material parameters (2.1), Q0 and the function wa(R) should be specified. Note that wa(R) also depends on the material parameters (2.1) and . Also note that the expression for comes from the tip asymptotic solution. Then, by specifying λ and taking an initial guess for α=4/9 (corresponds to zero toughness and zero leak-off solution), equation (3.4) can be solved for R(t) (e.g. by using Newton’s method). Once R(t) is obtained, the value of α is updated according to . After that, equation (3.4) is solved again. The iteration procedure is performed until the desired level of convergence is reached, which is typically achieved quickly owing to relatively modest dependence on α. Once R(t) is calculated, the width can be calculated using (3.2). Efficiency, which is defined as the ratio between the current fracture volume and the total amount of injected fluid, can be calculated as To calculate fluid pressure, one can substitute (3.2) into (2.3) and evaluate the integral, so that where the function can be evaluated numerically, and the kernel M(ρ,s) is given in (2.4). The accuracy of such pressure calculation is examined later in the paper.

It should be noted here that λ is taken as a parameter at this point, and its value can be varied to achieve a better approximation. At the same time, in order to capture the elliptic solution for a uniformly pressurized crack, one should set λ=1/2. As is shown later, the values of λ close to 1/2 provide the most accurate results.

Solution in scaled variables

To solve (3.4) numerically, it is useful to introduce scaling, in which the tip asymptotic solution wa is expressed. In particular, as shown in [14,17], the solution wa is given implicitly by the equation where g is a known function (see appendix A), and is used for velocity of the crack tip. By introducing the following scales and the dimensionless quantities the global fluid balance equation (3.4) can be rewritten as where the parameter is specified later, and (see appendix A). The first equation in (3.7), the first relation in (3.9) and (3.10) can be combined to yield the system of equations To obtain the solution, the system of equations (3.12) is solved numerically for and using Newton’s method. Note that is a parameter, and the initial value of α is taken as 4/9, which is then updated iteratively using (so that the system (3.12) is solved a few times until α converges).

Once the time histories and are obtained, equation (3.7) is used to find , and the relations in (3.9) are used to calculate and . The values of and λ (see §4.5) are also expressed in terms of and . The pressure and the fracture opening solutions can be obtained from (3.6) and (3.2), respectively, as whereas the efficiency can be expressed using (3.5), (3.10) and (3.11) as Finally, the solution is given in terms of the fracture radius , fracture width , fluid pressure and efficiency . Equations (3.8) and (3.9) can be used to convert the solution to the unscaled form if needed.

Comparison with vertex solutions

The problem of a penny-shaped fracture with fluid loss and no lag has four limiting regimes of propagation (or vertex solutions) [44,27]: the storage viscosity (M), storage toughness (K), leak-off viscosity () and leak-off toughness (). These limiting regimes are determined by dominance of one of the dissipative mechanisms associated with fluid viscosity and rock toughness, and by dominance of either the fluid storage in the fracture or fluid leak-off in the surrounding rock. The limiting regimes can be expressed using the two variables and 0≤η≤1 as Indeed, the case corresponds to the situation of no fracture toughness, while implies that the tip asymptotic solution follows the LEFM square root solution and hence, the effect of fluid viscosity is negligible. At the same time, the efficiency η=1 corresponds to the situation in which all fluid is stored in the fracture, whereas η=0 signifies the case in which all the injected fluid leaks into the formation. It is important to note that η=1 corresponds to and η=0 corresponds to , as can be seen from (3.11) and (3.14).

To evaluate accuracy of the approximate solution, and to specify the values for λ, the remainder of this section is devoted to a comparison between the limiting regimes (4.1) of the approximate solution and the existing reference solutions for these cases. In addition, setting some of the parameters, such as K′ and C′ (or and ), to zero makes the scaling (3.8) and (3.9) deficient. It is still possible, however, to use a combination of variables, such as , because it may not depend on and, consequently, on C′. This fact will be used for the remaining part of this section to write the limiting solutions in terms of scaling (3.8) and (3.9).

M vertex limit of the solution

Before proceeding with the M vertex limit, it is useful to rewrite equations (3.9) and (3.12) in the form Because the M vertex solution corresponds to the case of no toughness and no leak-off, it is possible to set in (4.2), in which case α=4/9. Note that and (see appendix A) are constants, where β=21/335/6.

By introducing scaling that is associated with the M vertex solution [13],41,44] equations (4.2) can be reduced to where α=4/9 (because L∝t4/9) is used to obtain the result. The solutions for the fracture width and pressure variations can be written using (3.13) as where the fact that is used. Comparison of the fracture width at the wellbore Ω(0) and fracture radius γ that are calculated using (4.4) and (4.5) with those obtained in [41] for the M vertex case implies that λ=0.487 provides the best approximation, for which the error for both Ω(0) and γ is below 0.5%.

The M vertex limit of the approximate solution can be summarized as which can be written in the dimensional form using (4.3) as where the function is defined in (3.6).

vertex limit of the solution

In order to obtain the limit, one should consider the limiting case of no fracture toughness, i.e. , and zero efficiency, i.e. η=0. In this situation, equation (3.14) implies that Because the zero efficiency case corresponds to large values of , one can use the asymptotic behaviour of g at large , namely (see appendix A), where . As a result, equations (3.9) and (3.12) can be rewritten in the form The combination of (4.8) and (4.9) finally leads to in which case it is clear that , and hence the beta function can be evaluated as .

By introducing scaling that is associated with the vertex solution [44,13] relations (4.10) can be reduced to and solutions for the fracture width and pressure variations can be written using (3.13) as where is used, which corresponds to the case and . Comparison of the fracture width at the wellbore and the fracture radius in (4.13) with the reference solution for the vertex [44,13] indicates that provides the most accurate approximation, in which the error of is below 0.1%, whereas the value of is captured precisely.

The vertex limit of the approximate solution can be summarized from (4.12) and (4.13) as and can be written in the dimensional form using (4.11) as where the function is defined in (3.6).

K vertex limit of the solution

The K vertex solution corresponds to the case when and η=1 (or ). In this situation, and equations (3.9) and (3.12) can be rewritten in the form which implies that α=2/5.

By introducing scaling that is associated with the K vertex solution [13,41,44] relations (4.16) reduce to in which case the width and pressure spatial variations can be obtained from (3.13) as where was used. Because the fracture width for the K vertex solution is an ellipse [41], then λ=1/2 leads to the exact solution.

Given that and , the solution can be summarized as which coincides with the K vertex solution in [41,44]. The above results can be written in the unscaled form as Note that because the fracture width for the K vertex limit is an ellipse and that the approximation (3.1) captures elliptic width profile precisely, the obtained limited solution is the exact solution for the K vertex limit and should not be tested against a reference solution.

The vertex limit of the solution corresponds to the situation for which and η=0. As a result, the efficiency relation (3.14) reduces to and equations (3.9) and (3.12) can be rewritten as As for the solution, , and the beta function can be evaluated as .

By introducing scaling that is associated with the vertex solution [13,42,44] relations (4.23) can be reduced to and the width and pressure spatial variations can be obtained using (3.13) as where, similarly to the K vertex limit, was used. As for the K vertex case, the fracture width for the vertex solution is an ellipse [42,44], and leads to the exact solution.

Because the results in (4.23) and (4.26) can be summarized as which coincides with the vertex solution in [42,44]. The unscaled form of (4.27) is Similar to the K vertex solution, the limit of the approximate solution coincides with the limit of the exact solution, because the approximation (3.1) captures elliptic width profile precisely.

Parameter λ interpolation

This section aims to describe an interpolation scheme to determine values of the parameter λ, knowing its values for the limiting cases. Comparison between the developed approximate solution and the vertex solutions indicates that λ=0.487, , λ=0.5 and . Owing to a small variation of λ, one may consider the following interpolation procedure: where is the approximation for the efficiency (3.14) for which λ=0.5 is used. The interpolation (4.29) captures the correct values of λ at all four vertices, and provides an approximation elsewhere. It is important to note that the choice of in lieu of e.g. is due to the fact that the function g depends on when (transition between and regimes); see appendix A. Note that the transition from M to K (η=0) is less important for the interpolation purposes, because λ≈λ and any interpolation would work. The transitions from M to and from K to are determined by the efficiency, hence the efficiency is used as a parameter for the interpolation in (4.29). Because even a constant approximation λ≡0.5 is able to provide relatively accurate approximation owing to a small variation of λ, the accuracy of the interpolation (4.29) is expected to be sufficient. This statement will be verified later in §6, where the approximate solution is compared with the reference numerical solution.

Structure of the solution

Solution in the original scaling is given by the quantities , , and η, which depend on the dimensionless time and the parameter . At this point, it is useful to rewrite all the variables in terms of the scaling that will be used in the reference numerical solution in §6. In particular, the scaling of the numerical solution is determined by the following quantities [44] whereas the leak-off parameter ϕ is defined as In this setting, the relationship between the new scaling and the original scaling (3.8) and (3.9) can be summarized as where the solution in new scaling is given by the fracture radius γ(τ), fracture width Ω(ρ,τ), fluid pressure Π(ρ,τ) and efficiency η(τ).

To visualize the structure of the solution, figure 1 plots the approximate width solution evaluated at the wellbore Ω(0,τ) versus τ for a range of values ϕ. Blue, red, green and magenta regions indicate, respectively, validity zones of the M vertex solution (4.6)–(4.7), K vertex solution (4.20)–(4.21), vertex solution (4.14)–(4.15) and vertex solution (4.27)–(4.28). These validity zones are defined as where i corresponds to either M, K, or vertex solutions. Note that all quantities in formula (5.4) are written in terms of the scaling (5.3).

Figure 1.

Variation of the approximate fracture width solution Ω(0,τ) versus dimensionless time τ and leak-off parameter ϕ. Blue, red, green and magenta regions indicate, respectively, validity zones of the M vertex solution (4.6)–(4.7), K vertex solution (4.20)–(4.21), vertex solution (4.14)–(4.15) and vertex solution (4.27)–(4.28), according to (5.4). Black lines indicate boundaries of applicability of the vertex solutions.

Figure 1 is a map of the solutions and allows one to determine the regime of propagation of a radial fracture for a given problem parameters. For instance, for the parameters E′=10 GPa, K′=1 MPa⋅m1/2, Q0=10−2 m3 s−1, t=103 s, μ′=1 Pa⋅s and C′=10−6 m s−1/2, the dimensionless time and leak-off parameters, respectively, are τ=10−5 and ϕ=1, in which case, according to the figure 1, the M vertex solution still can be used to estimate the fracture parameters. At the same time, if the viscosity is changed to μ′=0.01 Pa⋅s (other parameters are kept the same), then τ=1 and ϕ=10−6, so that none of the vertex solutions can be used. Instead, the developed approximation can be used. In addition, figure 1 can be used to indicate possible time evolution paths. The fracture starts at the M vertex and ultimately reaches the vertex. Along its way, the fracture may pass through the K vertex, through the vertex, or be in the transition region between the two, as was described in [44,27]. Black lines in figure 1 show the boundaries of applicability of the vertex solutions and are determined next. Note that all the timescales that are associated with the transitions from one vertex solution to another (obtained in the remainder of this section using limits of the approximate solution) are consistent with the results in [44], which were obtained using a different method. This indicates that the developed approximate solution precisely captures the scaling that is associated with the transitions from one vertex solution to another.

Transition along the MK edge. To quantify the boundaries between the M and K vertices, it is necessary to consider the MK edge transition. This MK edge solution corresponds to the situation of no leak-off, or η=1. As a result, and equations (3.12) and (3.14) reduce to Solution of equation (5.5) determines the MK edge solution. It is important to note that the solution depends on the parameter , which can be rewritten using (5.2) and (5.3) as τ. As a result, the dimensionless time that quantifies the MK transition is The associated boundaries of the MK transition, shown in figure 1, are calculated numerically as where corresponds to the beginning of the MK transition (near the M vertex), whereas corresponds to the end of the transition (near the K vertex).

Transition along the edge. To quantify the boundaries between the M and vertices, it is necessary to consider the edge transition. This edge solution corresponds to the situation in which . As a result, equations (3.12) reduce to which determines the edge solution. The solution depends on the parameter , which can be rewritten using (5.2) and (5.3) as τϕ9/14. In this case, the dimensionless time that quantifies the transition is The associated boundaries of the transition, shown in figure 1, are calculated numerically as where corresponds to the beginning of the transition (near the M vertex), whereas corresponds to the end of the transition (near the vertex).

Transition along the edge. To quantify the boundaries between the K and vertices, it is necessary to consider the edge transition. This edge solution corresponds to the situation in which . As a result, equations (3.12) reduce to which determines the edge solution. The solution depends on the parameter , which can be rewritten using (5.2) and (5.3) as τϕ5/6. In this case, the dimensionless time that quantifies the transition is The associated boundaries of the transition, shown in figure 1, are calculated numerically as where corresponds to the beginning of the transition (near the K vertex), whereas corresponds to the end of the transition (near the vertex).

Transition along the edge. To quantify the boundaries between the and vertices, it is necessary to consider the edge transition. This edge solution corresponds to the situation in which η=0 or . Given that for large (where ), equations (3.12) and (3.14) can be reduced to which determines the edge solution. The solution depends on the parameter , which can be rewritten using (5.2) and (5.3) as τϕ−1/2. In this case, the dimensionless time that quantifies the transition is The associated boundaries of the transition, shown in figure 1, are calculated numerically as where corresponds to the beginning of the transition (near the vertex), whereas corresponds to the end of the transition (near the vertex).

Comparison with numerical solution

To check the accuracy of the developed approximate solution, it is necessary to compare predictions of the approximation to a reference solution. In this study, the reference solution is calculated numerically by solving the original equations governing the problem of a radial hydraulic fracture (2.2) and (2.3). One of the key features of the scheme is the ability to capture multiscale tip behaviour by using the tip asymptotic solution as a propagation condition. Brief description of the numerical scheme is given in appendix B.

Because §4 examined accuracy of the approximation at the limiting cases of vertex solutions, the aim of this section is to check the accuracy in the transition regions. For this reason, the choice of the parameters for comparison is ϕ={10−18,10−15,10−12,10−9,10−6,10−3,1, 103,106} and 10−7<τ<107, which cover the transition region according to figure 1. It is important to note that the numerical solution is more difficult to obtain near the K and vertices, because the pressure is nearly constant, in which case the fluid flow within the fracture becomes irrelevant. At the same time, because the approximate width solution is exact at the K and vertices (and also transition), the approximation is expected to be very accurate in these regions. In addition, the numerical scheme becomes less computationally efficient (i.e. requires finer mesh for convergence) near the vertex solution. This occurs, because the efficiency becomes very small, in which case most of the injected fluid is leaks into the formation. The accuracy of the approximation at the vertex, however, was evaluated earlier and there is no need to check it again with the numerical solution.

Figure 2 shows comparison between the developed approximate solution (dashed grey lines) and the numerical solution (solid black lines) in terms of time histories of the fracture radius (panel (a)), efficiency (panel (b)), width at the wellbore (panel (c)) and pressure at ρ=0.5 (panel (d)), all for different values of ϕ (arrows indicate the direction in which the curves shift as ϕ increases). The fracture width and pressure are evaluated as where the pressure is not evaluated at the origin owing to a singularity located at that point. All the quantities are normalized by the M vertex solution (4.6), which can be written in terms of the scaling (5.1) and (5.3) as The purpose of such normalization is to achieve better visual separation between the curves in figure 2. The vertex solutions are shown by the blue (M vertex), red (K vertex), green ( vertex) and magenta ( vertex) lines. All the solutions originate near the M vertex (i.e. the blue lines), and then transition to either K, , or intermediate solution. Note that, eventually, all solutions (except the ϕ=0 case) reach the limiting case, as can be clearly seen from figure 1.

Figure 2.

Comparison between the numerical solution (solid black lines) and the approximation (dashed grey lines) in terms of time histories of: (a) fracture radius, (b) efficiency, (c) width at the wellbore (ρ=0) and (d) fluid pressure at ρ=0.5. The fracture radius, width and pressure are normalized by the M vertex solution (6.1). Results are shown for ϕ={10−18,10−15,10−12,10−9,10−6,10−3,1,103,106}, where arrows indicate the directions in which the curves shift as the leak-off parameter ϕ increases. Dashed blue, red, green and magenta lines indicate the M, K, and limiting solutions, respectively.

As can be seen from figure 2, there is no visual difference between the numerical solution and the approximation, which indicates that the developed approximate solution is accurate. Despite this fact, it is still important to quantify the level of accuracy of the approximation. To help doing this, the error is defined as where 〈⋅〉 indicates average over all time points, subscripts ‘apr’ and ‘num’ correspond to the approximation and numerical solution respectively, whereas A is either radius γ(τ), efficiency η(τ), wellbore width Ω0=Ω(0,τ) or pressure Π0=Π(0.5,τ).

Figure 3a shows the error that is calculated using (6.2) for different values of the parameter ϕ for the finest mesh (for numerical solution) considered. The magnitude of the error for γ, Ω0 and Π0 does not vary notably versus ϕ, which indicates that the interpolation procedure for λ (4.29) is sufficiently accurate. The accuracy of all quantities, but the pressure, is under a per cent. The pressure is less accurate because it is calculated using the elasticity integral (2.3), which is very sensitive to the fracture width profile. To better understand the effect of mesh size that is used to calculate numerical solution, figure 3b shows the maximum error versus the number of spatial discretization points of the numerical solution N. Here the maximum is calculated for different values of ϕ, which can be obtained from figure 3a and its analogues for different meshes. Results in figure 3b demonstrate that error for the fluid pressure Π0 and fracture length γ reached the plateau values, but the efficiency η and the wellbore width Ω0 did not reach their respective plateau values for N=200. Despite this fact, it is still clear that the error of the developed approximation is under a per cent (for the considered error measure), which is sufficient for most practical cases.

Figure 3.

(a) Time-average discrepancy between the developed approximation and the numerical solution calculated using (6.2) for the fracture radius γ(τ), efficiency η(τ), wellbore width Ω0=Ω(0,τ) and pressure Π0=Π(0.5,τ). (b) Variation of the maximum error (maximum is calculated over ϕ) versus number of mesh points N of the numerical solution for the fracture length γ, efficiency η, wellbore width Ω0 and pressure Π0.

In order to check accuracy of the pressure and width spatial variations, figure 4 plots pressure and width predictions stemming from the approximate solution (dashed grey lines) and those calculated numerically (solid black lines) for ϕ=10−15 (panel (a)), ϕ=10−6 (panel (b)) and ϕ=103 (panel (c)). Both the pressure and the width are normalized by the M vertex solution (6.1) for convenience. The vertex solutions are shown by the blue (M vertex), red (K vertex), green ( vertex) and magenta ( vertex) lines for the first and the last value of τ. For the ϕ=10−15 case (panel (a)), the solution transitions from the M vertex and almost reaches the K vertex. For the ϕ=10−6 case (panel (b)), the solution transitions from the M vertex and approaches the vertex. For the ϕ=103 case (panel (c)), the solution transitions from the M vertex to the vertex. These observations agree with the results shown in figure 1. Fracture width variations are consistently accurate, whereas the fluid pressure estimations deviate notably in the tip region. This is the weakest point of the approximate solution, and, therefore, it should not be used to evaluate pressure for . There might be a possibility to overcome this error by using (2.2) to estimate the pressure field. This, however, is a minor correction to the developed procedure and therefore is beyond the scope of this study.

Figure 4.

Fracture width and fluid pressure spatial variations calculated numerically (solid black lines) and obtained from the approximate solution (dashed grey lines) for: (a) ϕ=10−15 and τ={1,10,102,103,104,105}, (b) ϕ=10−6 and τ= {10,102,103,104,105,106}, (c) ϕ=103 and τ={10−4,10−3,10−2,10−1,1,10}. The fracture width and pressure are normalized by the M vertex solution (6.1). Dashed blue, red, green and magenta lines indicate the M, K, and limiting solutions, respectively.

Summary

This paper presents a closed-form approximate solution for a penny-shaped hydraulic fracture, whose behaviour is governed by a simultaneous interplay of fluid viscosity, fracture toughness and fluid leak-off into the formation. This development is made possible by combining the approximation for the fracture width profile, which uses the multiscale tip asymptotic solution in combination with the global fluid volume balance. First, the limiting regimes of the approximation are analyzed and the results are compared with the existing solutions. Then, a solution map is presented. This solution map indicates the areas of applicability of the limiting solutions and allows one to easily determine whether a problem with given parameters corresponds to one of the limiting cases or not. In addition, the developed solution captures the scaling that is associated with the transition from one limiting solution to another. In order to estimate accuracy of the approximate solution, a reference numerical solution is constructed. The latter numerical solution uses a moving mesh together with the implicit time stepping and a multiscale tip asymptotic solution to locate the moving fracture front. Results indicate that predictions of the fracture width and radius by the developed approximation are accurate to within a fraction of a per cent. The fluid pressure results are less accurate and lie within 2% error if measured half-way to the fracture tip. The importance of the developed approximate solution is twofold. First, it can be used to estimate fracture radius very quickly, in which case it can be used for a rapid fracture design calculations. Second, it can be used as a reference solution to test accuracy of more complex hydraulic fracture simulators (such as non-axisymmetric or non-planar), and can also be used as an initial condition for the aforementioned numerical codes.