Mathematical Model of a Hyperbolic Hydraulic Fracture with Tortuosity

M. R. R. Kgatle-Maseko, D. P. Mason

School of Computer Science and Applied Mathematics, DST-NRF Centre of Excellence in Mathematical and Statistical Sciences, University of the Witwatersrand, Johannesburg, South Africa.

DOI: 10.4236/jamp.2021.95078                

Abstract

The aim of the research is to study the propagation of a hydraulic fracture with tortuosity due to contact areas between touching asperities on opposite crack walls. The tortuous fracture is replaced by a model symmetric partially open fracture with a hyperbolic crack law and a modified Reynolds flow law. The normal stress at the crack walls is assumed to be proportional to the half-width of the model fracture. The Lie point symmetry of the nonlinear diffusion equation for the fracture half-width is derived and the general form of the group invariant solution is obtained. It was found that the fluid flux at the fracture entry cannot be prescribed arbitrarily, because it is determined by the group invariant solution and that the exponent n in the modified Reynolds flow power law must lie in the range 2 < n < 5. The boundary value problem is solved numerically using a backward shooting method from the fracture tip, offset by 0 < δ ≪ 1 to avoid singularities, to the fracture entry. The numerical results showed that the tortuosity and the pressure due to the contact regions both have the effect of increasing the fracture length. The spatial gradient of the half-width was found to be singular at the fracture tip for 3 < n < 5, to be finite for the Reynolds flow law n = 3 and to be zero for 2 < n < 3. The thin fluid film approximation breaks down at the fracture tip for 3 < n < 5 while it remains valid for increasingly tortuous fractures with 2 < n < 3. The effect of the touching asperities is to decrease the width averaged fluid velocity. An approximate analytical solution for the half-width, which was found to agree well with the numerical solution, is derived by making the approximation that the width averaged fluid velocity increases linearly with distance along the fracture.

Keywords

Hyperbolic Crack Law, Partially Open Fracture with Tortuosity, Modified Reynolds Flow Law, Backward Shooting Method, Approximate Analytical Solution

1. Introduction

In hydraulic fracturing, fluid is pumped at high pressure into a crack in a rock mass in order to open the crack. Hydraulic fracturing has many applications. In mining, high pressure water is used to break rocks instead of explosives which leave small rock fragments in the air that could damage the lungs of miners [1]. In geothermal reservoirs, in their natural state, the cracks allow only a small flow of water. In order to increase the flow, high pressure water is pumped into the crack network at one borehole and extracted from another borehole [2]. Hydraulic fracturing is also used to enhance the extraction of oil and gas in large underground shale deposits [3].

Spence and Sharp [4] were the first to show that the equations of hydraulic fracture admit a similarity solution. These authors considered the enlargement of a lens shaped crack and a two-dimensional crack by a viscous fluid modelled by lubrication theory. The fluid pressure and the crack shape were connected by a singular integral equation from linear elasticity. The theory was applied to magma-driven propagation of cracks in geophysics [5] [6] [7].

The Cautchy principal value integral in the singular integral equation relating pressure to the crack shape is difficult to analyze. An important simplification was the Perkins-Kern-Nordgren (PKN) approximation in which the normal stress at the fluid-rock interface is proportional to the half-width of the fracture [8] [9]. The PKN approximation puts the differential equations of hydraulic fracture in a form which can be analyzed using the theory of Lie point symmentries and conservation laws.

The scientific literature on hydraulic fracturing is now very large. We will only comment briefly on the application of Lie point symmetries to hydraulic fracturing, because that approach will be used in this paper. Lie point symmetry analysis was first applied to hydraulic fracturing by Fitt, Mason and Moss who considered the propagation of a fracture in impermeable rock [1]. Fareo and Mason analyzed the propagation of a hydraulic fracture in permeable rock with fluid leak-off into the rock mass and investigated the effect of leak-off on the rate of propagation of a fracture [10]. Fareo and Mason also considered a fracture driven by a power law fluid in an impermeable rock mass [11]. Anthonyrajah, Mason and Fareo [12] compared laminar and turbulent fluid driven fractures using the wall shear stress model of Emerman, Turcotte and Spence [6] for the fluid and the PKN approximation instead of the Cautchy principal value integral model to relate the fluid pressure to the crack shape.

The effects of tortuosity on hydraulic fracturing due to asperities or surface roughness at the fluid rock interface and contact regions caused by touching asperities will be investigated in this paper. The hydraulic fracture with tortuosity will be replaced by a symmetric two-dimensional hydraulic fracture without asperities but with a modified Reynolds flow law, which accounts for the effect that the presence of asperities at the fluid rock interface has on the fluid flow, and with a modified crack law, which models the effect of the contact regions on the stress at the fluid-rock interface [2].

In this paper, we model the contact regions by the hyperbolic crack law which was introduced by Goodman [13]. This is motivated by the discussion by Fitt et al. [2] that the hyperbolic crack law is generally considered to be a more realistic model to describe the presence of contact regions (deformations formed by touching asperities) in a fracture than the linear crack law. Although the concept of tortuosity was investigated in [14] [15], the contact regions in these papers were modelled by the linear crack law which was proposed by Pine and Cundall [16] and further discussed by Fitt et al. [2].

In this paper, we aim to solve the problem of the evolution of a hydraulic fracture with tortuosity described by the hyperbolic crack law. The fluid flow in the fracture is described by a modified Reynolds flow law and the fluid pressure and crack shape are related by the PKN approximation. The aim of this work is to derive numerical and approximate analytical solutions for the evolution of the half-width and length of the model symmetric fracture which replaces the fracture with tortuosity. The problem is formulated mathematically and in dimensionless form in Section 2.6.

In Section 2, we present the derivation of the model describing fluid flow in a tortuous hydraulic fracture with contact regions modelled by the hyperbolic crack law. Section 3 outlines the derivation of the group invariant solution of the problem. In Section 4 we consider possible operating conditions at the fracture entry and investigate the existence of analytical solutions. Section 5 describes a numerical solution for the problem. In Section 6, the investigation of the width averaged fluid velocity leads to the derivation of an approximate analytical solution for the problem. In this Section, a comparison of the numerical and the approximate analytical solutions is made. Finally in Section 7 we summarize important findings and conclude the paper.

2. Model formulation

We will first review the fluid flow results that were derived in [14]. We will then close the fluid flow model by using both the hyperbolic crack law, which describes the presence of contact regions in the fracture, and the Perkins-Kern [8] and Nordgren [9] approximation, which relates the normal stress at the fracture walls with the half-width of the fracture. Lastly, we will express the governing equations in dimensionless form.

A tortuous hydraulic fracture with asperities at the fluid-rock interface is modelled. The tortuous fracture is replaced by a symmetric two-dimensional hydraulic fracture without asperities, as shown in Figure 1, but with a modified Reynolds flow law and a modified crack law. The x-axis is along the length of the model symmetric fracture, the z-axis is along the width and the y-axis is along the breadth. All quantities are independent of y. The half-width of the fracture is

Figure 1. Hydraulic fracture without asperities on the crack walls or contact regions.

$h\left(t\mathrm{,}x\right)$. Since the fracture is long and thin, the lubrication approximation is imposed on the viscous flow in the fracture which gives that the fluid pressure in independent of the variable z. It follows that the fluid velocity components and the fluid pressure respectively have the form

$v_x=v_x\left(t,x,z\right),v_y=0,v_z=v_z\left(t,x,z\right),p=p\left(t,x\right).$ (2.1)

The fluid is incompressible with constant density $\rho$ and dynamic viscosity $\mu$. The body force due to gravity is neglected.

2.1. Review of the Flow Model

From lubrication theory for a two-dimensional hydraulic fracture the volume flux of fluid across the fracture per unit breadth, $Q\left(t\mathrm{,}x\right)$, satisfies the Reynolds flow law [14]:

$Q\left(t,x\right)={\displaystyle {\int }_{-h\left(t,x\right)}^{h\left(t,x\right)}}{v}_x\left(t,x,z\right)\text{d}z=-\frac{2}{3\mu }{h}^3\frac{\partial p}{\partial x}$ (2.2)

and the width averaged fluid velocity is

${\bar{v}}_x\left(t,x\right)=\frac{Q\left(t,x\right)}{2h\left(t,x\right)}=-\frac{h^2}{3\mu }\frac{\partial p}{\partial x}.$ (2.3)

The half-width of the fracture $h\left(t,x\right)$ satisfies

$\frac{\partial h}{\partial t}+\frac{1}{2}\frac{\partial }{\partial x}Q\left(t\mathrm{,}x\right)=\mathrm{0,}$ (2.4)

which when expanded with the aid of (2.2) gives the nonlinear diffusion equation relating half-width of the fracture $h\left(t,x\right)$ and the fluid pressure $h\left(t,x\right)$ :

$\frac{\partial h}{\partial t}=\frac{1}{3\mu }\frac{\partial }{\partial x}\left(h^3\frac{\partial p}{\partial x}\right).$ (2.5)

For hydraulic fractures with high tortuosity the Reynolds flow law has been shown to give theoretical results which diverge from experimental results [17] [18] [19] [20]. It follows that the Reynolds flow law is not satisfactory in describing very tortuous fractures. Several authors have considered various flow models with the aim of accounting for high tortuosity in the fracture. We will consider a law that governs one of these models, the modified Reynolds flow law, in which $h^3$ in $Q\left(t\mathrm{,}x\right)$ given in (2.2) is replaced by $a_n{h}^n$ where n is a dimensionless constant and $a_n$ has dimensions $L^{3-n}$. Both n and $a_n$ are obtained through experiments [2]. Equations (2.2), (2.3) and (2.5) become:

$Q\left(t,x\right)=-\frac{2}{3\mu }{a}_n{h}^n\frac{\partial p}{\partial x},$ (2.6)

${\bar{v}}_x\left(t,x\right)=-\frac{a_n}{3\mu }{h}^{n-1}\frac{\partial p}{\partial x}\mathrm{,}$ (2.7)

and

$\frac{\partial h}{\partial t}=\frac{a_n}{3\mu }\frac{\partial }{\partial x}\left(h^n\frac{\partial p}{\partial x}\right)\mathrm{.}$ (2.8)

2.2. Hyperbolic Crack Law

In [14] the linear crack law was used to describe the contribution made by the contact areas to support the compressive normal stress, $-{\sigma }_{zz}\left(t,x\right)$, at the fracture walls. In this paper we will use the hyperbolic crack law which was first introduced by Goodman [13] and further considered by Bandis et al. [21], Murphy et al. [22] and Fitt et al. [2].

A tortuous fracture can either be open without contact regions or partially open with contact regions. In an open fracture, the fluid pressure $p\left(t,x\right)$ is sufficient to support the compressive normal stress at the fracture walls, $-{\sigma }_{zz}\left(t,x\right)$, and therefore

$p\left(t,x\right)=-{\sigma }_{zz}\left(t,x\right).$ (2.9)

For a partially open fracture, the fluid pressure is insufficient to support the compressive normal stress at the fracture walls and therefore

$p\left(t,x\right)<-{\sigma }_{zz}\left(t,x\right).$ (2.10)

From (2.9) and (2.10), it is clear that the effective stress, $p\left(t,x\right)+{\sigma }_{zz}\left(t,x\right)$, vanishes for an open fracture and is negative for a partially open fracture. While both an open and a partially open fracture were discussed in [14] and [15], in this paper we will discuss only a partially open fracture. We consider the hyperbolic crack law [13],

${\sigma }_{zz}\left(t,x\right)+p\left(t,x\right)=k\left(\frac{h_{\max }-h\left(t,x\right)}{h\left(t,x\right)-h_{\min }}\right),k<0,$ (2.11)

where k is a negative constant, $h_{\min }$ is the minimum half-width of the fracture and $h_{\max }$ is the maximum half-width of the fracture. We make the approximation that $h_{\min }=0$ which simplifies (2.11) to

${\sigma }_{zz}\left(t,x\right)+p\left(t,x\right)=k\left(\frac{h_{\max }}{h\left(t,x\right)}-1\right),k<0.$ (2.12)

Although in practice $h_{\min }$ can never be zero even when the fracture is under the largest possible compressive stresses, $h_{\min }=0$ is a reasonable approximation since $h_{\min }\ll h_{\max }$ [2] [23]. It is clear that Equation (2.12) satisfies the negative effective stress condition, $p\left(t,x\right)+{\sigma }_{zz}\left(t,x\right)<0$. Equation (2.12) can be rearranged to give

$p\left(t,x\right)+p_1\left(t,x\right)=-{\sigma }_{zz}\left(t,x\right),$ (2.13)

where

$p_1\left(t,x\right)=-k\left(\frac{h_{\max }}{h\left(t,x\right)}-1\right)>0,k<0,$ (2.14)

which clearly shows that the compressive normal stress at the fluid-rock interface, $-{\sigma }_{zz}\left(t,x\right)$, is supported by both the fluid pressure $p\left(t,x\right)$ and the pressure due to contact regions $p_1\left(t,x\right)$.

2.3. The PKN Approximation

Lastly, it is necessary to specify a relation between the compressive normal stress at the fluid-rock interface, $-{\sigma }_{zz}\left(t,x\right)$, and the half-width of the fracture, $h\left(t,x\right)$. We will use the Perkins-Kern-Nordgren (PKN) approximation, which is widely used in the oil and gas industry to relate the normal stress at the crack walls to the half-width of the fracture, [8] [9]:

$-{\sigma }_{zz}\left(t,x\right)=-{\sigma }_{zz}^{\infty }+\Lambda h\left(t,x\right),$ (2.15)

where $-{\sigma }_{zz}^{\left(\infty \right)}>0$ is the compressive stress at infinity (far-field compressive stress) within the rock and

$\Lambda =\frac{E}{\left(1-{\nu }^2\right)B}\mathrm{,}$ (2.16)

where E is the Young’s modulus and $\nu$ is the Poisson ratio of the rock encasing the fracture and B is the breadth of the fracture. The PKN approximation has the disadvantage that the stress intensity factor vanishes at the fracture tip, $x=L\left(t\right)$, because the half-width is zero at the tip but it can be applied to a single-sided fracture with non-zero initial length. This is unlike the relation,

${\sigma }_{zz}\left(t\mathrm{,}x\right)={\sigma }_{zz}^{\left(\infty \right)}+\frac{G}{2\pi \left(1-\nu \right)}{{\displaystyle \sout{\int }}}_{-\infty }^{\infty }\frac{\partial h}{\partial s}\left(t\mathrm{,}s\right)\frac{\text{d}s}{s-x}\mathrm{,}$ (2.17)

which describes a two sided fracture expanding from a point source, where the bar on the integral sign denotes the Cauchy principal value and G and $\nu$ are the shear modulus and Poisson ratio of the rock mass respectively. Equation (2.17) was first introduced by Spence and Sharp [4] and used by Fitt et al. [2] when considering a geothermal reservoir. While it is possible to derive a similarity solution for the integro-differential equation for the half-width of the fracture using the relation (2.17), the resulting boundary value problem is difficult to solve [4] [5] [6] [7]. By expanding (2.17), Adachi and Peirce have shown that the PKN approximation is a good approximation away from an ε-neighbourhood of the fracture tip where $\epsilon \ll 1$ [24].

2.4. Model Closure

From equations (2.13) and (2.15), the fluid pressure in the fracture is

$p\left(t,x\right)=-{\sigma }_{zz}^{\left(\infty \right)}+\Lambda h\left(t,x\right)+k\left(\frac{h_{\max }}{h\left(t,x\right)}-1\right)$ (2.18)

and therefore

$\frac{\partial p}{\partial x}\left(t,x\right)=\Lambda \left(1-\frac{k{h}_{\max }}{\Lambda h^2\left(t,x\right)}\right)\frac{\partial h}{\partial x}.$ (2.19)

Since $k<0$ the magnitude of the pressure gradient in the fluid along the fracture is increased by the contact regions. Equations (2.6) and (2.7) for the volume flux of fluid in the fracture and the width averaged fluid velocity respectively become

$Q\left(t,x\right)=-\frac{2a_n\Lambda }{3\mu }\left(h^n-\frac{k{h}_{\max }}{\Lambda }{h}^{n-2}\right)\frac{\partial h}{\partial x},$ (2.20)

${\bar{v}}_x\left(t,x\right)=-\frac{a_n\Lambda }{3\mu }\left(h^{n-1}-\frac{k{h}_{\max }}{\Lambda }{h}^{n-3}\right)\frac{\partial h}{\partial x},$ (2.21)

and (2.8) becomes the nonlinear diffusion equation for the half-width of the fracture, $h\left(t,x\right)$ :

$\frac{\partial h}{\partial t}=\frac{a_n\Lambda }{3\mu }\frac{\partial }{\partial x}\left[\left(h^n-\frac{k{h}_{\max }}{\Lambda }{h}^{n-2}\right)\frac{\partial h}{\partial x}\right].$ (2.22)

When $k=0$, the effective stress vanishes and the condition (2.9) for an open fracture is satisfied. Equations (2.20) to (2.22) reduce to the governing equations for an open fracture which was discussed in [14] and [15]. In this paper we consider only partially open fractures ( $k<0$ ) for which the compressive normal stress at the fracture walls is supported by both the fluid pressure $p\left(t,x\right)$ and the pressure $p_1\left(t,x\right)$ due to contact regions given by (2.14).

Equation (2.22) is a nonlinear diffusion equation of the form

$\frac{\partial h}{\partial t}=a_1\frac{\partial }{\partial x}\left[D\left(h\right)\frac{\partial h}{\partial x}\right],$ (2.23)

where the diffusion coefficient $D\left(h\right)$ satisfies

$D\left(h\right)=h^n+a_2{h}^{n-2}\mathrm{,}$ (2.24)

and $a_1$ and $a_2$ are positive constants. Many problems in heat propagation, diffusion of fluids in porous media and ground water flow are of the form (2.23) but with diffusion coefficients different from (2.24) [25]. To the best of our knowledge there is no application to hydraulic fracturing with a diffusion coefficient of the form (2.24).

2.5. Boundary and Initial Conditions

Consider first the boundary conditions at the fracture tip, $x=L\left(t\right)$, where $L\left(t\right)$ is the length of the fracture. The half-width of the fracture vanishes at the tip,

$x=L\left(t\right):h\left(t,L\left(t\right)\right)=0,t\ge 0.$ (2.25)

The second boundary condition at the fracture tip is that the volume flux of fluid vanishes:

$x=L\left(t\right):Q\left(t,L\left(t\right)\right)=0,t\ge 0,$ (2.26)

and therefore using (2.20),

$x=L\left(t\right):\left(h^n\left(t,L\left(t\right)\right)-\frac{k{h}_{\max }}{\Lambda }{h}^{n-2}\left(t,L\left(t\right)\right)\right)\frac{\partial h}{\partial x}\left(t,L\left(t\right)\right)=0.$ (2.27)

The boundary conditions (2.25) and (2.27) are moving boundary conditions [26].

Consider next the condition at the fracture entry, $x=0$. We consider a partially open fracture. Let

$x=0,t=0:\frac{h\left(0,0\right)}{h_{\max }}=\beta ,0<\beta <1.$ (2.28)

The parameter $\beta$ is prescribed and describes the extent to which the initial fracture is partially open. The boundary condition at the fracture entry is

$-\frac{2a_n\Lambda }{3\mu }\left(h^n\left(t,0\right)-\frac{k{h}_{\max }}{\Lambda }{h}^{n-2}\left(t,0\right)\right)\frac{\partial h}{\partial x}\left(t,0\right)=S\left(t\right),$ (2.29)

where $S\left(t\right)=Q\left(t,0\right)$ is the flux of fluid per unit breadth into the fracture at the entry $x=0$. It is prescribed as far as the invariant solution will allow.

2.6. Dimensionless Governing Equations

In order to make all equations and boundary and initial conditions dimensionless we transform to dimensionless variables. Consider first the fluid pressure. Equation (2.18) can be written as

$P\left(t,x\right)=\Lambda h\left(t,x\right)+k\left(\frac{h_{\max }}{h\left(t,x\right)}-1\right)$ (2.30)

where

$P\left(t,x\right)=p\left(t,x\right)+{\sigma }_{zz}^{\left(\infty \right)}\mathrm{.}$ (2.31)

We will work with $P\left(t,x\right)$ instead of $p\left(t,x\right)$. The pressure $P\left(t,x\right)$ is the difference between the pressure of the fluid in the fracture, $p\left(t,x\right)$, and the background pressure, $-{\sigma }_{zz}^{\left(\infty \right)}$, in the rock. Since ${\sigma }_{zz}^{\left(\infty \right)}$ is a constant, equations (2.3) to (2.8) on which the theory is based, remain unchanged when expressed in terms of $P\left(t,x\right)$ instead of $p\left(t,x\right)$.

We choose characteristic quantities which are independent of n and $a_n$ so that the results for a range of values of n can be compared. We choose $h_{\max }$ as the characteristic half-width. From (2.30) the characteristic value of $P\left(t,x\right)$ is $\Lambda h_{\max }$. The characteristic value of the x-coordinate is the initial length of the fracture $L_o=L\left(0\right)$ which is also the characteristic length of the fracture. Equation (2.5) was derived using lubrication theory and gives for the characteristic velocity along the fracture:

$U=\frac{\Lambda h_{\max }^3}{\mu L_o}\left(1-\frac{{\sigma }_{zz}^{\left(\infty \right)}}{\Lambda h_{\max }}\right).$ (2.32)

We transform to the following dimensionless variables:

$\begin{array}{l}{x}^{*}=\frac{x}{L_o},t^{*}=\frac{Ut}{L_o},h^{*}=\frac{h}{h_{\max }},L^{*}=\frac{L}{L_o},\\ {\bar{v}}_x^{*}=\frac{{\bar{v}}_x}{U},Q^{*}=\frac{Q}{U{h}_{\max }},S^{*}=\frac{S}{U{h}_{\max }},\\ V^{*}=\frac{V}{h_{\max }{L}_o},P^{*}=\frac{P}{\Lambda h_{\max }},\end{array}$ (2.33)

where

$V\left(t\right)=2{\displaystyle {\int }_0^{L\left(t\right)}}h\left(t\mathrm{,}x\right)\text{d}x$ (2.34)

is the total volume of the fracture per unit breadth and

$L^{*}\left(0\right)=1.$ (2.35)

We define

$K_n=\frac{a_n{h}_{\max }^{n-3}}{3}{\left(1-\frac{{\sigma }_{zz}^{\left(\infty \right)}}{\Lambda h_{\max }}\right)}^{-1},\varphi =-\frac{k}{\Lambda h_{\max }}>0.$ (2.36)

We suppress the stars on the dimensionless variables and parameters to keep the notation simple. Expressed in dimensionless form the half-width of the fracture, $h\left(t,x\right)$, satisfies the nonlinear diffusion equation

$\frac{\partial h}{\partial t}=K_n\frac{\partial }{\partial x}\left[\left(h^n+\varphi h^{n-2}\right)\frac{\partial h}{\partial x}\right],$ (2.37)

subject to the boundary conditions

$x=L\left(t\right):h\left(t,L\left(t\right)\right)=0,$ (2.38)

$x=L\left(t\right):\left(h^n\left(t,L\left(t\right)\right)+\varphi h^{n-2}\left(t,L\left(t\right)\right)\right)\frac{\partial h}{\partial x}\left(t,L\left(t\right)\right)=0,$ (2.39)

$x=0:-2K_n\left(h^n\left(t,0\right)+\varphi h^{n-2}\left(t,0\right)\right)\frac{\partial h}{\partial x}\left(t,0\right)=S\left(t\right).$ (2.40)

and the initial/boundary condition

$x=0,t=0:h\left(0,0\right)=\beta ,0<\beta <1.$ (2.41)

In dimensionless variables, the volume flux of fluid in the fracture, the width averaged fluid velocity, the fluid pressure and the volume of the fracture are

$Q\left(t,x\right)=-2K_n\left(h^n\left(t,x\right)+\varphi h^{n-2}\left(t,x\right)\right)\frac{\partial h}{\partial x}\left(t,x\right),$ (2.42)

${\bar{v}}_x\left(t,x\right)=-K_n\left(h^{n-1}\left(t,x\right)+\varphi h^{n-3}\left(t,x\right)\right)\frac{\partial h}{\partial x}\left(t,x\right),$ (2.43)

$P\left(t,x\right)=h\left(t,x\right)-\varphi \left(\frac{1}{h\left(t,x\right)}-1\right),$ (2.44)

$V\left(t\right)=2{\displaystyle {\int }_0^{L\left(t\right)}}h\left(t,x\right)\text{d}x.$ (2.45)

None of the characteristic quantities depend on n or $a_n$. The solution depends on the two parameters $K_n$ and $\varphi$ and on $S\left(t\right)$, the fluid source at the fracture entry. The parameter $K_n$ contains all the dependence on n and $a_n$. The physical significance of $\varphi$ is

$\varphi =\frac{\text{pressure due to contact regions}}{\text{compressive normal stress at the fracture walls minus the background compressive stress}}$ (2.46)

This completes the formulation of the mathematical model of a tortuous fracture with the hyperbolic crack law.

3. Group Invariant Solution

We consider the propagation of a pre-existing hydraulic fracture with a non-zero initial length. Therefore methods used to derive a similarity solution for gravity currents [27] and for hydraulic fractures [4], both with a zero initial length, are not applicable in this work. Lie group analysis of partial differential equations, which can be used to analyze fractures with non-zero initial length, will be used to solve the problem [14] [10]. We derive the Lie point symmetries of the nonlinear diffusion Equation (2.37) which are then used to reduce the problem to a boundary value problem for an ordinary differential equation which describes the effects of surface roughness and contact regions on the evolution of a two-dimensional hydraulic fracture.

Lie Point Symmetries

The form of the Lie point symmetry generator of the nonlinear diffusion equation (2.37) is

$X={\xi }^1\left(t,x,h\right)\frac{\partial }{\partial t}+{\xi }^2\left(t,x,h\right)\frac{\partial }{\partial x}+\eta \left(t,x,h\right)\frac{\partial }{\partial h},$ (3.1)

which satisfies the determining equation

$X^{\left[2\right]}{\left[h_t-K_n\left(n{h}^{n-1}+\varphi \left(n-2\right)h^{n-3}\right)h_x^2-K_n\left(h^n+\varphi h^{n-2}\right)h_{xx}\right]|}_{\mathrm{<mo\; stretchy="false">\; (\; </mo>2.37\; <mo\; stretchy="false">\; )\; </mo>}}=\mathrm{0,}$ (3.2)

where the subscripts t and x denote partial differentiation with respect to t and x and

$X^{\left[2\right]}=X+{\zeta }_1\frac{\partial }{\partial h_t}+{\zeta }_2\frac{\partial }{\partial h_x}+{\zeta }_{11}\frac{\partial }{\partial h_{tt}}+{\zeta }_{12}\frac{\partial }{\partial h_{tx}}+{\zeta }_{22}\frac{\partial }{\partial h_{xx}}$ (3.3)

is the second prolongation of the generator X with [28]

${\zeta }_i=D_i\left(\eta \right)-h_k{D}_i\left({\xi }^k\right),i=1,2,$ (3.4)

${\zeta }_{ij}=D_j\left({\zeta }_i\right)-h_{ik}{D}_j\left({\xi }^k\right),i,j=1,2.$ (3.5)

There is summation over the repeated index k from 1 to 2 and

$D_1=D_t=\frac{\partial }{\partial t}+h_t\frac{\partial }{\partial h}+h_{tt}\frac{\partial }{\partial h_t}+h_{xt}\frac{\partial }{\partial h_x}+\cdots ,$ (3.6)

$D_2=D_x=\frac{\partial }{\partial x}+h_x\frac{\partial }{\partial h}+h_{tx}\frac{\partial }{\partial h_t}+h_{xx}\frac{\partial }{\partial h_x}+\cdots .$ (3.7)

We find that the Lie point symmetry generator X for the non-linear diffusion Equation (2.37) with $n>0$ and $\varphi >0$ is

$X=\left(c_1+c_{2}t\right)\frac{\partial }{\partial t}+\left(c_3+\frac{c_2}{2}x\right)\frac{\partial }{\partial x},=c_1{X}_1+c_2{X}_2+c_3{X}_3,$ (3.8)

where

$X_1=\frac{\partial }{\partial t},X_2=t\frac{\partial }{\partial t}+\frac{x}{2}\frac{\partial }{\partial x},X_3=\frac{\partial }{\partial x},$

and $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3$ are constants. We found that the coefficient $\eta \left(t,x,h\right)=0$. This is the first significant consequence of the from (2.24) for the diffusion coefficient.

The group invariant solution, $h=\Psi$, of the partial differential Equation (2.37) satisfies the condition

${X\left(h-\Psi \left(t,x\right)\right)|}_{h=\Psi }=0,$ (3.9)

which expands to the first order partial differential equation

$\left(c_1+c_{2}t\right)\frac{\partial \Psi }{\partial t}+\left(c_3+\frac{1}{2}{c}_{2}x\right)\frac{\partial \Psi }{\partial x}=0.$ (3.10)

Consider the general case for which $c_2\ne 0$. By solving the differential equations of the characteristic curves of (3.10) the general solution is readily derived and since $h=\Psi \left(t,x\right)$ we obtain

$h\left(t,x\right)=F\left(\xi \right),$ (3.11)

where $F\left(\xi \right)$ is an arbitrary function of the similarity variable $\xi$ given by

$\xi =\frac{2\frac{c_3}{c_2}+x}{{\left(\frac{c_1}{c_2}+t\right)}^{\frac{1}{2}}}.$ (3.12)

Unlike the group invariant solution for the linear fracture law, $h\left(t,x\right)$ depends only on the similarity variable $\xi$ and does not depend on an arbitrary parameter. Substituting the solution (3.11) and (3.12) into (2.37) reduces the second order partial differential equation to a second order ordinary differential equation for $F\left(\xi \right)$ :

$2K_n\frac{\text{d}}{\text{d}\xi }\left[\left(F^n\left(\xi \right)+\varphi F^{n-2}\left(\xi \right)\right)\frac{\text{d}F}{\text{d}\xi }\right]+\xi \frac{\text{d}F}{\text{d}\xi }=0.$ (3.13)

Since the constant $c_3$ does not appear in the ordinary differential Equation (3.13) we choose $c_3=0$ so that $\xi =0$ when $x=0$. Thus

$\xi =\frac{x}{{\left(\frac{c_1}{c_2}+t\right)}^{\frac{1}{2}}}.$ (3.14)

Consider next the boundary conditions. At the fracture tip the half-width of the fracture vanishes and (2.38) is satisfied. Expressed in terms of the group invariant solution (3.11), (2.38) becomes

$F\left(w\left(t\right)\right)=0$ (3.15)

where

$w\left(t\right)=\frac{L\left(t\right)}{{\left(\frac{c_1}{c_2}+t\right)}^{\frac{1}{2}}}.$ (3.16)

Differentiating (3.15) with respect to “t” gives

$\frac{\text{d}F}{\text{d}w}\frac{\text{d}w}{\text{d}t}=0.$ (3.17)

We consider the general case for which $F\left(w\right)$ is not a constant function. It follows that $\text{d}F/\text{d}w\ne 0$ and therefore

$\frac{\text{d}w}{\text{d}t}=\mathrm{0,}$ (3.18)

which implies that

$w\left(t\right)=w_o\mathrm{,}$ (3.19)

where $w_o$ is a constant. It follows that

$\frac{L\left(t\right)}{{\left(\frac{c_1}{c_2}+t\right)}^{\frac{1}{2}}}=w_o\mathrm{,}$ (3.20)

which can be written as

$L\left(t\right)=w_o{\left(\frac{c_1}{c_2}\right)}^{\frac{1}{2}}{\left(1+\frac{c_2}{c_1}t\right)}^{\frac{1}{2}},$ (3.21)

where we have assumed that $c_1\ne 0$. Since the characteristic length of the fracture is the initial length of the fracture, $L\left(0\right)=1$ and therefore

$w_o{\left(\frac{c_1}{c_2}\right)}^{\frac{1}{2}}=1.$ (3.22)

Thus

$L\left(t\right)={\left(1+\frac{c_2}{c_1}t\right)}^{\frac{1}{2}}$ (3.23)

and the boundary condition at the fracture tip (3.15) becomes

$F\left({\left[\frac{c_2}{c_1}\right]}^{\frac{1}{2}}\right)=0.$ (3.24)

Expressed in terms of $F\left(\xi \right)$ the boundary condition (2.39) becomes

$\left(F^n\left(\xi \right)+\varphi F^{n-2}\left(\xi \right)\right){\frac{\text{d}F}{\text{d}\xi }|}_{\xi ={\left(\frac{c_2}{c_1}\right)}^{\frac{1}{2}}}=\mathrm{0,}$ (3.25)

and the condition (2.41) becomes

$F\left(0\right)=\beta .$ (3.26)

Writing the source boundary condition (2.40) in terms of the group invariant solution (3.11) gives

$S\left(t\right)=-2K_n{\left(\frac{c_2}{c_1}\right)}^{\frac{1}{2}}\left(F^n\left(0\right)+\varphi F^{n-2}\left(0\right)\right)\frac{\text{d}F}{\text{d}\xi }\left(0\right){\left(1+\frac{c_2}{c_1}t\right)}^{-\frac{1}{2}}.$ (3.27)

Thus

$S\left(t\right)=A{\left(1+\frac{c_2}{c_1}t\right)}^{-\frac{1}{2}}$ (3.28)

where

$A=-2K_n{\left(\frac{c_2}{c_1}\right)}^{\frac{1}{2}}\left(F^n\left(0\right)+\varphi F^{n-2}\left(0\right)\right)\frac{\text{d}F}{\text{d}\xi }\left(0\right).$ (3.29)

We see that the time dependence of fluid flux into the fracture at the entry, $S\left(t\right)$, cannot be prescribed arbitrarily but must be of the form (3.28) determined by the invariant solution. The constants A and $c_2/c_1$ are prescribed if no restrictions are placed on them by the invariant numerical solution. It follows from (3.26) and (3.29) that

$\frac{\text{d}F}{\text{d}\xi }\left(0\right)=-\frac{1}{2K_n}{\left(\frac{c_1}{c_2}\right)}^{\frac{1}{2}}\frac{A}{{\beta }^n+\varphi {\beta }^{n-2}}.$ (3.30)

Since $\beta$ and $\varphi$ are prescribed, (3.30) is a boundary condition on $\text{d}F/\text{d}\xi$ at $\xi =0$. Expressed in terms of $F\left(\xi \right)$, the dimensionless fluid variables (2.42) to (2.45) become

$Q\left(t,x\right)=-\frac{2K_n}{L\left(t\right)}{\left(\frac{c_2}{c_1}\right)}^{\frac{1}{2}}\left(F^n\left(\xi \right)+\varphi F^{n-2}\left(\xi \right)\right)\frac{\text{d}F}{\text{d}\xi },$ (3.31)

${\bar{v}}_x\left(t,x\right)=-\frac{K_n}{L\left(t\right)}{\left(\frac{c_2}{c_1}\right)}^{\frac{1}{2}}\left(F^{n-1}\left(\xi \right)+\varphi F^{n-3}\left(\xi \right)\right)\frac{\text{d}F}{\text{d}\xi },$ (3.32)

$P\left(t,x\right)=F\left(\xi \right)-\varphi \left(\frac{1}{F\left(\xi \right)}-1\right),$ (3.33)

$V\left(t\right)=2{\left(\frac{c_1}{c_2}\right)}^{\frac{1}{2}}L\left(t\right){\displaystyle {\int }_0^{{\left(\frac{c_2}{c_1}\right)}^{\frac{1}{2}}}}F\left(\xi \right)\text{d}\xi .$ (3.34)

We simplify the formulation by transforming from variables $\xi$ and $F\left(\xi \right)$ to variables u and $f\left(u\right)$ where

$u={\left(\frac{c_1}{c_2}\right)}^{\frac{1}{2}}\xi =\frac{x}{L\left(t\right)},f\left(u\right)={\left(\frac{c_1}{c_2}\right)}^{\frac{1}{n}}F\left(\xi \right).$ (3.35)

We also define

$B=\frac{c_1}{c_2}.$ (3.36)

The boundary value problem transforms to the ordinary differential equation

$2K_n\frac{\text{d}}{\text{d}u}\left[\left(f^n\left(u\right)+\varphi B^{\frac{2}{n}}{f}^{n-2}\left(u\right)\right)\frac{\text{d}f}{\text{d}u}\right]+u\frac{\text{d}f}{\text{d}u}=0,$ (3.37)

subject to the boundary conditions

$f\left(0\right)=\beta B^{\frac{1}{n}},0<\beta <1,$ (3.38)

$\frac{\text{d}f}{\text{d}u}\left(0\right)=-\frac{B^{\frac{1}{n}}}{2K_n{\beta }^n}\frac{A}{1+\frac{\varphi }{{\beta }^2}},$ (3.39)

$f\left(1\right)=0,$ (3.40)

$\left(f^n\left(u\right)+\varphi B^{\frac{2}{n}}{f}^{n-2}\left(u\right)\right){\frac{\text{d}f}{\text{d}u}|}_{u=1}=0,$ (3.41)

where the fluid flux source at the fracture entry is

$S\left(t\right)=A{\left(1+\frac{t}{B}\right)}^{-\frac{1}{2}}.$ (3.42)

The remaining physical variables become

$L\left(t\right)={\left(1+\frac{t}{B}\right)}^{\frac{1}{2}}\mathrm{,}$ (3.43)

$h\left(t\mathrm{,}u\right)=B^{-\frac{1}{n}}f\left(u\right)\mathrm{,}$ (3.44)

$Q\left(t,u\right)=-\frac{2K_n}{B^{\frac{n+1}{n}}}\frac{1}{L\left(t\right)}\left[f^n\left(u\right)+\varphi B^{\frac{2}{n}}{f}^{n-2}\left(u\right)\right]\frac{\text{d}f}{\text{d}u},$ (3.45)

${\bar{v}}_x\left(t,u\right)=-\frac{K_n}{B}\frac{1}{L\left(t\right)}\left[f^{n-1}\left(u\right)+\varphi B^{\frac{2}{n}}{f}^{n-3}\left(u\right)\right]\frac{\text{d}f}{\text{d}u},$ (3.46)

$P\left(t,u\right)=B^{-\frac{1}{n}}f\left(u\right)-\varphi B^{\frac{1}{n}}\left(\frac{1}{f\left(u\right)}-B^{-\frac{1}{n}}\right),$ (3.47)

$V\left(t\right)=2B^{-\frac{1}{n}}L\left(t\right){\displaystyle {\int }_0^1}f\left(u\right)\text{d}u.$ (3.48)

The Lie point symmetry which generates the invariant solution is

$X=\left(B+t\right)\frac{\partial }{\partial t}+\frac{1}{2}x\frac{\partial }{\partial x}\mathrm{.}$ (3.49)

The physical significance of the source (3.42) is that the fluid pressure remains constant at the fracture entry for all time $t\ge 0$. From (3.47) and the boundary condition (3.38),

$P\left(t,0\right)=\beta -\varphi \left(\frac{1}{\beta }-1\right),t\ge 0.$ (3.50)

From (3.44) and (3.38), it is clear that the half-width of the fracture also remains constant at the fracture entry for all time $t\ge 0$ :

$h\left(t,0\right)=\beta ,t\ge 0.$ (3.51)

4. Asymptotic Solution

We now derive the asymptotic solution of the ordinary differential Equation (3.37) in an ε-neighbourhood of the fracture tip where $\epsilon \ll 1$. This asymptotic solution will be used in the numerical solution to offset the boundary condition at the fracture tip from $u=1$ to $u=1-\delta$, in order to avoid singularities at the fracture tip ( $u=1$ ) [14] [15].

We consider an asymptotic solution of the form

$f\left(u\right)\sim C{\left(1-u\right)}^p\mathrm{,}\text{as}u\to \mathrm{1,}$ (4.1)

where p and C are constants to be determined with $p>0$ to satisfy the boundary condition $f\left(1\right)=0$ and $C>0$ since $f\left(u\right)>0$ for $0\le u<1$. While the resulting asymptotic solution is expected to satisfy the boundary condition at the fracture tip, (3.40), it is not expected to satisfy the boundary conditions at the fracture entry.

Substituting (4.1) into the ordinary differential Equation (3.37) and dividing the result by Cp gives:

$\begin{array}{l}2K_n{C}^n\left[\left(n+1\right)p-1\right]{\left(1-u\right)}^{\left(n+1\right)p-2}\\ +2K_n\varphi C^{n-2}{B}^{\frac{2}{n}}\left[\left(n-1\right)p-1\right]{\left(1-u\right)}^{\left(n-1\right)p-2}\\ +{\left(1-u\right)}^p-{\left(1-u\right)}^{p-1}\sim \mathrm{0,}u\to 1.\end{array}$ (4.2)

The dominant terms as $u\to 1$ are the terms with the smallest values for the exponent of $\left(1-u\right)$. To obtain an asymptotic solution these terms must balance and the remaining terms must vanish as $u\to 1$. The dominant terms in (4.2) are the second and the fourth terms with exponents $\left(n-1\right)p-2$ and $p-1$, respectively. Equating these exponents gives

$p=\frac{1}{n-2},n\ne 2.$ (4.3)

The special case of $n=2$ will be examined later. For $p>0$ we require $n>2$. Equating coefficients of the second and fourth terms gives

$C={\left[\frac{n-2}{2K_n\varphi B^{\frac{2}{n}}}\right]}^{\frac{1}{n-2}}\mathrm{,}$ (4.4)

and Equation (4.2) reduces to

$\frac{6K_n{C}^n}{n-2}{\left(1-u\right)}^{\frac{5-n}{n-2}}+{\left(1-u\right)}^{\frac{1}{n-2}}\sim \mathrm{0,}u\to \mathrm{1,}$ (4.5)

for $n>2$. The second term will vanish as $u\to 1$ since $n>2$ but the first term will vanish only if $2<n<5$. It follows that the asymptotic solution for the ordinary differential Equation (3.37) is

$2<n<5:f\left(u\right)\sim {\left[\frac{n-2}{2K_n\varphi B^{\frac{2}{n}}}\right]}^{\frac{1}{n-2}}{\left(1-u\right)}^{\frac{1}{n-2}}\mathrm{,}u\to 1.$ (4.6)

The asymptotic behaviour close to the fracture tip described by (4.6) is valid for all solutions of the hyperbolic hydraulic fracture described in this work. Since

$\begin{array}{l}\left(f^n\left(u\right)+\varphi B^{\frac{2}{n}}{f}^{n-2}\left(u\right)\right)\frac{\text{d}f}{\text{d}u}\\ \sim -\frac{C^{n+1}}{n-2}\left[{\left(1-u\right)}^{\frac{3}{n-2}}+\frac{\varphi B^{\frac{2}{n}}}{C^2}{\left(1-u\right)}^{\frac{1}{n-2}}\right]\text{as}u\to \mathrm{1,}\end{array}$ (4.7)

it follows that the boundary condition (3.41) at the fracture tip, $u=1$, is identically satisfied.

Consider now the special case of $n=2$. When $n=2$, (4.2) becomes

$\begin{array}{l}2K_2{C}^2\left[3p-1\right]{\left(1-u\right)}^{3p-2}+2K_2\varphi B\left[p-1\right]{\left(1-u\right)}^{p-2}\\ +{\left(1-u\right)}^p-{\left(1-u\right)}^{p-1}\sim \mathrm{0,}u\to 1.\end{array}$ (4.8)

The second term can no longer be balanced with the fourth term and it cannot be balanced with the first term for $p>0$. We balance the first term with the fourth term to give $p=1/2$. Equation (4.8) becomes

$K_2{C}^2{\left(1-u\right)}^{-\frac{1}{2}}-K_2\varphi B{\left(1-u\right)}^{-\frac{3}{2}}+{\left(1-u\right)}^{\frac{1}{2}}-{\left(1-u\right)}^{-\frac{1}{2}}\sim \mathrm{0,}u\to 1.$ (4.9)

Since we are considering $B\ne 0$, condition (4.9) requires

$\varphi =0,C=\frac{1}{\sqrt{K_2}}.$ (4.10)

Thus we have the result

$n=2:\varphi =\mathrm{0,}f\left(u\right)\sim \frac{1}{\sqrt{K_2}}{\left(1-u\right)}^{\frac{1}{2}}\text{as}u\to 1.$ (4.11)

which is the asymptotic solution for the linear crack law with $n=2$ [14].

Since the asymptotic solution for the hyperbolic crack law exists only for $2<n<5$ we conjecture that the solution of the boundary value problem, (3.37) to (3.41), exists only for $2<n<5$. This conjecture will be investigated further when the numerical solution is derived in Section 5.

Consider now the behaviour of $\text{d}f/\text{d}u$ as $u\to 1$. From (4.6) for $2<n<5$,

$\frac{\text{d}f}{\text{d}u}\left(u\right)\sim -\frac{1}{n-2}{\left[\frac{n-2}{2K_n\varphi B^{\frac{2}{n}}}\right]}^{\frac{1}{n-2}}{\left(1-u\right)}^{\frac{3-n}{n-2}}\text{as}u\to \mathrm{1,}$ (4.12)

and therefore

$\frac{\text{d}f}{\text{d}u}=(\begin{array}{ll}0,\hfill & 2<n<3,\hfill \\ -\frac{1}{2K_3\varphi B^{\frac{2}{3}}},\hfill & n=3,\hfill \\ -\infty ,\hfill & 3<n<5,\hfill \end{array}\text{as}u\to 1.$ (4.13)

The asymptotic results, (4.6) and (4.13), will be used in the numerical solution in Section 5 to offset the boundary condition at $u=1$ to $u=1-\delta$ where $0<\delta \ll 1$.

Since $u=x/L\left(t\right)$ the half-width of the fracture (3.44), with $f\left(u\right)$ given by (4.6), is

$h\left(t\mathrm{,}x\right)\sim B^{-\frac{1}{n}}{\left[\frac{n-2}{2K_n\varphi B^{\frac{2}{n}}}\right]}^{\frac{1}{n-2}}{\left(1-\frac{x}{L\left(t\right)}\right)}^{\frac{1}{n-2}}\mathrm{,}$ (4.14)

as $x\to L\left(t\right)$, provided $2<n<5$. Differentiating (4.14) with respect to x and evaluating the limit of the result as $x\to L\left(t\right)$ gives

$\frac{\partial h}{\partial x}=(\begin{array}{ll}0,\hfill & 2<n<3,\hfill \\ -\frac{1}{2K_n\varphi B}\frac{1}{L\left(t\right)},\hfill & n=3,\hfill \\ -\infty ,\hfill & 3<n<5,\hfill \end{array}\text{as}x\to L\left(t\right).$ (4.15)

It is clear that the spatial gradient of the half-width of a two-dimensional hyperbolic model fracture is singular at $x=L\left(t\right)$ for $3<n<5$ and the thin fluid film approximation breaks down at the fracture tip. However, we see that for the hyperbolic crack law, the spatial gradient of the half-width at the fracture tip is not singular for $2<n\le 3$. In the elasticity problem for an infinite body with a plane cut subject to symmetric loads Barenblatt describes similar behaviour near the cut tip [29].

5. Numerical Solution

In this Section, we present the numerical solution for the boundary value problem (3.37) to (3.42). The boundary value problem is solved with the backward shooting method from the neigbourhood of the fracture tip, $u=1-\delta$ where $\delta \ll 1$, to the fracture entry, $u=0$. Formulating a numerical scheme so that its step-size is adjusted according to the behaviour of the solution and the resulting truncating error in each step can greatly improve the accuracy of numerical solutions [30]. The boundary value problem presented in this work was solved using the Matlab ode45 solver, which applies to the problem the Runge Kutta method of order 4 and 5 with adaptive step-size [31]. A brief algorithm that shows how the problem is solved is provided below:

1) Specify parameters $\left(K_n\mathrm{,}\beta \mathrm{,}\varphi \mathrm{,}n\right)$,

2) Choose an initial guess for the shooting parameter B,

3) Solve the ODE (3.37) subject to the initial conditions $f_{asymp}\left(1-\delta \right)$ and ${\frac{\text{d}{f}_{asymp}}{\text{d}u}|}_{u=1-\delta }$,

4) While $\left|f_N\left(0\right)-f_T\left(0\right)\right|>tol$

Adjust B using the bisection method

Solve the ODE (37) subject to the initial conditions $f_{asymp}\left(1-\delta \right)$ and ${\frac{\text{d}{f}_{asymp}}{\text{d}u}|}_{u=1-\delta }$,

end

where $f_{asymp}\left(u\right)$ is the asymptotic solution (4.6), $\text{d}{f}_{asymp}\left(u\right)/\text{d}u$ is its derivative (4.12), $0<\delta \ll 1$, $f_N\left(u\right)$ is the numerical solution obtained, $f_T\left(u\right)$ is the moving target (3.38) due to the changing parameter B and the tolerance within which a solution is accepted is $tol={10}^{-3}$.

From steps 1 and 2 of the algorithm, it is clear that all parameters in the problem are specified. We fix $K_n=1$ for all values of n although $K_n$ depends on n. We consider analysis of partially open fractures, therefore we choose the parameter $\beta$ in the range $\beta \in \left[\mathrm{0\; <mo>\; ;\; </mo>1}\right)$. We expect the pressure difference at the fracture entry, $P\left(t\mathrm{,0}\right)$, given by Equation (3.50) to be positive, therefore

$\varphi <\frac{{\beta }^2}{1-\beta }.$ (5.1)

We let

$\varphi =\frac{{\beta }^2}{1-\beta }-\gamma ,$ (5.2)

where $\gamma <{\beta }^2/\left(1-\beta \right)$. A random initial guess of the shooting parameter $B>0$ is made.

In steps 3 and 4 of the algorithm the second order ODE is written as the system of two first order ordinary differential equations.

$\frac{\text{d}{f}_1}{\text{d}u}=f_2\left(u\right)\mathrm{,}$ (5.3)

$\frac{\text{d}{f}_2}{\text{d}u}=-\frac{1}{f_1^2+\varphi B^{\frac{2}{n}}}\left[\frac{u}{2K_n}\frac{f_2}{f_1^{n-2}}+n{f}_1{f}_2^2+\left(n-2\right)\varphi B^{\frac{2}{n}}\frac{f_2^2}{f_1}\right]$ (5.4)

where $f_1\left(u\right)=f\left(u\right)$. The boundary value problem is solved as an initial value problem by shooting a solution from the fracture tip to the fracture entry. The differential Equations (5.3) and (5.4) however have singularities at the fracture tip. We use (4.6) and (4.12) for $f_1\left(u\right)$ and $f_2\left(u\right)$ as $u\to 1$. From (4.13) we see that (4.12) is singular at $u=1$ for $3<n<5$. Thus (5.3) is singular at $u=1$ for $3<n<5$. Also in (5.4) we find that

$\frac{f_2\left(u\right)}{f_1^{n-2}\left(u\right)}\sim {\left(1-u\right)}^{\frac{5-2n}{n-2}}\text{as}u\to \mathrm{1,}$ (5.5)

which is singular at $u=1$ for $2.5<n<5$,

$f_1\left(u\right)f_2^2\left(u\right)\sim {\left(1-u\right)}^{\frac{7-2n}{n-2}}\text{as}u\to \mathrm{1,}$ (5.6)

which is singular at $u=1$ for $3.5<n<5$ and

$\frac{f_2^2\left(u\right)}{f_1\left(u\right)}\sim {\left(1-u\right)}^{\frac{5-2n}{n-2}}\text{as}u\to \mathrm{1,}$ (5.7)

which is singular at $u=1$ for $2.5<n<5$. It follows that the system of equations (5.3) and (5.4) is singular for $2.5<n<5$ as $u\to 1$. While there are no singularities in $\text{d}f/\text{d}{f}_1$ and $\text{d}f/\text{d}{f}_2$ for $2<n\le 2.5$, the algorithm is unable to solve the problem when the initial condition $f_1\left(1\right)=0$ is used in the entire range $2<n<5$. We therefore offset the boundary conditions at the fracture tip from $u=1$ to $u=1-\delta$ for the entire range $2<n<5$. Therefore the system of Equations (5.3) and (5.4) is solved by shooting a solution from $u=1-\delta$, where $0<\delta \ll 1$, to the fracture entry. The initial conditions on $f_1\left(u\right)$ and $f_2\left(u\right)$ are evaluated at $u=1-\delta$ using (4.6) and (4.12):

$f_1\left(1-\delta \right)={\left[\frac{n-2}{2K_n\varphi B^{\frac{2}{n}}}\right]}^{\frac{1}{n-2}}{\delta }^{\frac{1}{n-2}}\mathrm{,}$ (5.8)

$f_2\left(1-\delta \right)=-\frac{1}{n-2}{\left[\frac{n-2}{2K_n\varphi B^{\frac{2}{n}}}\right]}^{\frac{1}{n-2}}{\delta }^{\frac{3-n}{n-2}}\mathrm{.}$ (5.9)

The shooting method is applied iteratively, with the parameter B updated at each iteration by the bisection method, until the boundary condition at the fracture entry (3.38),

$f_1\left(0\right)=\beta B^{\frac{1}{n}}\mathrm{,}$ (5.10)

is satisfied. The boundary condition (3.41) is identically satisfied while the value of the constant A in the source flux of fluid at the entry (3.42) can now be determined from the boundary condition (3.39),

$f_2\left(0\right)=-\frac{B^{\frac{1}{n}}A}{2K_n{\beta }^n\left(1+\frac{\varphi }{{\beta }^2}\right)}\mathrm{.}$ (5.11)

The numerical solution for $f\left(u\right)$, obtained from solving the boundary value problem (3.37) to (3.42), is substituted into equations (3.43) to (3.48) which give the properties of the fracture as well as the behaviour of the fluid inside the fracture.

We now analyze the half-width of the model fracture $h\left(t\mathrm{,}x\right)$. Figure 2 shows the half-width $h\left(t\mathrm{,}x\right)$ of a partially open model fracture plotted against x for a range of dimensionless times and for $n=4$, $n=3$ and $n=2.5$. From (3.50) a tortuous fracture with the hyperbolic crack law admits only one working condition of constant pressure difference at the fracture entry. Since the far-field compressive stress $-{\sigma }_{zz}^{\left(\infty \right)}$ is also constant this implies for the hyperbolic fracture the pressure of fluid injection at the fracture entry, $p\left(t\mathrm{,0}\right)$, is constant. Due to the PKN approximation the half-width of the model fracture is also constant for $t\ge 0$ at the fracture entry and is given by (3.51). This is clearly shown in Figure 2 where $\beta =0.5$. Both the half-width $h\left(t\mathrm{,}x\right)$ for $x>0$ and the length $L\left(t\right)$ of the model fracture grow as time increases resulting in the volume of the model fracture per unit breadth, $V\left(t\right)$, which from (2.45) is the area under the half-width curve, growing with increasing time t.

The behaviour of $\frac{\partial h}{\partial x}$ at the fracture tip given by (4.15) is clearly illustrated in Figure 2. The chosen values, $n=4$, $n=3$ and $n=2.5$, lie in the ranges $3<n<5$, $n=3$ and $2<n<3$ and we see that $\frac{\partial h}{\partial x}$ is negative infinity, finite and zero at the tip.

The length, $L\left(t\right)$, of a partially open fracture is plotted against time for a range of values of $\beta$ in Figure 3 and for a range of values of $\varphi$ in Figure 4. The parameter $\beta$ describes the extent to which the fracture si partially open while $\varphi$ compares the contribution of the pressure due to the contact regions with the pressure due to the fluid in the fracture. From Figure 3 we observe that for a prescribed value of $\varphi$ as the value of $\beta$ increases, the fracture becomes longer. From Figure 3 and Figure 4 it is clear that for fixed values of $\beta$ and $\varphi$, tortuosity increases, that is as n decreases from 5 to 2, the length of the model fracture increases. We also see from Figure 4 that there is less change in the length due to change in $\varphi$ as the tortuosity increases.

The numerical solution of the boundary value problem (3.37) to (3.42) was obtained for a range of values of $\beta$, $\varphi$ and n with the diffusion constant prescribed as $K_n=1$. The values obtained by the shooting method for A and B in the source flux of fluid at the entry, (3.42), are presented in Table 1. We see from Table 1 that as n decreases, that is as tortuosity increases, for fixed values of $\beta$ and $\varphi$, the constant A increases. Thus a stronger fluid flux at the fracture entry, $S\left(t\right)$, is required to propagate a very tortuous fracture. From (3.43), the speed of propagation of the fracture is

$\frac{\text{d}L}{\text{d}t}=\frac{1}{2B}{\left(1+\frac{t}{B}\right)}^{-\frac{1}{2}}\mathrm{.}$ (5.12)

(a)(b)(c)

Figure 2. Numerical solution for the half-width $h\left(t\mathrm{,}x\right)$ of a partially open model fracture, $\beta =0.5$ with $\varphi =0.45$, plotted against x with the diffusion constant $K_n=1$ and for increasing values of time t with (a) $n=4$, (b) $n=3$, (c) $n=2.5$.

(a)(b)(c)

Figure 3. Numerical solution of the length $L\left(t\right)$ of a partially open model fracture plotted against time t for $\beta =0.5$, $\beta =0.5$, $\beta =0.7$, with $\varphi =0.45$, diffusion constant $K_n=1$ and (a) $n=4$, (b) $n=3$, (c) $n=2.5$.

(a)(b)(c)

Figure 4. Numerical solution of the length of the fracture $L\left(t\right)$ of a partially open model fracture plotted against time t for $\varphi =0.1$, $\varphi =0.3$, $\varphi =0.45$, with $\beta =0.5$, diffusion constant $K_n=1$ and (a) $n=4$, (b) $n=3$, (c) $n=2.5$.

It is clear that the constant B determines the magnitude of the speed of propagation of the fracture. From Table 1, for fixed values of $\beta$ and $\varphi$, the constant B decreases as n decreases. Thus the speed of propagation of the fracture increases as tortuosity of the fracture increases. This agrees with Figure 3 in which $L\left(t\right)$ is plotted against time t for a range of values of n, $\beta$ and $\varphi$.

We also see from Table 1 that for fixed values of n and $\varphi$, A decreases as $\beta$ decreases and therefore the flux of fluid at the entry is weaker when the width of the partially open fracture is smaller. For fixed values of n and $\varphi$, the constant B increases as $\beta$ decreases and therefore the speed of propagation of the fracture decreases for partially open fractures of smaller width, in agreement with Figure 3. For fixed values of $\beta$ and n, A decreases as $\varphi$ decreases and therefore the fluid flux at the entry is weaker as the pressure ratio $\varphi$ decreases, consistent with Figure 5. Also for fixed values of $\beta$ and n, B increases as $\varphi$ decreases and therefore the speed of propagation of the fracture decreases as $\varphi$ decreases, consistent with Figure 4 for $L\left(t\right)$.

6. Width Averaged Fluid Velocity

The ratio of the width averaged fluid velocity, ${\bar{v}}_x$, given by (3.46) to the speed of propagation of the fracture tip,

$\frac{\text{d}L}{\text{d}t}=\frac{1}{2B}{\left(1+\frac{t}{B}\right)}^{-\frac{1}{2}}=\frac{1}{2BL\left(t\right)}$ (6.1)

is

$\frac{{\bar{v}}_x\left(t,u\right)}{\text{d}L/\text{d}t}=-2K_n\left(f^2\left(u\right)+\varphi B^{\frac{2}{n}}\right)f^{n-3}\left(u\right)\frac{\text{d}f}{\text{d}u}.$ (6.2)

Using the asymptotic solutions (4.6) for $f\left(u\right)$ and (4.12) for $\text{d}f/\text{d}u$ it can be verified that

${\frac{{\bar{v}}_x\left(t,u\right)}{\text{d}L/\text{d}t}|}_{u=1}=1.$ (6.3)

The width averaged fluid velocity therefore tends to the speed of propagation of the fracture tip as $u\to 1$. Using the boundary conditions (3.38) and (3.39) for $f\left(u\right)$ and $\text{d}f/\text{d}u$ at $u=0$ we find that

${\frac{{\bar{v}}_x\left(t,u\right)}{\text{d}L/\text{d}t}|}_{u=0}=\frac{AB}{\beta }\mathrm{.}$ (6.4)

The curves of the width averaged fluid velocity ratio (6.2) plotted against u for $0\le u\le 1$ therefore have end points $\left(\mathrm{0,}AB/\beta \right)$ at $u=0$ and $\left(\mathrm{1,1}\right)$ at $u=1$.

In Figure 6, the velocity ratio curves (6.2) for a partially open fracture with $\beta =0.5$ and $\varphi =0.1$, $\varphi =0.3$ and $\varphi =0.45$ respectively are plotted against u for $0\le u\le 1$, with the diffusion constant $K_n=1$ and for $n=4$, $n=3$ and $n=2.5$. The curves all pass through the point $\left(\mathrm{1,1}\right)$. In Figure 6(c) for example, the end points at $u=0$ for $\varphi =0.1$, $\varphi =0.3$ and $\varphi =0.45$ are $\left(\mathrm{0,0.5215}\right)$,

Table 1. Numerical values of the parameters A and B in the source fluid flux at the fracture entry for $K_n=1$.

(a)(b)(c)

Figure 5. Numerical solution of the source fluid flux $S\left(t\right)$ at the entry of a partially open model fracture plotted against time t for $\beta =0.7$, diffusion constant $K_n=1$ and (a) $n=4$, (b) $n=3$, (c) $n=2.5$.

(a)(b)(c)

Figure 6. Velocity ratio curves ${\bar{v}}_x\left(t,u\right)/\left(\text{d}L/\text{d}t\right)$ of a partially open model fracture, $\beta =0.5$, plotted against u with the diffusion constant $K_n=1$ and for parameter values $\varphi =0.1$, $\varphi =0.3$ and $\varphi =0.45$ respectively with (a) $n=4$, (b) $n=3$, (c) $n=2.5$.

$\left(\mathrm{0,0.4422}\right)$ and $\left(\mathrm{0,0.4272}\right)$ respectively and are in good agreement with the exact result $\left(\mathrm{0,}AB/\beta \right)$ which, with the aid of Table 1, gives $\left(\mathrm{0,0.5215}\right)$, $\left(\mathrm{0,0.4598}\right)$ and $\left(\mathrm{0,0.4453}\right)$. We see from Figure 6 that the velocity ratio decreases as $\varphi$ increases, that is, as the effect of the pressure due to the contact regions increases. The effect of touching asperities is therefore to decrease the width averaged fluid velocity along the length of the tortuous fracture. We also

see from Figure 6 that the velocity ratio ${\bar{v}}_x\left(t,u\right)/\frac{\text{d}L}{\text{d}t}$ increases approximately

linearly with u along the fracture for the values of the parameters considered. The linear approximation is best for $n=4$ and less applicable as n decreases and the fracture becomes more tortuous.

In order to derive approximate analytical solutions we will make the approximation that the variation of ${\bar{v}}_x/\frac{\text{d}L}{\text{d}t}$ along the fracture is linear in u for

$0\le u\le 1$ [10] [14] [15]. In the approximate solution, we will replace $f\left(u\right)$, A and B by $\hat{f}\left(u\right)$, $\hat{A}$ and $\hat{B}$. By using the equation of the straight line between the points $\left(\mathrm{0,}\hat{A}\hat{B}/\beta \right)$ and $\left(\mathrm{1,1}\right)$ and using Equation (6.2) we obtain

$-2K_n\left({\hat{f}}^2\left(u\right)+\varphi {\hat{B}}^{\frac{2}{n}}\right){\hat{f}}^{n-3}\left(u\right)\frac{\text{d}\hat{f}}{\text{d}u}=\left(1-\frac{\hat{A}\hat{B}}{\beta }\right)u+\frac{\hat{A}\hat{B}}{\beta }\mathrm{,}$ (6.5)

for $2<n<5$. Equation (6.5) is a first order ordinary differential equation for $\hat{f}\left(u\right)$ subject to the boundary conditions

$\hat{f}\left(0\right)=\beta {\hat{B}}^{\frac{1}{n}}\mathrm{,}$ (6.6)

$\frac{\text{d}\hat{f}}{\text{d}u}\left(0\right)=-\frac{{\hat{B}}^{\frac{1}{n}}\hat{A}}{2K_n{\beta }^{n-2}\left({\beta }^2+\varphi \right)}\mathrm{,}$ (6.7)

$\hat{f}\left(1\right)=0.$ (6.8)

By using the differential Equation (6.5) and boundary condition (6.6) we see that the boundary condition (6.7) is identically satisfied. Integrating the ordinary differential Equation (6.5) gives

$-2K_n\left(\frac{{\hat{f}}^n}{n}\left(u\right)+\frac{\varphi {\hat{B}}^{\frac{2}{n}}{\hat{f}}^{n-2}\left(u\right)}{n-2}\right)=\frac{1}{2}\left(1-\frac{\hat{A}\hat{B}}{\beta }\right)u^2+\frac{\hat{A}\hat{B}}{\beta }u+C\mathrm{,}$ (6.9)

where C is an integration constant. Applying the boundary condition (6.8) gives

$C=-\frac{1}{2}\left(1+\frac{\hat{A}\hat{B}}{\beta }\right)$ (6.10)

and therefore

$\frac{4K_n}{n\left(n-2\right)}\left[\left(n-2\right){\hat{f}}^n\left(u\right)+\varphi {\hat{B}}^{\frac{2}{n}}{\hat{f}}^{n-2}\left(u\right)\right]=\left(1-u\right)\left[\left(1-\frac{\hat{A}\hat{B}}{\beta }\right)u+1+\frac{\hat{A}\hat{B}}{\beta }\right]\mathrm{.}$ (6.11)

Finally, imposing the boundary condition (6.6) on (6.11) gives a relation between $\hat{A}$ and $\hat{B}$ for prescribed values of n, $\varphi$, $\beta$ and $K_n$ :

$\hat{A}+\frac{\beta }{\hat{B}}=\frac{4K_n{\beta }^{n-1}}{n\left(n-2\right)}\left[\left(n-2\right){\beta }^2+n\varphi \right].$ (6.12)

Equation (6.12) can be re-expressed as

$\hat{A}+\frac{\beta }{\hat{B}}=E\left(n\mathrm{,}\beta \mathrm{,}\varphi \mathrm{,}{K}_n\right)\mathrm{,}$ (6.13)

where

$E\left(n\mathrm{,}\beta \mathrm{,}\varphi \mathrm{,}{K}_n\right)=\frac{4K_n{\beta }^{n-1}}{n\left(n-2\right)}\left[\left(n-2\right){\beta }^2+n\varphi \right]\mathrm{.}$ (6.14)

If $\hat{B}$ is prescribed then

$\hat{A}=E-\frac{\beta }{\hat{B}}.$ (6.15)

The graph of $\hat{A}$ against $\hat{B}$ is illustrated in Figure 7(a). We see that $\hat{A}<E$ and that for $\hat{A}>0$, which from (3.42) describes fluid flux into the fracture at the entry ( $u=0$ ), then it is necessary that $\hat{B}>\frac{\beta }{E}$. If $\hat{A}$ is prescribed then

$\hat{B}=\frac{\beta }{E-\hat{A}}.$ (6.16)

The graph of $\hat{B}$ against $\hat{A}$ is illustrated in Figure 7(b). The inequalities $\hat{A}<E$ and $\hat{B}>\frac{\beta }{E}$ are again apparent.

From (6.11)

$\hat{f}\left(u\right)=\frac{n\left(n-2\right)}{4K_n}\left[\frac{-\left(1-\frac{\hat{A}\hat{B}}{\beta }\right)u^2-\frac{2\hat{A}\hat{B}}{\beta }u+\left(1+\frac{\hat{A}\hat{B}}{\beta }\right)}{\left(n-2\right){\hat{f}}^{n-1}\left(u\right)+n\varphi {\hat{B}}^{\frac{2}{n}}{\hat{f}}^{n-3}\left(u\right)}\right],$ (6.17)

where $2<n<5$ and $0\le u\le 1$. Equation (6.17) is implicit because $\hat{f}\left(u\right)$ appears on both sides. Therefore in order to find the approximate analytical solution $\hat{f}\left(u\right)$, a numerical method is required in which (6.17) is iterated until convergence is reached. In (6.17) the numerator has the factor $\left(1-u\right)$ and the denominator has the factor $f^{n-3}\left(u\right)$. However, the iterative method used to solve for $\hat{f}\left(u\right)$ requires that the iterated function (6.17) is neither factorised nor a power law. The parameters n, $K_n$, $\beta$ and $\varphi$ are prescribed. Unlike the numerical solution, either $\hat{A}$ or $\hat{B}$ can be prescribed. If $\hat{A}$ is prescribed then $\hat{B}$ is obtained from (6.16) while if $\hat{B}$ is prescribed then $\hat{A}$ is obtained from (6.15). In this work we prescribe $\hat{B}$ and therefore the expression $\hat{A}\hat{B}/\beta$ appearing in (6.17) is

$\frac{\hat{A}\hat{B}}{\beta }=\frac{E\left(n,\beta ,\varphi ,K_n\right)\hat{B}}{\beta }-1.$ (6.18)

For $n=2.5$, $\beta =0.5$, $\varphi =0.45$ and $K_n=1$ we found that if we let

(a)(b)

Figure 7. Graphs of parameters $\hat{A}$ and $\hat{B}$ of the approximate analytical solution with $n=2.5$, $\beta =0.5$, $\varphi =0.45$ and $K_n=1$ : (a) Graph of $\hat{A}$ against $\hat{B}$, (b) Graph of $\hat{B}$ against $\hat{A}$.

$B=\hat{B}=0.49447$ then the approximate analytical value $\hat{A}=0.403031$ is in reasonable agreement with the numerical value $A=0.450311$ since the fracture is very tortuous with $n=2.5$.

In order to solve for the approximate analytical solution $\hat{f}\left(u\right)$ we use a non-orthogonal linesearch method described in [32]. We define the function (6.17), with the aid of (6.18), as

$H\left(\hat{f}\right)=\frac{n\left(n-2\right)}{4K_n}\left[\frac{-\left(2-\frac{E\hat{B}}{\beta }\right)u^2-\frac{2\hat{A}\hat{B}}{\beta }u+\frac{E\hat{B}}{\beta }}{\left(n-2\right){\hat{f}}^{n-1}\left(u\right)+n\varphi {\hat{B}}^{\frac{2}{n}}{\hat{f}}^{n-3}\left(u\right)}\right]\mathrm{.}$ (6.19)

We then consider two coordinates, $\left({\hat{f}}_{i+1}\mathrm{,}{\hat{f}}_{i+1}\right)$ and $\left({\hat{f}}_i\mathrm{,}H\left({\hat{f}}_i\right)\right)$, with which to calculate a gradient m:

$m=\frac{{\hat{f}}_{i+1}-H\left({\hat{f}}_i\right)}{{\hat{f}}_{i+1}-{\hat{f}}_i}\mathrm{.}$ (6.20)

We re-express (6.20) to give ${\hat{f}}_{i+1}$ :

${\hat{f}}_{i+1}=\frac{m{\hat{f}}_i-H\left({\hat{f}}_i\right)}{m-1},i=1,\cdots ,k,$ (6.21)

where k is a constant. Therefore in order to find the approximate analytical solution $\hat{f}\left(u\right)$ we proceed as follows:

1) Choose a value of m.

2) Let the first guess for ${\hat{f}}_1$ be a vector given by the asymptotic solution (4.6).

3) Iterate with each solution for ${\hat{f}}_{i+1}$ used as ${\hat{f}}_i$ in the next iteration,

4) The iteration is continued until convergence is reached at the k^th step where the solution ${\hat{f}}_{i+1}={\hat{f}}_i$.

We choose $m=-1$ which in a few iterations gives an approximate analytical solution $\hat{f}\left(\hat{u}\right)$. The gradient m is required to be negative because the spatial gradient of the upper half-width of the fracture is negative.

For a very tortuous partially open fracture with $2<n\le 3$ we found that more iterations are required to obtain the approximate analytical solution for the half-width than for a less tortuous fracture with $3<n<5$. The computational time required to find the approximate analytical solution for the half-width of a partially open fracture with $\beta =0.5$, and $\varphi =0.45$ is approximately 0.7 seconds. Figure 8 shows that the graphs of the approximate analytical solution for the half-width approximately overlap the graphs of the numerical solution. We found that the approximate analytical solution is a better approximation for $n\ge 3$ than for $2<n\le 3$. We also observe from Table 2 that with a specified value of $B=\hat{B}$ the parameter A obtained by the numerical method and the parameter $\hat{A}$ obtained from Equation (6.15) are very closely matched with an error below 6% for $n\ge 3$ and an increased error below 13% for $2<n<3$. This is in agreement with Figure 6 in which the approximation that the curves are linear is clearly better for the less tortuous fracture with $n=4$ and $n=3$ than for the more tortuous fracture with $n=2.5$. Results of the approximate analytical solution and the numerical solution are in good agreement as observed both in Figure 8 and in Table 2. In the absence of an exact analytical solution, it is clear that the approximate analytical solution is a reasonable alternative.

7. Conclusions

There are several novel features in this investigation. In the mathematical modelling of the tortuous hydraulic fracture, the hyperbolic crack law was used with the PKN approximation. Previous research [2] with the hyperbolic crack law used the elasticity model with the Cautchy principal value integral (2.5) for the stress in a two-sided fracture which has the advantage that the stress intensity

(a)(b)(c)

Figure 8. Numerical solution () and approximate analytical solution () for the half-width of a partially open model fracture, $\beta =0.5$, plotted against x with

diffusion constant $K_n=1$ and for $\varphi =0.45$, with (a) $n=4$, (b) $n=3$, (c) $n=2.5$.

Table 2. Numerical values A compared to approximate analytical values $\hat{A}$ for $B=\hat{B}$ with $K_n=1$, and various values of $\beta$ and $\varphi$.

factor at the fracture tip is non-zero but is difficult to analyze mathematically. With the PKN approximation, a nonlinear diffusion equation for a one-sided fracture was derived with a diffusion coefficient (2.24) not considered before in hydraulic fracturing. Although an exact analytical solution could not be derived, an existing model was developed based on the width-averaged fluid velocity which led to a new approximate analytical solution. In the mathematical analysis of the problem, the Lie point symmetry of the nonlinear diffusion equation had the special property that the coefficient $\eta \left(t\mathrm{,}x\mathrm{,}h\right)$ vanished which had a significant effect on the solution. A shooting method was developed to solve the problem numerically and a recently established non-orthorgonal line search method was applied to obtain an iterative solution for the implicit approximate analytical solution. The following discussion gives a summary of the results obtained.

For the hyperbolic crack, law the only working condition at the fracture entry is constant pressure. From the PKN approximation, the half-width of the fracture remains constant at the fracture entry. Mathematically these results follow because in the Lie point symmetry, which generates the group invariant solution, the coefficient $\eta \left(t\mathrm{,}x\mathrm{,}h\right)$ vanishes. The exponent n in the power law in the modified Reynolds flow law lies in the finite range $2<n<5$. While the linear crack law yielded a full range of working conditions of both fluid injection and extraction at the fracture entry and $0<n<\infty$, it was observed that for a very tortuous fracture, all solutions for a range of working conditions converged towards the solution of the constant pressure working condition [14]. A brief comparison of the results for the linear crack law and the hyperbolic crack law suggests that for very tortuous fractures, the constant pressure working condition is a dominant condition at the fracture entry.

The effect of the pressure due to contact regions is described by the pressure ratio $\varphi$. The length of the fracture after time t increases as $\varphi$ increases for all values n which shows that the effect of touching asperities is to grow the length of the fracture. The effect of $\varphi$ on the length of the fracture is more prominent in a very tortuous fracture with $2<n<3$ than in a less tortuous fracture with $n\ge 3$. The width averaged fluid velocity decreases as $\varphi$ increases near the fracture entry but is almost independent of $\varphi$ as the half-width decreases near the fracture tip. This shows that the effect of touching asperities is to decrease the width averaged fluid velocity in regions where the half-width is not small.

As n decreases from $n=5$ to $n=2$, the tortuosity of the fracture increases. The value $n=3$ describes the Reynolds flow law and even then there is some surface roughness. We saw that as n decreases from $n=4$ to $n=2.5$ the length of the model fracture after time t increases due to an increase in tortuosity. The tortuosity and the pressure due to contact regions both have the effect of increasing the fracture length.

For the hyperbolic crack law the spatial gradient of the half-width is singular at the fracture tip for $3<n<5$ and the thin fluid film approximation breaks down at the tip. However, for the Reynolds flow law $n=3$ the spatial gradient of the half-width is finite while for $2<n<3$ it is zero and the thin fluid film approximation remains valid over the whole length of the fracture. As with the linear crack law, increasing tortuosity removes the singularity at the fracture tip.

The fluid flux at the fracture entry cannot be prescribed arbitrarily. It is determined by the invariant solution. The dependence on time must be of the form (3.42) and the constants A and B are determined numerically by the shooting method.

The approximate analytical solution is useful because there is no exact analytical solution for a hydraulic fracture with contact regions described by the hyperbolic crack law. There are no special cases which could lead to an analytical solution because constant pressure at the fracture entry is the only working condition. The approximate analytical solution is implicit in form and has to be solved iteratively and numerically in the final step. It is in good agreement with the numerical solution although the error increases slightly as the tortuosity increases. For the numerical solution neither of the constants, A or B, in the fluid flux at the fracture entry can be specified. For the approximate analytical solution either $\hat{A}$ can be specified and $\hat{B}$ is obtained from (6.16) or $\hat{B}$ can be specified and $\hat{A}$ is obtained from (6.15). This is because, built into the approximate analytical solution is the approximate result that the ratio of the width averaged fluid velocity to the speed of propagation of the fracture changes linearly along the fracture.

Acknowledgements

D.P.M. acknowledges financial support from the National Research Foundation, South Africa, under Grant No 132189.

M.R.R.K. acknowledges financial support from the National Research Foundation, Pretoria, South Africa.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Fitt, A.D., Mason, D.P. and Moss, E.A. (2007) Group Invariant Solution for a Pre-Existing Fluid-Driven Fracture in Impermeable Rock. Zeitschrift fur Angewandte Mathematik und Physik, 58, 1049-1067. https://doi.org/10.1007/s00033-007-7038-2
[2]	Fitt, A.D., Kelly, A.D. and Please, C.P. (1995) Crack Propagation Models for Rock Fracture in a Geothermal Energy Reservoir. SIAM Journal on Applied Mathematics, 55, 1592-1608. https://doi.org/10.1137/S0036139993260241
[3]	Geertsma, J. and de Klerk, F. (2013) A Rapid Method of Predicting Width and Extent of Hydraulic Induced Fractures. Journal of Petroleum Technology, 246, 1571-1581. https://doi.org/10.2118/2458-PA
[4]	Spence, D.A. and Sharp, P. (1985) Self-Similar Solutions for Elastohydrodynamic Cavity Flow. Proceedings of the Royal Society of London, 400, 298-313. https://doi.org/10.1098/rspa.1985.0081
[5]	Spence, D.A. and Turcotte, D.L. (1985) Magma-Driven Propagation of Cracks. Journal of Geophysical Research, 90, 575-580. https://doi.org/10.1029/JB090iB01p00575
[6]	Emerman, S., Turcotte, D.L. and Spence, D.A. (1986) Transport of Magma and Hydrothermal Solutions by Laminar and Turbulent Fluid Fracture. Physics of the Earth and Planetary Interiors, 36, 276-284.
[7]	Spence, D.A., Sharp, P. and Turcotte, D.L. (1987) Buoyancy-Driven Crack Propagation: A Mechanism for Magma Migration. Journal of Fluid Mechanics, 174, 135-153. https://doi.org/10.1017/S0022112087000077
[8]	Perkins, T. and Kern, L. (1961) Widths of Hydraulic Fractures. Journal of Petroleum Technology, 222, 937-949. https://doi.org/10.2118/89-PA
[9]	Nordgren, R. (1972) Propagation of Vertical Hydraulic Fractures. Journal of Petroleum Technology, 253, 306-314. https://doi.org/10.2118/3009-PA
[10]	Fareo, A.G. and Mason, D.P. (2011) Group Invariant Solutions for a Pre-Existing Fluid-Driven Fracture in Permeable Rock. Nonlinear Analysis: Real World Applications, 12, 767-779. https://doi.org/10.1016/j.nonrwa.2010.08.004
[11]	Fareo, A.G. and Mason, D.P. (2013) Group Invariant Solution for a Pre-Existing Fracture Driven by a Power Law Fluid in Impermeable Rock. Communications in Nonlinear Science and Numerical Simulation, 18, 3298-3316. https://doi.org/10.1016/j.cnsns.2013.04.019
[12]	Anthonyrajah, M., Mason, D.P. and Fareo, A.G. (2013) Propagation of a Turbulent Fluid Fracture. International Journal of Non-Linear Mechanics, 54, 105-114. https://doi.org/10.1016/j.ijnonlinmec.2013.03.006
[13]	Goodman, R.E. (1947) The Mechanical Properties of Joints. Proceedings 3rd Congress ISRM, Denver, 1A, 127-140.
[14]	Kgatle, M.R.R. and Mason, D.P. (2014) Propagation of a Linear Two-Dimensional Hydraulic Fracture with Tortuosity. International Journal of Non-Linear Mechanics, 63, 39-53. https://doi.org/10.1016/j.ijnonlinmec.2014.01.013
[15]	Kgatle, M.R.R. and Mason, D.P. (2016) Linear Hydraulic Fracture with Tortuosity: Conservation Laws and Fluid Extraction. International Journal of Non-Linear Mechanics, 79, 11-25. https://doi.org/10.1016/j.ijnonlinmec.2015.10.014
[16]	Pine, R.J. and Cundall, P.A. (1985) Applications of the Fluid-Rock Interaction Programme (FRIP) to the Modelling of Hot Dry Geothermal Energy Systems. Proceedings of the International Symposium on Fundamentals of Rock Joints, Björkliden, 15-20 September 1985, 293-302.
[17]	Ge, S. (1997) A Governing Equation for Fluid Flow in Rough Fractures. Water Resources Research, 33, 53-61. https://doi.org/10.1029/96WR02588
[18]	Oron, A.P. and Berkowitz, B. (1998) Flow in Rock Fractures: The Local Cubic Law Assumption Re-Examined. Water Resources Research, 34, 2811-2825. https://doi.org/10.1029/98WR02285
[19]	Qian, J., Chen, Z., Zhan, H. and Guan, H. (2011) Experimental Study of the Effect of Roughness and Reynolds Number on Fluid Flow in Rough-Walled Single Fractures: A Check of Local Cubic Law. Hydrological Processes, 25, 614-622. https://doi.org/10.1002/hyp.7849
[20]	Witherspoon, P.A., Wang, J.S.Y., Iwai, K. and Gale, J.E. (1980) Validity of Cubic Law for Fluid Flow in a Deformable Rock Fracture. Water Resources Research, 16, 1016-1024. https://doi.org/10.1029/WR016i006p01016
[21]	Bandis, S., Lumsden, A.C. and Barton, N.R. (1983) Fundamentals of Rock Joint Deformation. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 20, 249-268. https://doi.org/10.1016/0148-9062(83)90595-8
[22]	Murphy, H., Dash, Z., Zyvoloski, G., White, A. and Wing, M. (1983) Fundamentals of Rock Joint Deformation. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 20, 249-268. https://doi.org/10.1016/0148-9062(83)91304-9
[23]	Cook, N.G.W. (1992) Natural Joints in Rock: Mechanical, Hydraulic and Seismic Behaviour and Properties under Normal Stress. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 29, 198-223. https://doi.org/10.1016/0148-9062(92)93656-5
[24]	Adachi, J.I. and Peirce, A.P. (2008) Asymptotic Analysis of an Elasticity Equation for a Finger-Like Hydraulic Fracture. Journal of Elasticity, 90, 43-69. https://doi.org/10.1007/s10659-007-9122-4
[25]	Vazquez, J.L. (2017) The Porous Medium Equation: Chapter 21—Further Applications. Oxford University Press, Oxford.
[26]	King, J.R. and Bowen, M. (2001) Moving Boundary Problems and Non-Uniqueness for Thin Film Equations. European Journal of Applied Mathematics, 12, 321-356. https://doi.org/10.1017/S0956792501004405
[27]	Huppert, H.E. (1982) The Propagation of Two-Dimensional and Axisymmetric Viscous Gravity Currents over a Rigid Horizontal Surface. Journal of Fluid Mechanics, 121, 43-58. https://doi.org/10.1017/S0022112082001797
[28]	Ibragimov, N.H. and Anderson, R.L. (1994) One-Parameter Transformation Groups. In: Ibragimov, N.H., Ed., CRC Handbook of Lie Group Analysis of Differential Equations, Volume 1: Symmetries, Exact Solutions and Conservation Laws, CRC Press, Boca Raton, 7-14.
[29]	Barenblatt, G.I. (2014) Flow Deformation and Fracture. Cambridge University Press, Cambridge, Chapter 6.
[30]	Press, W.H. and Teukolsky, S.A. (1992) Adaptive Step Size Runge-Kutta Integration. Computers in Physics, 6, 188. https://doi.org/10.1063/1.4823060
[31]	Dormand, J.R. and Prince, P.J. (1980) A Family of Embedded Runge-Kutta Formulae. Computational and Applied Mathematics, 6, 19-26. https://doi.org/10.1016/0771-050X(80)90013-3
[32]	Pallares, M.M.R. and Rodriguez, C.W. (2014) Methodology for Solving the Divergence of Fixed Point Method for the Solution of a Nonlinear Equation. ARPN Journal of Engineering and Applied Sciences, 9, 2177-2182.