Remarks on application of different variables for the PKN model of hydrofracturing: various fluid-flow regimes

Abstract

The problem of hydraulic fracture for the PKN model is considered within the framework presented recently by Linkov (Doklady Phys 56(8):436–438, 2011). The modified formulation is further enhanced by employing an improved regularized boundary condition near the crack tip. This increases solution accuracy especially for singular leak-off regimes. A new dependent variable having clear physical sense is introduced. A comprehensive analysis of numerical algorithms based on various dependent variables is provided. Comparison with know numerical results has been given.

Appendices

Appendices

1.1 A Carter’s leak-off function in the normalised formulation

Consider the transformation of the Carter law described by (4) when applying the normalization (18). Assume that:

$$\begin{aligned} \frac{1}{\sqrt{t-\tau (x)}}=\frac{D(t)}{\sqrt{1-\tilde{x}}}+R(t,\tilde{x}), \end{aligned}$$

(74)

where function $D(t)$ is defined in (22) while the remainder $R$ is estimated later in (77).

To find function $D(t)$, and thus to obtain an exact form of Eq. (22), it is enough to compute the limit

$$\begin{aligned} D^2(t)=\lim _{\tilde{x}\rightarrow 1}\frac{1-\tilde{x}}{t-\tau (x)}. \end{aligned}$$

(75)

This can be done by utilising L’Hopital’s rule with taking into account that $x\rightarrow L(t)$ as $\tilde{x}\rightarrow 1$,

$$\begin{aligned} \tau (x)=\tau \left( L(t)\tilde{x}\right) =L^{-1}(L(t)\tilde{x}), \end{aligned}$$

(76)

and that the crack length is a smooth function of time ($L\in C^1$ at least). The last fact immediately follows from the problem formulation in terms of evolution system (32).

Having the value of $D(t)$ we can estimate the remainder $R(t,\tilde{x})$ when $\tilde{x}\rightarrow 1$, or, what it is equivalent to when $x\rightarrow l(t)$ (or $t\rightarrow \tau (x)$). For this reason, we search for a parameter $\xi \ne 0$ which guarantees that the limit

$$\begin{aligned} A&= \lim _{\tilde{x}\rightarrow 1}\frac{R(t,\tilde{x})}{(1-\tilde{x})^{\xi }}\\&= \lim _{\tilde{x}\rightarrow 1} \frac{1}{2\xi (1\!-\!\tilde{x})^{\xi -1}} \left( \frac{D(t)}{(1\!-\!\tilde{x})^{3/2}}\!-\!\frac{L(t)\tau ^{\prime }(x)}{(t\!-\!\tau (x))^{3/2}} \right) \end{aligned}$$

does not turn to zero or infinity. Due to this assumption, we can write

$$\begin{aligned} \frac{1}{\sqrt{t\!-\!\tau (x)}}\!=\! \frac{D(t)}{\sqrt{1\!-\!\tilde{x}}}\!+\!A (1\!-\!\tilde{x})^{\xi }\!+\!o\left( (1\!-\!\tilde{x})^{\xi }\right) , \end{aligned}$$

(77)

when $\tilde{x} \rightarrow 1$, or equivalently $x \rightarrow l(t)$. Taking the last estimate into account $A$ can be expressed as:

$$\begin{aligned} A&= \lim _{\tilde{x}\rightarrow 1}\frac{1}{2\xi (1\!-\!\tilde{x})^{\xi -1}} \left( \frac{D(t)}{(1\!-\!\tilde{x})^{3/2}}\!-\!\frac{L(t)\tau ^{\prime }(x)}{t\!-\!\tau (x)} \frac{D(t)}{\sqrt{1\!-\!\tilde{x}}} \right) \\&\quad -\frac{AL(t)}{2\xi }\lim _{\tilde{x}\rightarrow 1}\frac{(1-\tilde{x})\tau ^{\prime }(x)}{t-\tau (x)}\big (1+o(1)\big ). \end{aligned}$$

Now, on substitution of $\tau ^{\prime }(x)=1/L^{\prime }(t)$ at $x=L(t)$ and (75) into the limit one has:

$$\begin{aligned} A\!=\!\lim _{\tilde{x}\rightarrow 1}\frac{D(t)}{2\xi (1\!-\!\tilde{x})^{\xi -1/2}} \left( \frac{1}{1\!-\!\tilde{x}}\!-\!\frac{L(t)\tau ^{\prime }(x)}{t\!-\!\tau (x)} \right) \!-\!\frac{AL(t)D^2(t)}{2\xi L^{\prime }(t)}. \end{aligned}$$

Applying (75) and (22) here gives:

$$\begin{aligned}&\frac{1+2\xi }{2\xi }A\!=\!\lim _{\tilde{x}\rightarrow 1}\frac{D(t)}{2\xi (1\!-\!\tilde{x})^{\xi -1/2}} \left( \frac{1}{1\!-\!\tilde{x}}\!-\!\frac{L(t)\tau ^{\prime }(x)}{\sqrt{t\!-\!\tau (x)}} \frac{D(t)}{\sqrt{1\!-\!\tilde{x}}} \right) \\&\quad -\frac{AD(t)L(t)}{2\xi }\lim _{\tilde{x}\rightarrow 1}\frac{\tau ^{\prime }(x)\sqrt{1-\tilde{x}}}{\sqrt{t-\tau (x)}}. \end{aligned}$$

By repeating the same process one more time we have:

$$\begin{aligned} (2+2\xi )A\!=\!\lim _{\tilde{x}\rightarrow 1}\frac{D(t)}{(1\!-\!\tilde{x})^{\xi }} \left( \frac{1}{\sqrt{1\!-\!\tilde{x}}}\!-\!\frac{L(t)\tau ^{\prime }(x)D(t)}{\sqrt{t\!-\!\tau (x)}} \right) . \end{aligned}$$

Finally by eliminating the square root with use of (77) we obtain (after some algebra)

$$\begin{aligned} (3+2\xi )A=D(t)\lim _{\tilde{x}\rightarrow 1} \frac{1-L(t)\tau ^{\prime }(x)D^2(t)}{(1-\tilde{x})^{\xi +1/2}}. \end{aligned}$$

This relationship gives a finite value of $A$ if and only if $\xi =1/2$ and, as a result, we find:

$$\begin{aligned} A=\frac{1}{4}D^3(t)L^2(t)\tau ^{\prime \prime }(L(t))= -\frac{1}{4}\frac{L^{\prime \prime }(t)}{L^{\prime }(t)}\sqrt{\frac{L(t)}{L^{\prime }(t)}}. \end{aligned}$$

1.2 B Asymptotics of the solutions for different leak-off functions

Asymptotic expansion for the crack opening and the fluid velocity near the crack tip in the normalised variables (18) can be written in the following general forms:

$$\begin{aligned} w(t, x)\!=\!\sum _{j=0}^Nw_j(t)(1\!-\!x)^{\alpha _j}\!+\!O((1\!-\! x)^{\varrho _w}), \quad x\rightarrow 1,\nonumber \\ \end{aligned}$$

(78)

and

$$\begin{aligned} V(t, x)\!=\!\sum _{j=0}^NV_j(t)(1\!-\!x)^{\beta _j}\!+\!O((1\!-\! x)^{\varrho _V}),\quad x\rightarrow 1,\nonumber \\ \end{aligned}$$

(79)

with $\varrho _w>\alpha _n$, $\varrho _V>\beta _n$, $\alpha _0=1/3$, $\beta _0=0$ and some increasing sequences $\alpha _0,\alpha _1,\ldots ,\alpha _n$ and $\beta _0,\beta _1,\ldots ,\beta _n$. Note that the asymptotics are related to each other by the speed Eq. (19) and thus, regardless of the chosen leak-off function, we can write

$$\begin{aligned}&\sum _{j=0}^NV_j(t)(1-x)^{\beta _j}+\ldots =\frac{1}{3L(t)}\sum _{k=0}^N\sum _{m=0}^N\sum _{j=0}^N\nonumber \\&\quad \times \alpha _jw_j(t) w_m(t)w_k(t)(1-x)^{\alpha _j+\alpha _m+\alpha _k-1}.\nonumber \\ \end{aligned}$$

(80)

In line with the discussion after Eq. (16), we are interested only in the terms such that $\beta _j\le 1$, restricting ourselves to the smallest $\varrho _V>1$, since the values of $\beta _j$ are combinations of a sum of three consequent components of the exponents $\alpha _j$. However, since $\alpha _0$ is known $(\alpha _0=1/3)$, one can write (compare with (17)):

$$\begin{aligned}&V_{0}(t)=\frac{1}{3L(t)}w_{0}^{3}(t), \end{aligned}$$

(81)

$$\begin{aligned}&V_{1}(t)\!=\!\frac{1}{L(t)}\left( \alpha _1\!+\!\frac{2}{3}\right) w_{0}^{2}(t) w_{1}(t),\quad \beta _1\!=\!\alpha _1-\frac{1}{3}.\nonumber \\ \end{aligned}$$

(82)

To continue the process one now needs to compute the value of the exponent $\alpha _1$ as it is not clear a priori which value determining the next exponent $\beta _2=\min \{2/3+\alpha _2,1/3+2\alpha _1\}$ is larger. To do so let us rewrite the continuity Eq. (20) in the form:

$$\begin{aligned} \frac{\partial w}{\partial t}\!+\!\frac{V_0(t)}{L( t)}(1\!-\!x)\frac{\partial w}{\partial x}\!=\!\frac{1}{L( t)}\frac{\partial \big (w(V_0\!-\!V)\big )}{\partial x}\!-\!q_l(t,x).\nonumber \\ \end{aligned}$$

(83)

Here, the terms on the left-hand side of the equation are always bounded near the crack tip, while those on the right-hand side can behave differently depending on the chosen leak-off function.

Consider the following three cases of $q_l$ behaviour.

• (i) Assume first that

  $$\begin{aligned} q_l(t,x)=o\big (w(t,x)\big ),\quad x\rightarrow 1. \end{aligned}$$

  This case naturally includes the impermeable rock formation. Analysing the leading order terms in the Eq. (83), it is clear that $w(V_0-V)=O((1-x)^{4/3})$, as $x\rightarrow 1$. This, in turn, is only possible for $\beta _1=1$ and, therefore, $\alpha _1=4/3$. Finally, comparing the left-hand side and the right-hand side of the equation we obtain:

  $$\begin{aligned}&w_{0}^{\prime }(t)= \frac{w_{0}(t)}{3L(t)} \big (V_0(t)+4V_1(t)\big ),\nonumber \\&V_{1}(t)=\frac{2}{L(t)}w_{0}^{2}(t) w_{1}(t). \end{aligned}$$

  (84)

  This case has been considered in Linkov (2011d) and Mishuris et al. (2012).

• (ii) If we assume that the leak-off function is estimated by the solution as $O\big (w(t,x)\big )$, or equivalently;

  $$\begin{aligned} q_l(t,x)\sim \Upsilon (t)w_0(t)(1-x)^{1/3},\quad x\rightarrow 1, \end{aligned}$$

  then the previous results related to the values of $\alpha _1$ and $\beta _1$ and, therefore, the Eq. (84)$_2$ remain the same, while the first one changes to

  $$\begin{aligned} w_{0}^{\prime }(t)= \frac{1}{3L(t)} w_{0}(t)\big (V_0(t)\!+\!4V_1(t)\big )\!-\!\Upsilon (t)w_0(t). \nonumber \\ \end{aligned}$$

  (85)

  This case corresponds to (21)$_3$ when $C_{32}=0$ and $\Upsilon (t)=kC_{31}(t)$.

• (iii) The leak-off function in a general form:

  $$\begin{aligned} q_l(t,x)\!=\!\Phi (t)(1\!-\!x)^{\theta }\!+\!o((1\!-\!x)^{1/3}),\quad x\rightarrow 1, \end{aligned}$$

  where $-1/2\le \theta <1/3$. Here, one can conclude that $w(V_0-V)=O((1-x)^{1+\theta })$, as $x\rightarrow 1$ or equivalently, $\beta _1=\theta +2/3$, and $\alpha _1=1+\theta$. Moreover, in this case:

  $$\begin{aligned}&(1+\theta )w_0V_1=L(t)\Phi (t), \nonumber \\&\quad V_{1}(t)=\frac{1}{L(t)}\left( \theta +\frac{4}{3}\right) w_{0}^{2}(t) w_{1}(t), \end{aligned}$$

  (86)

  and, thus

  $$\begin{aligned} w_{1}(t)=\frac{3L^2(t)\Phi (t)}{(4+3\theta )(1+\theta )w_{0}^{3}(t)}. \end{aligned}$$

  (87)

  Note, that as one would expect, the particle velocity function is not smooth in this case near the crack tip, its derivative is unbounded and exhibits the following behaviour:

  $$\begin{aligned} \frac{\partial V}{\partial x}=O\big ((1-x)^{\theta -1/3}\big ), \quad x\rightarrow 1. \end{aligned}$$

  To formulate the equation similar to (84)$_1$ or (85), one needs to continue asymptotic analysis of the Eq. (83) incorporating the available information. Apart from the fact that the analysis can be done in the general case, we restrict ourselves only to three variants used from the beginning (compare (4)), respectively: $\theta =0,\, \theta =1/3-1/2=-1/6$ and $\theta =-1/2$.

When $\theta =0,\, \alpha _1=1$ and $\beta _1=2/3$, returning to the Eq. (80), one concludes that $\beta _2>1$ and, therefore,

$$\begin{aligned} w_{0}^{\prime }(t)= \frac{1}{3L(t)} w_0(t)V_0(t). \end{aligned}$$

(88)

This case corresponds to (21)$_3$ when $\Phi (t)=C_3^{(2)} (t)w_0(t)$ and $C_3^{(1)}=0$.

If $\theta =-1/6$, then $\alpha _1=5/6$ and $\beta _1=1/2$. In this case the function $\Phi (t)$ can be written as $\Phi (t)=C_2 D(t)w_0(t)$ [compare to (21)$_2$] and again Eq. (80) gives $\beta _2>1$, while Eq. (83) leads to

$$\begin{aligned} w_{0}^{\prime }(t)= \frac{1}{3L(t)} \big (w_0(t)V_0(t)+4w_1(t)V_1(t)\big ). \end{aligned}$$

(89)

Summarizing, in both mentioned above cases, there exists a single term in asymptotics of the particle velocity which has singular derivative near the crack tip. Moreover, those terms ($w_1$ and $V_1$, respectively) are fully defined by the leak-off function $\Phi (t)$ and the coefficient $w_0$ in front of the leading term for the crack opening in (87) and (86)$_1$.

The situation changes dramatically when $\theta =-1/2$ (Carter law). We now have $\alpha _1=1/2$ and $\beta _1=1/6$ and $\Phi (t)=C_1D(t)$. In this case, however, $\beta _2<1$ and we need to continue the asymptotic analysis further to evaluate all terms of the particle velocity which exhibit non-smooth behaviour near the crack tip. We omit the details of the derivation, presenting only the final result in a compact form. The first six exponents in the asymptotic expansions (78) and (79), that introduce the singularity of $w_x$, are:

$$\begin{aligned}&\alpha _j=\frac{1}{2}+\frac{j}{6},\quad \beta _j=\frac{j}{6},\quad j=1,2,\ldots ,6.\\&w_j(t)=\kappa _j\frac{\Phi ^j(t)L^{2j}(t)}{w_0^{4j-1}(t)},\quad V_j(t)=\psi _j\frac{\Phi ^j(t)L^{2j-1}(t)}{w_0^{4j-3}(t)}, \end{aligned}$$

where $j=1,2,\ldots ,5$ and

$$\begin{aligned} \begin{array}{l} \kappa _1=\frac{12}{7},\quad \psi _1=2,\quad \kappa _2=-\frac{270}{49},\quad \psi _2=-\frac{24}{7},\\ \kappa _3=\frac{9768}{343}, \,\, \psi _3=\frac{828}{49},\,\,\kappa _4=-\frac{2097252}{12005},\,\, \psi _4=-\frac{5136}{49},\\ \kappa _5=\frac{1081254096}{924385},\quad \psi _5=\frac{1234512}{1715}. \end{array} \end{aligned}$$

1.3 C Benchmark solutions

There are several benchmarks in the literature to be utilized for investigation of the numerical algorithms. Benchmark solutions for impermeable rock have been constructed in Kemp (1989); Linkov (2011d), while that corresponding to the non-zero leak-off model with $q_l$ vanishing at a crack tip has been analyzed in Mishuris et al. (2012).

In this paper, we introduce three different analytical benchmark solutions corresponding to the representations (21). Moreover, for each of the leak-off functions under consideration we take two different relationships between the injection flux rate $q_0$ and the leak-off to formation $q_l$. In this way six different benchmark solutions are analyzed.

In order to formulate the benchmark solutions let us assume the following form of the crack opening function:

$$\begin{aligned} w(t,x)=W_0(1+t)^\gamma h(x),\quad W_0=\root 3 \of {\frac{3}{2}(3\gamma +1)},\nonumber \\ \end{aligned}$$

(90)

where $\gamma$ is an arbitrary parameter, and the function $h(x)$ ($0<x<1$) is given by:

$$\begin{aligned} h(x)=(1- x)^\frac{1}{3}+b_1(1-x)^{\lambda _1}+b_2(1-x)^{\lambda _2}. \end{aligned}$$

(91)

The choice of the next powers $1/3<\lambda _1<\lambda _2$ will depend on the leak-off variant from (4). On consecutive substitutions of (90)–(91) into the relations (19), (24), (29) and (31) one obtains the remaining benchmark quantities:

$$\begin{aligned}&L(t)=(1+t)^{\frac{3\gamma +1}{2}},\nonumber \\&V(t,x)=-W_0^3(1+t)^{\frac{3\gamma -1}{2}}h^2(x)\frac{\partial h}{\partial x}.\end{aligned}$$

(92)

$$\begin{aligned}&q_0(t)=-W_0^4(1+t)^{\frac{5\gamma -1}{2}}\left( h^3\frac{\partial h}{\partial x}\right) |_{x=0}.\end{aligned}$$

(93)

$$\begin{aligned}&q_l(t,x)=W_0(1+t)^{\gamma -1}\times \nonumber \\&\Big (\frac{3}{2}(3\gamma \!+\!1)\Big [\frac{1}{3}x\frac{\partial h}{\partial x} \!+\!3h^2\left( \frac{\partial h}{\partial x}\right) ^2 \!+\!h^3\frac{\partial ^2 h}{\partial x^2}\Big ]\!-\!\gamma h\Big ).\nonumber \\ \end{aligned}$$

(94)

It can be easily checked that for $\lambda _1=1/2$ and $\lambda _2=4/3$ the leak-off function incorporates a square root singular term of type (21)$_1$. By setting $\lambda _1=5/6$ and $\lambda _2=4/3$ we comply with representation (21)$_2$. Although in both of these cases $q_{1(2)}^*$ exhibits a singular behaviour at the crack tip, it does not detract from the applicability of our benchmarks. Finally, when using $\lambda _1=4/3$ and $\lambda _2=7/3$, the benchmark gives a non-singular leak-off function in the form (21)$_3$.

Note also, that by manipulating with the value of $\gamma$ one can simulate some very specific regimes of crack propagation. For example $\gamma =1/5$ corresponds to the constant injection flux rate, while $\gamma =1/3$ gives a constant crack propagation speed. For our computations we always set the value of $\gamma =1/5$.

Choosing appropriate values $b_1$ and $b_2$ one can change the relation between the amount of fluid loss to formation and the injection rate. This ratio can be defined by the measure, $Q_l/q_0$, where $Q_l$ is the total volume of leak-off $\int _0^1 q_l dx$. It is important to note that this measure decreases in time, from its maximum value to zero, for all chosen benchmarks. Thus, taking the maximal value high enough and tracing the solution accuracy in time, one can analyse performance of the algorithm for any possible value of the parameter. We consider two variants of $Q_l/q_0$, one where fluid injection doubles the size of total fluid loss, and a second where the total fluid loss is close to injection rate. The values of the corresponding parameters $b_1,\, b_2$ are presented in Table 9.

Table 9 The values of parameters $b_1$ and $b_2$ for different benchmark solutions modelling desired leak-off to fluid injection ratios

Additionally one can compute a parameter $\gamma _v$ defined in Mishuris et al. (2012) as a measure of the uniformity of fluid velocity distribution:

$$\begin{aligned} \gamma _v=\left[ \max (V(t,x))-\min (V(t,x))\right] \left[ \int \limits _0^1V(t,\xi )d\xi \right] ^{-1}.\nonumber \\ \end{aligned}$$

(95)

Interestingly, this measure is directly correlated with the leak-off ratio $Q_l/q_0$.

In Fig. 14 the distributions of the leak-off functions and the corresponding particle velocities for the respective benchmarks are presented. It shows that the velocity near the crack tip depends strongly on the benchmark variant. To highlight this fact, a zoom picture is placed in the Fig. 14b.

Fig. 14

Distributions of the leak-off functions $q_l(t,x)$ and the respective particle velocity $V(t,x)$ over $\tilde{x} \in (0,1)$ at initial time $t=0$

Note that the benchmark $q_l^{(1)}$ is worse, in a sense, than the original Carter’s model as it contains additional singular terms of the leak-off function. These terms are absent in the normalised Carter’s law as it follows from “Appendix B”.