An approximate solution for a plane strain hydraulic fracture that accounts for fracture toughness, fluid viscosity, and leak-off

Abstract

The goal of this paper is to develop an approximate solution for a propagating plane strain hydraulic fracture, whose behavior is determined by a combined interplay of fluid viscosity, fracture toughness, and fluid leak-off. The approximation is constructed by assuming that the fracture behavior is primarily determined by the three-process (viscosity, toughness, and leak-off) multiscale tip asymptotics and the global fluid volume balance. First, the limiting regimes of propagation of the solution are considered, that can be reduced to an explicit form. Thereafter, applicability regions of the limiting solutions are investigated and transitions from one limiting solution to another are analyzed. To quantify the error of the constructed approximate solution, its predictions are compared to a reference numerical solution. Results indicate that the approximation is able to predict hydraulic fracture parameters for all limiting and transition regimes with an error of under one percent. Consequently, this development can be used to obtain a rapid solution for a plane strain hydraulic fracture with leak-off, which can be used for quick estimations of fracture geometry or as a reference solution to evaluate accuracy of more advanced hydraulic fracture simulators.

Appendices

Appendix 1: Functions $g_\delta ({\hat{K}},{\hat{C}})$ and $\varDelta ({\hat{K}},{\hat{C}})$

This appendix provides expressions for the functions $g_\delta \left( \hat{K},{\hat{C}}\right)$ and $\varDelta \left( {\hat{K}},{\hat{C}}\right)$ that are used in the paper to approximate the solution for a semi-infinite hydraulic fracture. Note that these functions were obtained in (Dontsov and Peirce 2015, 2017).

With the reference to the scaling (12) and the fact that $\dot{l}=\alpha l/t$, the function f can be introduced as

$$\begin{aligned} {\hat{s}}= & {} \dfrac{1}{3C_1(\delta ) }\left[ 1-{\hat{K}}^3-\frac{3}{2}{\hat{C}}{\hat{b}} \left( 1-{\hat{K}}^2\right) +3 {\hat{C}}^2{\hat{b}}^2\left( 1-{\hat{K}}\right) \right. \nonumber \\&\left. -\,3{\hat{C}}^3{\hat{b}}^3\ln \left( \dfrac{{\hat{C}}{\hat{b}}+1}{{\hat{C}}\hat{b}+{\hat{K}}}\right) \right] \equiv f\left( {\hat{K}},{\hat{C}}{\hat{b}},C_1\right) , \end{aligned}$$

(62)

where ${\hat{b}}={C_2(\delta )}/{C_1(\delta )}$ and

$$\begin{aligned} C_1(\delta )= & {} \dfrac{4(1-2\delta )}{\delta (1-\delta )} \tan \left( \pi \delta \right) ,\nonumber \\ C_2(\delta )= & {} \dfrac{16(1-3\delta )}{3\delta (2-3\delta )} \tan \left( \dfrac{3\pi }{2}\delta \right) . \end{aligned}$$

(63)

The zeroth-order approximation for the solution can be written as

$$\begin{aligned} {\hat{s}}=f\left( {\hat{K}},\dfrac{3\beta _{{\tilde{m}}}^4}{4\beta _m^3}\hat{C},\frac{\beta _m^3}{3}\right) \equiv g_0\left( {\hat{K}}, {\hat{C}}\right) , \end{aligned}$$

(64)

where $\beta _{{\tilde{m}}}={4}/{\left( 15^{1/4}\left( \sqrt{2}-1\right) ^{1/4}\right) }$ and $\beta _{m}=2^{1/3}3^{5/6}$. As mentioned in (Dontsov and Peirce 2015), the solution varies spatially as $w_a(s)\propto s^{{\bar{\delta }}}$, where ${\bar{\delta }} =\tfrac{1}{2}(1+\delta )$ and the power $\delta$ is given by

$$\begin{aligned} \delta =\dfrac{\beta _m^3 }{3} \left( 1+\dfrac{3\beta _{\tilde{m}}^4}{4\beta _m^3}{\hat{C}}\right) g_0\left( {\hat{K}}, {\hat{C}}\right) \equiv \varDelta \left( {\hat{K}},{\hat{C}}\right) ,\nonumber \\ \end{aligned}$$

(65)

which defines the function $\varDelta ({\hat{K}},{\hat{C}})$ and leads to the relation ${\bar{\delta }} =\tfrac{1}{2}\left( 1+\varDelta \left( {\hat{K}},{\hat{C}}\right) \right)$. By substituting (65) into (62), the $\delta$-corrected solution (12) can be written as

$$\begin{aligned} {\hat{s}}=f\left( {\hat{K}},{\hat{C}}{\hat{b}}\left( \varDelta ({\hat{K}},\hat{C})\right) ,C_1\left( \varDelta ({\hat{K}},{\hat{C}})\right) \right) \equiv g_\delta ({\hat{K}}, {\hat{C}}),\nonumber \\ \end{aligned}$$

(66)

which defines the function $g_\delta ({\hat{K}}, {\hat{C}})$.

Appendix 2: Numerical scheme

To construct the numerical scheme, Eq. (2) is rewritten using $\xi =x/l(t)$ and the scaling (47)–(49) as

$$\begin{aligned}&\dfrac{\partial \varOmega }{\partial \tau }- \dfrac{\xi V}{\gamma }\dfrac{\partial \varOmega }{\partial \xi } -\dfrac{1}{\gamma ^2}\dfrac{\partial }{\partial \xi }\left( \varOmega ^3\dfrac{\partial \varPi }{\partial \xi } \right) \nonumber \\&\quad +\dfrac{1}{\sqrt{\tau -\tau _0(\gamma \xi )}}=\dfrac{1}{\gamma }\delta (\xi ),\nonumber \\&\qquad \varOmega \rightarrow K_m\gamma ^{-1/2} (1-\xi )^{1/2},~~~ \xi \rightarrow 1, \end{aligned}$$

(67)

where $V={\dot{\gamma }}$, while the elasticity Eq. (3) is reduced to

$$\begin{aligned} \varPi (\xi ,\tau )= & {} -\dfrac{1}{2\pi \gamma }\int _0^1 M(\xi ,s) \dfrac{\partial \varOmega (s,\tau )}{\partial s}\,{\hbox {d}}s,\nonumber \\ M(\xi ,s)= & {} \dfrac{\xi }{\xi ^2-s^2}. \end{aligned}$$

(68)

The spatial coordinate $\xi$ is discretized as $\xi _j=(\tfrac{1}{2}+j) \varDelta \xi$, $j=1\ldots N$, in which case $\xi _1=\tfrac{1}{2}\varDelta \xi$ and $\xi _{N}=1-\tfrac{1}{2}\varDelta \xi$, and the temporal coordinate $\tau$ is discretized uniformly on a logarithmic scale. Piecewise constant approximation for $\varOmega$ is used, in which case $\varOmega ^i_j=\varOmega (\xi _j,\tau _i)$ and the vector $\varvec{\varOmega }^i$ represents an array of values of $\varOmega ^i_j$ for all j. In this situation, the elasticity Eq. (68) is discretized as

$$\begin{aligned} \varvec{\varPi }^i= & {} \varvec{C} \varvec{\varOmega }^i,\nonumber \\ C_{mn}= & {} \dfrac{1}{2\pi \gamma ^i} \left[ M\left( \xi _m,\xi _n+\tfrac{1}{2}\varDelta \xi \right) -M\left( \xi _m,\xi _n-\tfrac{1}{2}\varDelta \xi \right) \right] .\nonumber \\ \end{aligned}$$

(69)

Fluid balance (67), on the other hand, is discretized using backward time differencing

$$\begin{aligned} \dfrac{ \varvec{\varOmega }^i - \varvec{\varOmega }^{i-1}}{\varDelta \tau }=\varvec{B} \varvec{\varOmega }^i+\varvec{A}\left( \varvec{\varOmega }^i\right) \varvec{\varPi }^i+\varvec{S}^i, \end{aligned}$$

(70)

where the term that captures the moving mesh is discretized as

$$\begin{aligned} \left[ \varvec{B} \varvec{\varOmega }^i\right] _j= & {} \dfrac{\xi _j V}{\gamma }\dfrac{{\varOmega }_{j+1}^i-{\varOmega }^i_{j-1}}{2\varDelta \xi },\qquad j=2\ldots N-1,\\ \left[ \varvec{B} \varvec{\varOmega }^i\right] _1= & {} \dfrac{\xi _1 V}{\gamma }\dfrac{{\varOmega }_{2}^i-{\varOmega }^i_{1}}{2\varDelta \xi }, \\ \left[ \varvec{B} \varvec{\varOmega }^i\right] _N= & {} -\dfrac{\xi _N V}{\gamma }\dfrac{\varOmega ^i_N+{\varOmega }^i_{N-1}}{2\varDelta \xi }, \end{aligned}$$

the lubrication term is discretized as

$$\begin{aligned} \left[ \varvec{A}\left( \varvec{\varOmega }^i\right) \varvec{\varPi }^i\right] _j= & {} \dfrac{1}{\gamma ^2 \varDelta \xi }\left[ \left( \varOmega _{j+1/2}^i\right) ^3\dfrac{\varPi ^i_{j+1}-\varPi ^i_{j}}{\varDelta \xi }\right. \\&\left. - \left( \varOmega _{j-1/2}^i\right) ^3\dfrac{\varPi ^i_{j}-\varPi ^i_{j-1}}{\varDelta \xi } \right] ,~~ j=2\ldots N-1,\\ \left[ \varvec{A}\left( \varvec{\varOmega }^i\right) \varvec{\varPi }^i\right] _1= & {} \dfrac{\left( \varOmega _{3/2}^i\right) ^3}{\gamma ^2 \varDelta \xi }\dfrac{\varPi ^i_{2}-\varPi ^i_{1}}{\varDelta \xi } ,\\ \left[ \varvec{A}\left( \varvec{\varOmega }^i\right) \varvec{\varPi }^i\right] _N= & {} -\dfrac{ \left( \varOmega _{N-1/2}^i\right) ^3}{\gamma ^2 \varDelta \xi }\dfrac{\varPi ^i_{N}-\varPi ^i_{N-1}}{\varDelta \xi } , \end{aligned}$$

where the widths at the mid points $j\pm 1/2$ are calculated as an average between the corresponding values of the widths, and the source/leak-off term is

$$\begin{aligned} {S}_j^i= & {} -\dfrac{2}{\varDelta \tau } \left[ \sqrt{\tau _i+\varDelta \tau -\tau _0(\gamma \xi _j)} - \sqrt{\tau _i-\tau _0\left( \gamma \xi _j\right) } \right] \nonumber \\&+\dfrac{\delta _{1j}}{2\gamma \varDelta \xi }, \end{aligned}$$

where $j=1\ldots N-1$ and $\delta _{1j}$ is the Kronecker delta. The values of $\left[ \varvec{B} \varvec{\varOmega }^i\right] _N$ and $\left[ \varvec{A}( \varvec{\varOmega }^i) \varvec{\varPi }^i\right] _N$ are obtained by integration of (67) over the last element and using the no flux condition at the tip.

Similar to (Dontsov and Peirce 2017; Peirce 2016), in order to capture the multiscale behavior near the fracture tip, the following propagation condition is used

$$\begin{aligned} \varOmega ^i_{N-1}=\varOmega _a\left( \tfrac{3}{2}\gamma \varDelta \xi \right) , \end{aligned}$$

where $\varOmega _a$ is the scaled tip asymptotic solution. The latter equation implies that the numerical solution follows the asymptotic solution from the penultimate element to the tip, and allows one to determine the propagation velocity V. To successfully use this condition, pressure at the tip element $\varPi ^i_N$ is treated as an unknown. Since the tip asymptotic solution satisfies (12), it follows that

$$\begin{aligned} {\hat{s}}= & {} g_\delta \left( {\hat{K}},{\hat{C}}\right) ,\qquad {\hat{K}}=\dfrac{K_m d^{1/2}}{\varOmega ^i_{N-1}},\nonumber \\ {\hat{C}}= & {} \dfrac{2 d^{3/2} }{\left( \varOmega ^i_{N-1}\right) ^{5/2} \hat{s}^{1/2}},\qquad V=\dfrac{{\hat{s}} \,\left( \varOmega ^i_{N-1}\right) ^3}{d^2}, \end{aligned}$$

(71)

where $d=\tfrac{3}{2}\gamma \varDelta \xi$ signifies the distance from the center of the penultimate element to the tip. The above equation is solved for ${\hat{s}}$ using Newton’s method, and the velocity of propagation is calculated. Since the solution in the tip element follows the asymptotic solution, it is possible to determine its average opening as

$$\begin{aligned} \varOmega ^i_{N}=\left( \dfrac{2}{3}\right) ^{(3+\delta )/2}\, \dfrac{2\, \varOmega ^i_{N-1} }{3+\delta } \dfrac{d}{\varDelta \xi }, \end{aligned}$$

(72)

where it is used that $\varOmega _s\propto d^{(1+\delta )/2}$ and $\delta =\varDelta ({\hat{K}},{\hat{C}})$. The factor 2 / 3 ensures that the calculation is performed only within the tip element (and does not continue to the middle of the penultimate element). Leak-off in the tip element is calculated by taking $\tau -\tau _0=\gamma (1-\xi )/V$, in which case

$$\begin{aligned} S^i_{N}=-2 \sqrt{\dfrac{V}{\gamma \varDelta \xi }}. \end{aligned}$$

(73)

The numerical scheme consists of solving (69)–(73) for $\varOmega ^i_j$ (for $j=1\ldots N-1$) and $\varPi ^i_N$, which is done iteratively for each time step. The fracture length is then updated as $\gamma ^i=\gamma ^{i-1}+V\varDelta \tau$.