Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A non-singular integral equation formulation of permeable semi-infinite hydraulic fractures driven by… Gomez, Daniel 2016

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2016_november_gomez_daniel.pdf [ 5.32MB ]
JSON: 24-1.0308780.json
JSON-LD: 24-1.0308780-ld.json
RDF/XML (Pretty): 24-1.0308780-rdf.xml
RDF/JSON: 24-1.0308780-rdf.json
Turtle: 24-1.0308780-turtle.txt
N-Triples: 24-1.0308780-rdf-ntriples.txt
Original Record: 24-1.0308780-source.json
Full Text

Full Text

A Non-Singular Integral Equation Formulation ofPermeable Semi-Infinite Hydraulic FracturesDriven by Shear-Thinning FluidsbyDaniel GomezB.Eng. and B.Sc. (Hons.), University of Saskatchewan, 2014A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Mathematics)The University of British Columbia(Vancouver)August 2016c© Daniel Gomez, 2016AbstractThis thesis considers the problem of semi-infinite hydraulic fractures driven byshear-thinning power-law fluids through a permeable elastic medium. In the recent workof Dontsov and Peirce [6], the authors reformulated the governing equations in a way thatavoids singular integrals for the case of a Newtonian fluid. Moreover, the authors con-structed an approximating ODE whose solutions accurately describe the fracture openingat little to no computational cost. The present thesis aims to extend their work to themore general case where the fracture is driven by a shear-thinning power-law fluid.In the first two chapters of this thesis we outline the relevant physical modelling,and discuss the asymptotic propagation regimes typically encountered in hydraulic frac-turing problems. This is followed by Chapters 4 and 5 where we reformulate the gov-erning equations as a non-singular integral equation, and then proceed to construct anapproximating ODE. In the final chapter we construct a numerical scheme for solving thenon-singular integral equation. Solutions obtained in this way are then used to gauge theaccuracy of solutions obtained by solving the approximating ODE.The most important results of this thesis center on the accuracy of using theapproximating ODE. In the final chapter we find that when the fluid’s power-law index isin the range 0.4 ≤ n ≤ 1, an appropriate method of solving the approximating ODE yieldssolutions whose relative errors are less than 1%. However, this relative error increaseswith decreasing values of n so that in the range 0 ≤ n < 0.4 it reaches a maximumvalue of approximately 6%. Thus, at least for values of 0.4 ≤ n ≤ 1 the approximatingODE presents an accurate and computationally fast alternative to solving the semi-infiniteproblem. The same can’t be said for values of 0 ≤ n < 0.4, but the methods presented inthis thesis may be used as a starting point for future work in this direction.iiPrefaceThis dissertation is original, unpublished, independent work by the author, D. Gomez.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Elasticity Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.1.1 Tip Asymptotics and Propagation Criterion . . . . . . . . . . . . . . 72.1.2 Dislocation Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.2 Fluid Mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.2.1 Poiseuille’s Law for Power-Law Fluids . . . . . . . . . . . . . . . . . 142.2.2 Carter’s Leak-Off Model . . . . . . . . . . . . . . . . . . . . . . . . . 152.2.3 Mass Conservation and The Lubrication Equation . . . . . . . . . . 162.3 Summary of Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . 17iv3 Problem Scaling and Asymptotic Regimes . . . . . . . . . . . . . . . . . . 193.1 Viscous Regime (M) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.2 Toughness Regime (K) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.3 Leak-Off Regime (C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.4 Comments on Negative Pressure . . . . . . . . . . . . . . . . . . . . . . . . 243.5 Comment on the n = 0 Case . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.6 Spatial Ordering of Parameter Regimes . . . . . . . . . . . . . . . . . . . . 254 Non-Singular Integral Equation Formulation . . . . . . . . . . . . . . . . 264.1 Recovering Asymptotic Regimes from the Non-Singular Integral Equation . 295 Reduction to an Approximating First-Order ODE . . . . . . . . . . . . . 315.1 The Constant δ˜ Approximation . . . . . . . . . . . . . . . . . . . . . . . . . 355.2 Solutions of the Constant δ˜ ODE . . . . . . . . . . . . . . . . . . . . . . . . 405.3 The δ˜-Corrected Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 416 Numerical Solution of Non-linear Integral Equation . . . . . . . . . . . . 436.1 Performance of Non-Singular Integral Equation Solver . . . . . . . . . . . . 476.2 Qualitative Properties of the Numerical Solution . . . . . . . . . . . . . . . 496.3 Comparison to Approximating ODE . . . . . . . . . . . . . . . . . . . . . . 537 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65Appendix A Closed Formulae for two Integrals . . . . . . . . . . . . . . . . . 68A.1 Proof of Closed Form Formula for∫∞0G(t)tα dt . . . . . . . . . . . . . . . . . . 68A.2 Proof of Closed Form Formula for∫∞0G′(t)tκ dt . . . . . . . . . . . . . . . . . 70vAppendix B Matlab Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72B.1 Newton’s Method Solver for Non-Singular Integral Equation . . . . . . . . . 72B.2 Zero Delta-Corrected Solution . . . . . . . . . . . . . . . . . . . . . . . . . . 75B.3 First Delta-Corrected Solution . . . . . . . . . . . . . . . . . . . . . . . . . 76viList of Figures2.1 Geometry of a semi-infinite hydraulic fracture. . . . . . . . . . . . . . . . . 52.2 Infinite notch modelling the fracture tip. . . . . . . . . . . . . . . . . . . . . 72.3 Plot of λ versus opening half-angle ϕ. . . . . . . . . . . . . . . . . . . . . . 92.4 Schematic of propagating hydraulic fracture. . . . . . . . . . . . . . . . . . 112.5 Edge dislocation schematic. . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.1 Fracture width power α in Ω ∼ βξα for each parametric regime: αk tough-ness, αcm leak-off, and αm viscous regime. . . . . . . . . . . . . . . . . . . . 243.2 Schematic showing ordering of regions where each parameter regime domi-nates (K: toughness, C: leak-off, M: viscosity). . . . . . . . . . . . . . . . . 255.1 Plots of kernel function G(t) and its first derivative G′(t). . . . . . . . . . . 325.2 Variation in C0 as δ˜ varies between 0 ≤ δ˜ ≤ 1. . . . . . . . . . . . . . . . . . 365.3 Variation in coefficients C1 and C2 as δ˜ varies between 0 ≤ δ˜ ≤ 2−n2+n forseveral 0 < n ≤ 1 values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37vii6.1 Convergence of partial sums making up coefficient of X2−n2+n2N+1 in (6.10). Theleft pane shows values of the partial sums (solid lines) and infinite sumsapproximated with M = 104 (dashed lines) for values of n = 0 at the topand n = 1 at the bottom in 0.10 steps. The right pane shows the absoluteerror of partial sums (note that this plot is actually showing a small bandthat includes all the n values plotted in the left pane). . . . . . . . . . . . . 466.2 Maximum total computation time (a) and ODE computation time (b) mea-sured in seconds, where the maximum is taken over 10−5 ≤ χ ≤ 105. . . . . 486.3 Maximum number of iterations (a) and final residual (b) for Newton’smethod, where the maximum is taken over 10−5 ≤ χ ≤ 105. . . . . . . . . . 486.4 Numerical results for n = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 506.5 Numerical results for n = 0.60. . . . . . . . . . . . . . . . . . . . . . . . . . 516.6 Numerical results for n = 0.20. . . . . . . . . . . . . . . . . . . . . . . . . . 516.7 Plots of solutions w˜(x˜) for various n values (see legend) and χ = 10−5, 103, 105. 526.8 Relative error between full numerical solution and solutions to the approx-imating ODE for n = 0. For ode45 and zero δ˜-corrected solution themaximum relative error is ≈ 0.062, while for the first δ˜-corrected solutionit is 0.048. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 546.9 Relative difference of approximating ODE solutions with regime solutions(toughness and viscosity). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 556.10 Relative error between full numerical solution and solutions to the approx-imating ODE for n = 0.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 566.11 Relative error between full numerical solution and solutions to the approx-imating ODE for n = 0.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58viii6.12 Relative error between full numerical solution and solutions to the approx-imating ODE for n = 0.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 586.13 Relative error between full numerical solution and solutions to the approx-imating ODE for n = 0.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 596.14 Relative error between full numerical solution and solutions to the approx-imating ODE for n = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 596.15 Relative error of zero- and first-delta corrected solution for χ = 105 withdashed vertical lines showing regime boundaries. . . . . . . . . . . . . . . . 606.16 Relative error of zero- and first-delta corrected solution for χ = 100 withdashed vertical lines showing regime boundaries. . . . . . . . . . . . . . . . 606.17 Relative error of zero- and first-delta corrected solution for χ = 10−5 withdashed vertical lines showing regime boundaries. . . . . . . . . . . . . . . . 61A.1 Schematic of the integration contour. . . . . . . . . . . . . . . . . . . . . . . 69ixAcknowledgementsI want to thank my supervisors Dr. Anthony Peirce and Dr. Michael Ward for theirguidance and support both inside and outside of the classroom. I would also like to thankDr. Egor Dontsov, whose work with Dr. Peirce at UBC motivated much of the work inthis thesis. Finally, I am grateful to NSERC for funding my first year in this programthrough a one-year NSERC CGS-M fellowship.xChapter 1IntroductionSince its inception and first implementation in the early 19th century [1], the methodof hydraulic fracturing as a means of stimulating oil reserves has become ubiquitous inthe oil and gas industry. Hydraulic fractures are best described as brittle fractures thatpropagate due to the injection of a highly viscous fluid. Outside of the oil and gas industryhydraulic fractures arise naturally as a central mechanism in the formation of kilometrelong dikes for magma transport, or they are artificially induced for applications in minecaving operations or soil remediation projects (see [4] and the references therein). In eachcase the composition of propagating fluid and scale of the fractures may differ, but thephenomenon of hydraulic fracturing remains the same. The quantitative modelling of hy-draulic fractures dates back to the late 1950s [1] and has been motivated both by naturallyoccurring phenomena, as well as by the prospect of better understanding and optimizingthe use of hydraulic fracturing in industrial applications. The study of hydraulic fractur-ing as it pertains to the oil and gas industry has been of particular interest due to its wideuse and its engineering, economic, and environmental implications. For this reason thisthesis will focus primarily on the study of hydraulic fractures in the context of reservoirstimulation.1To motivate the need for quantitative models of hydraulic fractures it is instructiveto first understand their role in the oil and gas industry. While oil may be extractedby conventional methods from a variety of reservoirs (e.g. sandstone formations), otherlow permeability rock formations, such as as shale, are not susceptible to conventionalmethods [5]. However, reservoir stimulation techniques have been developed which allowaccess to oil locked within these otherwise inaccessible reservoirs. In general reservoirstimulation projects share the following key steps. First a hole, known as a wellbore, isdrilled into the ground to access deep underground oil reservoirs. After the wellbore hasbeen established numerous explosive charges on the wellbore walls are discharged. Thiscreates weak regions in the surrounding rock that serve as the starting location of thehydraulic fractures. Next, a highly viscous fluid is pumped into the wellbore at a highpressure, which results in the propagation of hydraulic fractures at each of the dischargesites. This is followed by pumping in a polymer mixture responsible for the “filter cake”building mechanism that mediates leak-off of the fracturing fluid into the surroundingrock. After the hydraulic fractures have reached a sufficient width, a proppant mixtureis pumped into the fractures with the purpose of keeping them open while the fracturingfluid is removed. With the fracturing fluid removed, a negative pressure gradient drivesthe surrounding oil into the hydraulic fractures and up the wellbore for collection [1, 5].With the above procedure in mind we can motivate the study of quantitative mod-els of hydraulic fractures by considering two problems. The first is known as prematureproppant bridging and it occurs when the proppant is pumped while the hydraulic frac-ture is not sufficiently wide. This creates a barrier that does not allow the proppant toadequately fill the hydraulic fracture thereby causing issues when the fracturing fluid isremoved. Current technologies are unable to make fracture width measurements during ahydraulic fracturing treatment so having quantitative models that describe how the frac-2ture width evolves with time are crucial to prevent premature bridging. Another problemwe can consider is that of fracture efficiency which is typically defined as the ratio of thevolume of fluid stored in the fracture to the total volume pumped in. Ideally this ratiowould equal one but due to fluid loss through leak-off it is typically smaller. A quanti-tative understanding of how the problem parameters affect this ratio is clearly of greatimportance both in terms of optimizing resource usage as well as reducing potentiallynegative environmental effects.The highly complex nature of hydraulic fractures means that, with the exception ofa few simplified scenarios, solutions are obtained almost exclusively by numerical methods.However, even to obtain a numerical solution there are still significant hurdles that mustbe overcome. The most notable challenges include moving boundaries as well as themulti-scale behaviour inherent in hydraulic fracture propagation. Fortunately, substantialprogress can be made in addressing both these challenges by furthering our understandingof the hydraulic fracture’s tip behaviour. Moreover it can be shown (see Appendix B of[19]) that near its tip, a hydraulic fracture behaves like a semi-infinite hydraulic fractureunder plane strain conditions. In contrast to the general hydraulic fracturing problem, thesemi-infinite case is more amenable to analytical techniques and has therefore been thefocus of numerous studies. These studies typically involve finding asymptotic solutionsin distinguished parameter regimes and using this to then calculate the full semi-infinitesolution structure numerically as is done in [9, 10] for Newtonian fracturing fluids and[3, 15, 16] for power-law fluids.One of the key difficulties in the numerical solution of the semi-infinite problem isthe coupling between a singular integral equation and a non-linear ordinary differentialequation. The recent work by Dontsov and Peirce [6] overcomes this difficulty by refor-mulating this problem in terms of a single non-singular and non-linear integral equation,3which is readily solved using standard numerical methods such as Simpson’s rule andNewton’s method. Moreover Dontsov and Peirce obtained an approximating separableordinary differential equation from which they were able to obtain an approximating so-lution, which is both fast to calculate and is in excellent agreement with full numericalsolutions of the semi-infinite problem. The most important feature of this approximationis the speed with which its solutions can be computed. Indeed, multiscale hydraulic frac-ture simulators require the rapid computation of the tip asymptotes in order to run ina reasonable time. Motivated by these results as well as those obtained by Adachi andDetournay [3] as well as Linkov [15, 16], the present thesis considers the non-singular refor-mulation of the problem governing semi-infinite hydraulic fractures driven by a power-lawfluid. Furthermore, as was done by Dontsov and Peirce for Newtonian fluids, one of ourprimary goals is to obtain an approximate solution that is both fast to calculate andreliably captures the solution of the semi-infinite fracture. Future sections will demon-strate that we have partially succeeded in this goal by obtaining a solution that is fast tocalculate but has a limited range of reliability.The proceeding sections are organized in the following way. In Chapter 2 wewill discuss the relevant mathematical modelling, making more precise the underlyingassumptions of our simplified model. Chapter 3 discusses in detail the relevant scalingand resulting asymptotic regimes of propagation for hydraulic fractures driven by power-law fluids. In Chapter 4 we will extend the formulation of Dontsov and Peirce [6] to thepower-law fluid case, and this is followed by a discussion of the approximating ODE andits solution structure in Chapter 5. Finally, before the concluding remarks of Chapter7, in Chapter 6 we will discuss the numerical implementation of the reformulation fromChapter 4.4Chapter 2Mathematical ModelAs discussed in the previous chapter, our study of the semi-infinite hydraulic fractureproblem is largely motivated by its connection and practical relevance to more generalfracturing problems (see for example [4] and Appendix B of [19]). We therefore considerthe case of a hydraulic fracture under plane strain conditions driven by a power-law fluid.The fracture is assumed to propagate in the positive X direction at a constant velocity V(see Figure 2.1) so that Xtip = V t. The distance from the fracture tip as measured fromwithin the fracture is denoted by x and is given by x = Xtip−X. Furthermore we assumethat the fracture velocity is sufficiently slow so that elastic deformations depend solely onXYx XtipVwpσ0σ0Figure 2.1: Geometry of a semi-infinite hydraulic fracture.5the distance from the fracture tip, x, and are otherwise independent of time. We makereference to the fixed coordinate X only when describing the behaviour of the fracturingfluid, which in the absence of fluid lag is also assumed to travel in the positive X directionat the same constant velocity V .The modelling of hydraulic fractures involves the coupling of two sub-fields of con-tinuum mechanics – elasticity theory and fluid mechanics. We are primarily interestedin obtaining the fracture width w(x) (or w(X, t) in the fixed coordinates) from the mod-elling, but as we will soon see the fluid pressure pf (x) can be obtained as a by-product.To facilitate the modelling of hydraulic fractures, we can isolate four distinct but stronglycoupled physical phenomena. From the point of view of elasticity this involves the sin-gular behaviour of the stress field at the tip and the related criterion for propagation, aswell as the coupling between the internal pressure of the fracture and the resulting elasticdeformation. These two phenomena are discussed in the Elasticity Theory section belowand involve the use of the stress function approach as well as ideas from dislocation theory.The remaining two phenomena involve the behaviour of the fluid, the first modelling thepropagation of the fluid using lubrication theory, and the second modelling the loss offluid by leak-off using Carter’s leak-off model.2.1 Elasticity TheoryWe begin our discussion by exploring how the surrounding elastic material respondsto hydraulic fractures. This interaction occurs via two mechanisms – the geometric dis-turbance caused by the fracture, and by interactions with the fracturing fluid. Both ofthese interactions are aptly described using ideas from elasticity theory. In particular thegeometric effect is found to be concentrated at the fracture tip and can therefore be de-scribed using asymptotic approximations. Furthermore, as we will see the tip behaviour6is closely linked with the propagation criterion for the fracture by what is known as thestress intensity factor and the material toughness. Finally, using ideas from dislocationtheory the interaction between the elastic material and the fracturing fluid leads to arelationship between the fluid pressure and the fracture width.2.1.1 Tip Asymptotics and Propagation CriterionTo describe the tip asymptotics we follow the method of Williams [20] as set out in theexcellent textbook on dislocations by Hills et. al. [12] by first considering the problem ofa notch of angle 2ϕ within an infinite elastic material under plane strain conditions andwith symmetric uni-axial remote loading (see Figure 2.2). Furthermore we assume thatthe notch is stationary (i.e. ϕ is fixed) and therefore satisfies the steady Navier equations.Such a notch asymptotically describes the opening of a fracture near its tip. We seek anAiry stress function U defined byσθθ =∂2U∂r2, σrθ = − ∂∂r(1r∂U∂θ), σrr =1r2∂2U∂θ2+1r∂U∂r. (2.1)Then the Airy stress function satisfies the biharmonic equation∆2U ≡(∂2∂r2+1r∂∂r+1r2∂2∂θ2)2U = 0,for which we seek a solution of the form U = rλf(θ). Note that with U of this formσ = O(rλ−2) and the displacement field is u = O(rλ−1). We are interested in the dominantφφθrσ0σ0Figure 2.2: Infinite notch modelling the fracture tip.7behaviour for small r and therefore seek the smallest value of λ which, in order to ensureu remains finite, satisfies λ > 1. By defining g(θ) = f ′′(θ) + λ2f(θ) we findg′′(θ) + (λ− 2)2g(θ) = 0 =⇒ g(θ) = a cos(λ− 2)θ + b sin(λ− 2)θ,where a and b are arbitrary constants. By symmetry of the loading we expect U to besymmetric and therefore set b = 0. Thusf ′′(θ) + λ2f(θ) + a cos(λ− 2)θ =⇒ f(θ) = A cosλθ +B cos(λ− 2)θ,where B = a4(λ−1) is an arbitrary constant and we have again used symmetry to removethe sinusoidal solution sinλθ. Now, since the notch is stationary we impose the tractionfree boundary conditions σrθ = σθθ = 0 at θ = ±ϕ (by symmetry we only need to considerone angle, say θ = ϕ) which yields the homogeneous system of equationsAλ(λ− 1) sinλϕ+B(λ− 1)(λ− 2) sin(λ− 2)ϕ = 0,Aλ(λ− 1) cosλϕ+Bλ(λ− 1) cos(λ− 2)ϕ = 0.(2.2)Assuming λ 6= 0, 1 we deduce that for this system to have a non trivial solution it’sdeterminant must vanish, which yieldssin 2(λ− 1)ϕ+ (λ− 1) sin 2ϕ = 0.The solution λ = λ(ϕ) is shown in Figure 2.3 where we have chosen the range pi2 < ϕ < piwith the intention of describing the fracture opening by the wedge. This figure revealsthat the smallest value taken by λ > 1 is λ = 3/2 which occurs at ϕ = pi. This results inσ = O(r−1/2) which signifies a singularity in the stress field as the apex of the notch isapproached. Note that the next solution at ϕ = pi is λ = 5/2 which gives the correctionterm for U and consequently for the stress and displacement fields.8100 120 140 160 180' (degrees) 2.3: Plot of λ versus opening half-angle ϕ.With ϕ = pi and λ = 3/2 we find that U = r3/2(A cos 32θ + B cos12θ) and using(2.1) we deduceσrθ =KI√2pir(14sinθ2+14sin3θ2), (2.3)σrr =KI√2pir(54cosθ2− 14cos3θ2), (2.4)σθθ =KI√2pir(34cosθ2+14cos3θ2), (2.5)where KI = 3√2piA is known as the mode-I stress intensity factor. The square rootsingularity present in the stress field is ubiquitous in two-dimensional fracture problems.In particular the same singularity arises for both mode II and III fractures, which resultin the definition of stress intensity factors KII and KIII respectively (see Hills et. al. [12]for a detailed description).We now turn to the question of the propagation criterion. The stress field formulae(2.3)-(2.5) reveal that regardless of the fracture geometry, there will be a square rootsingularity at the tip. However, not all fractures propagate, which suggests that it is the9strength of the singularity, characterized by KI , that determines whether a fracture willpropagate. Indeed Griffith proposed that a fracture propagates only when the energyreleased by the fracture exceeds that required to create two new surfaces [11]. Thiscriterion takes the form [12]G ≥ 2γ, (2.6)where G is the energy release rate per unit area of the crack whereas 2γ represents theenergy required to create two new surfaces. Furthermore the energy release rate is relatedto the stress intensity factor via [12]G =4(1− ν)8µK2I ,whence we deduce that a fracture propagates only when the stress intensity factor exceedssome value KIc known as the material toughness.Assuming the fracture is slowly propagating, we replace KI with KIC in (2.3)-(2.5)whence a lengthy calculation yields the displacement fieldux(r, θ) =KIc2µ√r2picosθ2(3− 4ν − cos θ), uy(r, θ) =KIc2pi√r2pisinθ2(3− 4ν − cos θ),(2.7)where µ is the shear modulus and ν is Poisson’s ratio. The fracture opening width w(x) isrelated to the displacement uy(r, θ) by w(x) = uy(x, pi)−uy(x,−pi) where x is the distancefrom the tip. Thus, both the propagation criterion and the near field behaviour of thefracture opening can be encoded in the single equationw(x) ∼ 8√2piKIcE′x1/2 +O(x3/2) (x→ 0+),(2.8)where E′ = E1−ν2 and we have used that in plane strain the shear stress and Young’smodulus are related by E = 2(1 + ν)µ.10σ0σ0pf(x,t)Vw(x,t) xFigure 2.4: Schematic of propagating hydraulic fracture.2.1.2 Dislocation TheoryWe focus now on the interaction between the elastic material and the fracturing fluid.As depicted in Figure 2.4 we expect the fracturing fluid pressure pf to exert outward forceson the fracture walls, driving the fracture ever wider and therefore increasing the widthw. In this section we will be seeking a quantitative formulation of this interaction.The first step is to determine the effect of a fracture on the elastic material’s stressfield. This is done by first considering the effect on the stress field due to a single edgedislocation. We will keep our discussion of dislocations brief and limited in generality,taking only what we need for subsequent sections. Some excellent resources on the theoryof dislocations include the Book by Hills, et. al. [12] as well as the collection of papers editedby Mura [17] and in particular the papers by Hirth [13], Mura [18], and Dundurs [7] foundtherein. An edge dislocation represents a defect within an elastic material that results in ajump discontinuity b in the displacement field u as a clockwise path is followed about thedislocation core (see Figure 2.5a). We remark that the quantity b is known as the Burger’svector. We consider an edge dislocation centred at the origin with Burger’s vector (0, by)and such that there is a jump discontinuity in the y-component of the displacement fieldas the positive real axis is crossed (see Figure 2.5b). With the “+” and “−” symbols in11bdislocationcore+ ++ +-- --discontinuitycut(a) Arbitrary edge dislocation with burgers - - - -xy(b) Edge dislocation with cut along positive xaxis and burgers vector b = (0, by).Figure 2.5: Edge dislocation schematic.Figure 2.5b imposing an orientation, the jump in the displacement is uy(+)− uy(−) = by(note that the Burgers vector points downwards). Seeking a stress function that gives thisjump discontinuity in the displacement field we obtain the stress field [12]σyy(x) =E′4pibyx. (2.9)Having found the response to a single dislocation we now consider a distributionof infinitesimal dislocations, each centred at ξ ≥ 0 and having burgers vector δby(ξ).Introducing the dislocation density By(ξ) by defining δby(ξ) = By(ξ) dξ we now use (2.9)to obtain the stress field’s response to each infinitesimal dislocationδσyy(x; ξ) =E′4piBy(ξ)x− ξ dξ.Integrating over all dislocations we obtain the total responseσyy(x) =E′4pi−∫ ∞0By(ξ)x− ξ dξ, (2.10)where −∫is the Cauchy Principal value.It remains to determine the dependence of (2.10) on w. This is done by observingthat every infinitesimal dislocation centred at ξ < x contributes an amount −By(ξ)dξ to12the total jump discontinuity in the y-displacement at x and hencew(x) = −∫ x0By(ξ)dξ ⇐⇒ By(x) = −dwdx,whence (2.10) becomesσyy = −E′4pi−∫ ∞0∂w∂ξ(ξ, t)dξx− ξ .Finally we relate the stress field to the pressure p(x, t) = pf (x, t) − σ0 by noting thatthe stress field on the boundary σyy calculated above is related to p(x, t) via p(x, t) =−σyy(x, t) and hencep(x, t) =E′4pi−∫ ∞0∂w∂ξ(ξ, t)dsx− s ds (2.11)2.2 Fluid MechanicsHaving established the elastic material’s response to the hydraulic fracture, we nowexplore the behaviour of the propagating fluid within the hydraulic fracture. This isaccomplished in three parts. First we determine the appropriate form of Poiseuille’s lawfor a power-law fluid, which relates the pressure gradient to the fluid flux. Secondly wediscuss Carter’s leak-off model, which describes how some of the fluid escapes the hydraulicfracture and enters the elastic material. Finally the resulting equations from these twodiscussions are combined by considering mass-conservation which leads to the governingequation for the propagating fluid. Throughout this section we will begin by consideringthe fracturing fluid in a fixed reference frame with the X-Y axes oriented in such a waythat the fracture propagates along the positive X direction and is centred on the Y axis(see Figure 2.1). After establishing the relevant equations we reintroduce the distancefrom the tip x by changing to the moving coordinates x = V T −X.132.2.1 Poiseuille’s Law for Power-Law FluidsOur leading assumptions for the hydraulic fracture and the fluid it confines are that thefluid moves sufficiently slowly that the inertial forces are negligible in comparison to theviscous forces, the confining fracture is much longer than it is wide, the fluid’s velocity isrestricted to the X − Y plane, and the fluid satisfies a no-slip boundary condition alongthe fracture walls (i.e. along Y = ±w/2). In light of these assumptions, the formulation oflubrication theory applies. The derivation leading to Poiseuille’s law is standard, and forthis reason we omit many of the details, pointing out only those areas where the power-lawrheology of the fluid plays a critical rule. We recall that the power-law rheology takes theformτij = M(2ε˙ij)n = M(∂Vi∂Xj+∂Vj∂Xi)n, (2.12)where n is the power-law index, M the consistency index, and Vi the components of thevelocity field. We restrict the power-law index to 0 < n ≤ 1, which reflects the shear-thinning rheology of fluids typically used in reservoir stimulation projects [8]. In additionto assuming the fluid is laminar we also assume that the confining fracture is much longerthan it is wide. The fluid therefore satisfies the conditions of lubrication theory, whichleads to the following relationship between the pressure gradient and the X-component ofthe velocity∂p∂X= nM(∂VX∂Y)n−1∂2VX∂Y 2, in − w2< Y <w2; v = 0, on Y = ±w2.Using symmetry to eliminate the constant resulting from the first integration, and thenintegrating again yields the velocity profileVX(X,Y, t) = − nn+ 1(1M∣∣∣∣ ∂p∂X∣∣∣∣)1/n(|Y |n+1n − (w2)n+1n).14Integrating the velocity profile from Y = −w/2 to +w/2 yields the fluid flow rateq(X, t) ≡∫ w/2−w/2V (X,Y, t) dy =(1M ′∣∣∣∣ ∂p∂X∣∣∣∣)1/nw 2n+1n ,whereM ′ ≡ 2n+1(2n+ 1n)nM. (2.13)Finally, we use the fact that ∂p/∂X < 0 and q > 0 to solve for the pressure gradient,thereby obtaining Poseuille’s law for a power-law fluid∂p∂X= −M′qnw2n+1. (2.14)2.2.2 Carter’s Leak-Off ModelThe remaining phenomenon to be described is that of fluid leak-off. This phenomenonis described by Carter’s leak-off model as follows. First, as the fluid leaks some of it isdeposited in a region immediately surrounding the fracture making up what is knownas the filter cake. The fluid that makes it past the filter cake is known as the filtrateand the region it propagates in is known as the invaded zone. Beyond the invaded zonewe assume that the elastic material is naturally saturated with a reservoir fluid, therebyputting up a pressure barrier which resists the filtrate fluid. With this process in mind, theleading assumptions of Carter’s leak-off are that (1) the leak-off fluid flows perpendicularto the fracture boundary thereby rendering it a one-dimensional process, (2) the reservoirpressure pr is much larger than the fluid pressure pf , (3) the rate at which the filter cakegrows is directly proportional to the the amount of fluid passing through it, and (4) theleaking fluid satisfies Darcy’s law in both the invaded zone and the filter cake. In additionto these four assumptions, we add a fifth which addresses the power-law behaviour of thefracturing fluid and states that (5) the added polymers that gives the fracturing fluid it’spower-law rheology are used up to build the filter cake and therefore do not otherwise15contribute to the leak-off process. This final assumption effectively allows us to treatthe leaking fluid as a Newtonian fluid which greatly simplifies the resulting leak-off termand moreover allows us to keep the same leak-off term found in the hydraulic fracturingliterature.With the fifth assumption in mind we refer to the excellent exposition of Carter’sleak-off model found in Adachi’s Ph.D. thesis [2, p.18-23] which leads to the expressiong(X, t) =2CL√t− t0(X, t), t > t0(X, t), (2.15)where CL is Carter’s leak-off coefficient and t − t0 represents the elapsed time since thefracture tip first passed the point currently at a distance x away from the fracture tip. Weremark here that the paper by Linkov [16] considers the problem of a hydraulic fracturewith a more general leak-off term of the formg(X, t) =2CL[t− t0(X, t)]β .Furthermore, by reproducing the derivation found in [2] but using the power-law rheology,it can be shown that β = 1n+1 , though the utility of this added generality is lost in viewof the fifth assumption.2.2.3 Mass Conservation and The Lubrication EquationThe continuity equation for the fracturing fluid including leak-off takes the form∂w∂t+∂q∂X+ g = 0. (2.16)Changing to the moving coordinate x = V t−X we obtain∂∂t= Vddx,∂∂X= − ddx,16Variable x w p E′ V C ′ M ′ K ′Units L L FL−2 FL−2 LT−1 LT−1/2 FTnL−2 FL−3/2Table 2.1: Units of problem variables and parameters (L = length, T = time, F = force).and moreover the elapsed time t− t0(X, t) in Carter’s leak-off model (2.15) is easily seento equal x/V whence the continuity equation becomesVdwdx− dqdx+2CLV1/2√x= 0.Integrating and substituting (2.14) in for q we obtainV w −(w2n+1M ′dpdx)1/n+ 4CLV1/2x1/2 = 0,from which we readily obtain the lubrication equationdpdx=M ′V nwn+1(1 +4CLV 1/2x1/2w)n. (2.17)2.3 Summary of Governing EquationsWe conclude this section by summarizing the governing equation and relevant parameterdefinitions. For convenience, we introduce the parametersM ′ ≡ 2n+1(2n+ 1n)nM, K ′ ≡ 8√2piKIc, C′ ≡ 2CL.With these parameters the governing equations are summarized below.Lubrication Equation:dpdx=M ′V nwn+1(1 +2C ′V 1/2x1/2w)n. (2.18)Elasticity Equation:p(x) =E′4pi−∫ ∞0dwdsdsx− s. (2.19)17Fracture Tip Equation (with Propagation Criterion):w(x) ∼ K′E′x1/2 +O(x3/2) (x→ 0+). (2.20)Finally the units of all problem parameters and variables are summarized in Table2.1 above.18Chapter 3Problem Scaling and AsymptoticRegimesA ubiquitous step in the analysis of mathematical models is to non-dimensionalize itsgoverning equations. This procedure has many benefits such as reducing the number ofparameters and determining scale invariants of the underlying problem. In the context ofhydraulic fractures non-dimensionalization yields another benefit – it leads to a substantialsimplification of the governing equations in three distinct parameter regimes. In each ofthese parameter regimes the governing equations are simplified to an extent that closed-form solutions can be determined. Moreover, each regime can be associated with a distinctlength scale thereby allowing us to associate each of the closed form solutions with thebehaviour of the solution at a particular length scale. This final point is the most importantfor subsequent discussions. Although we will not use the scalings introduced in this chapterfor subsequent ones, we will be making use of the behaviour of solutions when x 1 andwhen x 1.Before proceeding with the non-dimensionalization it is instructive to keep track19of the units of each variable and parameter in our problem for which we refer the readerto Table 2.1. Next we introduce the dimensionless variables ξ, Ω, and Π defined byξ = x/l, Ω = w/l, Π = p/E′.Substituting this into the (rearranged) lubrication equation (2.18) and the tip constraint(2.20) we obtaindΠdξ=((M ′)1nV(E′)1n l1Ωn+1n+2(M ′)1nC ′V 1/2(E′)1n l3/2ξ1/2Ω2n+1n)n,Ω ∼ K′E′l1/2ξ1/2 +O(ξ3/2) (ξ → 0+).We note that the scaling leaves the elasticity equation (2.19) invariant except for cancellingout the Young’s modulus E′. In the above dimensionless equations we recognize threedimensionless groupsGm ≡ (M′)1nV(E′)1n l, Gk ≡ K′E′l1/2, Gc ≡ 2(M′)1nC ′V 1/2(E′)1n l3/2=2GmC ′V 1/2l1/2. (3.1)We make two observations. First the dimensionless groups Gm and Gk are directly propor-tional to M ′ and K ′ respectively. Second, we observe that the dimensionless group Gc isproportional not only to C ′ but also to Gm and therefore M ′. This second observation isgoing to play an important role when we consider different scaling regimes since Gm = 0will always imply Gc = 0. With these dimensionless groups defined as in (3.1) we obtainthe dimensionless equationsdΠdξ=( GmΩn+1n+ Gc ξ1/2Ω2n+1n)n(3.2)Π =14pi−∫ ∞0Ω′(ζ)ξ − ζ dζ, (3.3)andΩ ∼ Gkξ1/2 +O(ξ3/2) (ξ → 0+). (3.4)20At this stage we note that the length scale l remains undefined. This additionaldegree of freedom allows us to set one of the dimensionless groups appearing in (3.1) tounity and solve for l. Hence each parameter M ′, K ′, and C ′ is directly associated with alength scale. Individually setting one of Gm, Gk, and Gc to unity we therefore obtain thethree distinct length scales given bylm =(M ′)1nV(E′)1n, lk =(K ′E′)2, lc =(2(M ′)1nC ′V 1/2(E′)1n)2/3. (3.5)By setting one of the dimensionless groups in (3.1) to unity and the remaining onesto zero we can explore the behaviour of solutions when the propagation of the hydraulicfracture is dominated by one of three processes – viscosity (M), toughness (K), or leak-off(C). The remainder of this chapter explores the solution structure within each of theseparameter regimes.Before continuing to study each propagation regime we make note of the specialcollection of solutions to (3.3) for each 0 < α < 1 given by [14]Ω(ξ) = βξα, Π(ξ) = δξα−1, (3.6)whereδ =αβ4cotpiα. (3.7)3.1 Viscous Regime (M)In the viscous regime we have Gm = 1, Gk = 0, and Gc = 0. The natural length scale isthen given by l = lm. We seek solutions of the form (3.6) whence the lubrication equation(3.2) yields(α− 1)δξα−2 = 1βn+1ξ(n+1)α.21Collecting like powers of ξ we deduce that α = 2/(n+2) and then combining the resultingequation with (3.7) allows us to solve for β and δ. Summarizing the final result we haveΩm ∼ βmξαmm , Πm ∼ δmξαm−1m , (3.8)whereαm =2n+ 2, βm =(−2(n+ 2)2n cot(piαm)) 1n+2, δm = −n+ 2nβ−(n+1)m , (3.9)and where we have used a subscript m to emphasize the viscous scaling l = lm.3.2 Toughness Regime (K)For the toughness regime we have Gk = 1, Gc = 0, and Gm = 0. In this case the naturallength scale is then given by l = lk. Seeking solutions of the form (3.6) satisfying (3.7),the tip behaviour (3.4) immediately implies that β = 1 and α = 1/2. Then (3.7) impliesthat δ = 0 so to leading order the solution pair in the toughness regime is given byΩk0 = ξ1/2, Πk0 = 0.Assuming Gm,Gc  1 we can deduce the leading order correction to Π by substi-tuting Ω = Ωk0 into the lubrication equation (3.2) which givesdΠk1dξ= (Gm + Gc)nξ−n+12 .This equation is readily integrated, making special note of the logarithmic form of thesolution when n = 1. Summarizing, when Gc and Gm are small, the leading order solutionin the toughness regime is given byΩk ∼ ξ1/2k , Πk ∼21−n(Gm + Gc)nξ1−n2k , n 6= 1,(Gm + Gc) log ξk, n = 1,(Gc,Gm  1), (3.10)where we’ve used the subscript k to remind the reader we are using the scaling l = lk.223.3 Leak-Off Regime (C)In the leak-off regime we set Gc = 1 and use the length scale l = lc. As alluded toearlier, since Gm = 0 =⇒ Gc = 0, the present case requires that Gm 6= 0 and we thereforecan’t proceed as in the previous two cases. Instead we proceed by setting 0 < Gm  1and Gk = 0 and seeking a solution of the formΩc ∼ GγmΩcm + ..., Πc ∼ GγmΠcm + ...,where we have used the elasticity equation (3.3) to deduce the same power γ for the firstorder correction. In the absence of toughness the square-root asymptote predicted bythe propagation criterion (3.4) no longer applies so we instead focus on the lubricationequation (3.2), which, to leading order, becomesGγmdΠcmdξ=(G1−γn+1nm1Ωn+1ncm+ G−γ2n+1nmξ1/2Ω2n+1ncm)n. (3.11)The dominant balance in the small Gm is between the term on the left-hand-side and thesecond one on the right-hand-side which yields γ = 0 and the dominant equationΩ2n+1cmdΠcmdξ= ξn/2.Substituting a solution of the form (3.6) we immediately deduce the values of each pa-rameter α, β, and δ whence we obtain the solutionΩc ∼ βcmξαcmc + o(1), Πc ∼ δcmξαcm−1c + o(1) (Gm  1), (3.12)whereαcm =n+ 44(n+ 1), βcm =( −64(n+ 1)23n(n+ 4) cot(piαcm)) 12(n+1), (3.13)and δcm is determined by using (3.7). We remark that, as in the previous two cases, wehave used a subscript of c to remind the reader of the scaling l = lc.230 0.2 0.4 0.6 0.8 1n0.,Fracture Width Power (,) versus Power Law Index (n),k = 1=2,cm = (n+ 4)=(4n+ 4),m = 2=(n+ 2)Figure 3.1: Fracture width power α in Ω ∼ βξα for each parametric regime: αk toughness,αcm leak-off, and αm viscous regime.3.4 Comments on Negative PressureWe begin by discussing the form of the dimensionless pressure within the viscosity andleak-off regimes. Since in each case 1/2 < α < 1 we deduce that cot(piα) < 0 and thepressure is therefore negative. This is best understood by considering the scenario of adry semi-infinite fracture for which we have the classical solution Ω ∝ ξ1/2 for the fractureopening [20, 12]. Within the viscosity and leak-off regimes the bound α > 1/2 indicatesthat the fracture is narrower than in a dry fracture, meaning that the fluid is exerting anegative pressure along the opening boundaries [2].3.5 Comment on the n = 0 CaseIt is important to keep in mind that the n = 0 case does not represent a physicallyrealistic system, and should rather be interpreted as a limiting solution for n → 0+. Inthis case the χ dependence vanishes from the governing equations, and as a result theleak-off regime also vanishes. The toughness solution remains unchanged. On the other24hand, ignoring the tip-asymptote we obtain the viscous solution w =√4piM ′/E′x.3.6 Spatial Ordering of Parameter RegimesAnalysing each parameter regime is not only instructive in understanding the solutionstructure in the three limiting cases, but it also allows us to determine the length scalesat which each parameter regime dominates. In Figure 3.1 we observe the orderingαk < αcm < αm, 0 < n ≤ 1.This suggests that for ξ  1 the toughness (K) regime will dominate, whereas for ξ  1the viscosity (M) regime will dominate. Furthermore the leak-off regime (C) regimedominates in some intermediate region. The resulting ordering is neatly summarized inthe schematic of a fracture depicted in Figure 3.2.K C MxFigure 3.2: Schematic showing ordering of regions where each parameter regime dominates(K: toughness, C: leak-off, M: viscosity).25Chapter 4Non-Singular Integral EquationFormulationSolutions to hydraulic fracture problems are typically obtained by simultaneously solv-ing (2.18), (2.19), and (2.20) numerically. One of the biggest difficulties with this approachis dealing with the singular integral appearing in (2.19). Following the approach of Dontsovet. al. in [6] we will reformulate the current system of equations in the unknowns w(x)and p(x) as one non-singular and non-linear integral equation in the single unknown w(x).We proceed by inverting the elasticity equation (2.19) using the inversion formulafound in [9]w(x) =K ′E′x1/2 +4piE′∫ ∞0K(x, s)p(s) ds, (4.1)where the kernel K(x, s) is given byK(x, s) = log∣∣∣∣x1/2 + s1/2x1/2 − s1/2∣∣∣∣− 2x1/2s1/2 . (4.2)Observing that K(x, s) = ∂F (x,s)∂s whereF (x, s) = (s− x) log∣∣∣∣x1/2 + s1/2x1/2 − s1/2∣∣∣∣− 2x1/2s1/2, (4.3)26we may integrate the rightmost term in (4.1) by parts to get∫ ∞0K(x, s)p(s) ds = lims→∞F (x, s)p(s)− lims→0+ F (x, s)p(s)−∫ ∞0F (x, s)p′(s) ds.Provided we can evaluate the two limit terms resulting from the integration by parts,substituting the lubrication equation (2.18) in for p′(s) will yield a non-linear integralequation in the single unknown function w(x). We readily obtainF (x, s) ∼ −4x1/2s1/2 +O(s3/2), (s→ 0+),F (x, s) ∼ −43x3/2s−1/2 +O(s−3/2), (s→∞).Furthermore, the parameter regime ordering determined in the previous chapter and morespecifically the form of Π(ξ) in (3.10) and (3.9) imply thatp(s) ∼O(s1−n2 ), if n 6= 1,O(log s), if n = 1,(s→ 0+),p(s) ∼ O(s− nn+2 ), (s→∞).Therefore both of the limits from the integration by parts vanish and after substitutingthe lubrication equation (2.18) in for p′(s), the inverted elasticity equation (4.1) becomesw(x) =K ′E′x1/2 +4M ′V npiE′x1/2∫ ∞0s1/2G(s1/2/x1/2)w(s)n+1(1 +2C ′V 1/2s1/2w(s))nds, (4.4)where the kernel G(t) is given byG(t) =1− t2tlog∣∣∣∣1 + t1− t∣∣∣∣+ 2. (4.5)We can further simplify (4.4) by introducing an appropriate scaling. Accordinglywe introduce the non-dimensional variables w˜ = w/wˆ and x˜ = x/xˆ whence (4.4) becomesw˜(x˜) = gkx˜1/2 +4pix˜1/2∫ ∞0s˜1/2G(s˜1/2/x˜1/2)w˜(s˜)n+1(gm + gcs˜1/2w˜(s˜))ns˜,27where we have defined the non-dimensional groupsgm ≡ (M′)1nV xˆ2n(E′)1n wˆn+2n, gk ≡ K′xˆ1/2E′wˆ, gc ≡ 2(M′)1nC ′V12 xˆn+42n(E′)1n wˆ2(n+1)n. (4.6)Note that when wˆ = xˆ = l each of these dimensionless groups reduce to the correspondingoriginal dimensionless groups (3.1). Furthermore, from each dimensionless group we candetermine a self-similar form for the solution on a different length scale. Indeed settingeach dimensionless group gk, gc and gm to unity we respectively obtainwˆm =(M ′V nE′) 1n+2xˆ2n+2m , wˆk =K ′E′xˆ1/2k , wˆc =(2nM ′(C ′)nVn2E′) 12(n+1)xˆn+44(n+1)c .Each such expression in turn suggests, respectively, the self-similar formw˜ =E′M ′V nwn+2x2, w˜ =E′K ′wx1/2, w˜ =E′2nM ′(C ′)nVn2w2(n+1)xn+42.Motivated by the analysis found in [6] we seek a solution of the form implied by thesecond self-similar form and therefore introduce the change of variables w˜ = E′K′w/x1/2,x˜ = x1/2/l1/2. The other two self-similar forms, while also valid, provide no significantsimplification in the governing equations. Substituting this into (4.4) givesw˜(x˜) = 1 +8piM ′V n(E′)n+1(K ′)n+2l2−n2∫ ∞0G(s˜/x˜)s˜n−1w˜(s˜)n+1(1 +2C ′E′V 1/2K ′1w˜(s˜))nds˜,within which we can choose l and χ to reduce the number of parameters. Summarizingthe results, we introduce new dimensionless variables w˜ and x˜ as well as a dimensionlessparameter χ defined byw˜ =E′K ′wx1/2, x˜ =x1/2l1/2, l ≡[(K ′)n+2M ′V n(E′)n+1] 22−n, χ ≡ 2C′E′V 1/2K ′, (4.7)so that (4.1) becomesw˜(x˜) = 1 +8pi∫ ∞0G(s˜/x˜)s˜1−nw˜(s˜)1+n(1 +χw˜(s˜))nds˜. (4.8)28Moreover with this scaling the propagation criterion (2.20) becomes w˜(0) = 1. This non-singular and non-linear integral equation will be the focus of the remaining discussions. Inparticular we will discuss how to obtain a first order ODE approximating (4.8), followedby brief discussion of the numerical solution to this ODE as well as some special cases forwhich an exact solution can be found.4.1 Recovering Asymptotic Regimes from the Non-SingularIntegral EquationIn this section we demonstrate that the solutions in each parameter regime obtained inChapter 3 can be recovered directly from the non-singular integral equation (4.8). Thatthis is the case for the toughness regime is immediately clear from the scaling (4.7) andthe fact that w˜(0) = 1. To deduce the remaining asymptotic solutions we will make useof the identity (see Appendix A.1)∫ ∞0G(t)tαdt =2piα(2− α) tan(piα2), (−1 < α < 1). (4.9)We consider first the case of w˜  χ. Then the integral equation (4.8) becomesw˜(x˜) ∼ 8pi∫ ∞0G(s˜/x˜)s˜1−nw˜(s˜)1+nds˜.Seeking a monomial solution (reminiscent of the solution pairs used in 3) of the formw˜(x˜) = β˜x˜α˜ and using (4.9) we deduceα˜ =2− n2 + n, β˜ =[−2(2 + n)2ntan(2pi2 + n)] 12+n= βm. (4.10)After reversing the scaling (4.7) we find that this solution corresponds to the viscoussolution described by (3.8) and (3.9).29Next we consider the case of w˜  χ so that, similarly to the the previous case, theintegral equation (4.8) becomesw˜(x˜) ∼ 8piχn∫ ∞0G(s˜/x˜)s˜1−nw˜(s˜)1+2nds˜.This time the monomial solution w˜(x˜) = β˜x˜α˜ yieldsα˜ =2− n2(1 + n), β˜ =[−64(1 + n)23n(4 + n)tan(pi(4 + n)4(1 + n))] 12(1+n)χn2(1+n) = βcmχn2(1+n) ,(4.11)which we can once again equate to the leak-off solution given by (3.12) and (3.13).30Chapter 5Reduction to an ApproximatingFirst-Order ODEIn this chapter we will follow the procedure set out by Dontsov and Peirce in [6] toconstruct an ODE that approximates solutions to the non-singular integral equation (4.8).As briefly mentioned in the introduction, we are motivated to seek an approximate solutionin order to cut down on computation time. While we may solve the full equation (4.8)numerically, the time required to do so results in multiscale hydraulic fracture simulatorsthat are impractical. We will consider the cases 0 < n ≤ 1 and n = 0 separately, startingwith the former. The resulting approximating ODE will be seemingly more complicatedthan the original non-linear integral equation, but we will be able to successively constructexplicit solutions for it of the form x˜ = x˜(w˜) where x˜(w˜) will be given by an integral thatis fast to calculate using standard quadrature rules.The key feature of the non-singular integral equation (4.8) that will allow us toapproximate it by an ODE is the sharp transition from 4 to 1 that G(t) undergoes neart = 1 on a logarithmic scale (see left pane of Figure 5.1). This suggests that G′(t) will3110-10 100 1010t01234G(t)10-10 100 1010t-6-5-4-3-2-101G0(t)Figure 5.1: Plots of kernel function G(t) and its first derivative G′(t).have a sharp peak near t = 1 on a logarithmic scale. Indeed, differentiating we obtainG′(t) = −1 + t2t2log∣∣∣∣1 + t1− t∣∣∣∣+ 2t , (5.1)which is plotted in the right pane of Figure 5.1 and reveals a sharp peak at t = 1. Thus,if f(t) varies slower than G′(t) near t = 1, then the preceding discussion motivates theapproximation ∫ ∞0G′(t)f(t)g(t) dt ≈ f(1)∫ ∞0G′(t)g(t) dt, (5.2)where g(t) is any function for which the integral on the right exists.Returning to the construction of an approximating ODE, we begin by differentiat-ing both sides of (4.8) and rearranging terms to obtaindw˜dx˜= − 8pix˜2∫ ∞0G′(s˜/x˜)(1 + χ/w˜(s˜))1−ns˜2−nw˜(s˜)1+nds˜− 8χpix˜2∫ ∞0G′(s˜/x˜)(1 + χ/w˜(s˜))1−ns˜2−nw˜(s˜)2+nds˜.(5.3)Similarly to the discussion in §4.1 we deduce that when w˜  χ, the first term on theright-hand-side will dominate, whereas when w˜  χ it is the second term that dominates.Restricting first our attention to those value of x˜ where w˜(x˜) χ, the sharp peak of G′(t)32at t = 1 implies that the dominant contribution to the first integral on the right-hand-sidewill come from s˜ ≈ x˜ whence within that integral w˜(s˜)  χ. Hence by the discussionin §4.1 we will have w˜(s˜) = β˜s˜δ˜ where β˜ and δ˜ are both functions of s˜ that don’t varysignificantly from the corresponding values found in (4.10). In particular, quantities ofthe form w˜(s˜)(s˜/x˜)δ˜will not vary significantly for s˜ ≈ x˜. In addition since (1 + χ/w˜(s˜)) ≈ 1we can use (5.2) to obtain the following approximation∫ ∞0G′(s˜/x˜)(1 + χ/w˜(s˜))1−ns˜2−nw˜(s˜)1+nds˜ ≈ x˜2−n(1 + χ/w˜(x˜))1−n∫ ∞0G′(s˜/x˜)(s˜/x˜)n−2+(1+n)δ˜[(s˜/x˜)δ˜w˜(s˜)]1+nds˜≈ x˜3−n/w˜(x˜)1+n(1 + χ/w˜(x˜))1−n∫ ∞0G′(t)tn−2+(1+n)δ˜dt.(5.4)To determine an appropriate approximation of the second integral appearing in(5.3) we proceed similarly. We restrict x˜ to the region where w˜(x˜)  χ and use thebehaviour of G′(t) to deduce that the dominant contribution in the second integral comesfrom s˜ ≈ x˜ where as a result w˜(s˜)  χ. In contrast to the previous case we now have(1 + w˜(s˜)/χ)1−n ≈ 1 so the approximation becomes∫ ∞0G′(s˜/x˜)(1 + χ/w˜(s˜))1−ns˜2−nw˜(s˜)2+nds˜ ≈ χn−1x˜2−n(1 + w˜(x˜)/χ)1−n∫ ∞0G′(s˜/x˜)(s˜/x˜)n−2+(1+2n)δ˜[(s˜/x˜)δ˜w˜(s˜)]1+2nds˜≈ x˜3−n/w˜(x˜)2+n(1 + χ/w˜(x˜))1−n∫ ∞0G′(t)tn−2+(1+2n)δ˜dt.(5.5)Combining (5.4) and (5.5) with (5.3) leads to the approximating ODE for 0 < n ≤ 1dw˜dx˜=C1(δ˜, n)(1 + χ/w˜)1−nx˜1−nw˜1+n+ χC2(δ˜, n)(1 + χ/w˜)1−nx˜1−nw˜2+n, w˜(0) = 1, (5.6)whereC1(δ˜, n) ≡ − 8pi∫ ∞0G′(t)tn−2+(1+n)δ˜dt, (5.7)C2(δ˜, n) ≡ − 8pi∫ ∞0G′(t)tn−2+(1+2n)δ˜dt. (5.8)33We proceed similarly for n = 0 with the added simplification that the problem isnow independent of χ. We differentiate both sides of (4.8) when n = 0 to obtaindw˜dx˜= − 8pix˜2∫ ∞0G′(s˜/x˜)s˜2w˜(s˜)ds˜.This time the approximation (5.2) takes the form∫ ∞0G′(s˜/x˜)s˜2w˜(s˜)ds˜ = x˜2∫ ∞0G′(s˜/x˜)(s˜/x˜)δ˜−2[(s˜/x˜)δ˜w˜(s˜)]ds˜≈ x˜3w˜(x˜)∫ ∞0G′(t)tδ˜−2dt,where we have assumed the term in square brackets does not undergo a significant changenear s˜ = x˜. The resulting approximating ODE is therefore given bydw˜dx˜= C0(δ˜)x˜w˜, w˜(0) = 1, (5.9)whereC0(δ˜) ≡ − 8pi∫ ∞0G′(t)tδ˜−2dt, (5.10)At this stage it is worthwhile to discuss the method used to obtain the approxi-mating ODEs, as well as the ODEs themselves. First, it is important to remark that ineach of the definitions (5.7), (5.8), and (5.10) the variable δ˜ is a function of tx˜. Moreover,writing w˜(x˜) = β˜(x˜)x˜δ˜(x˜) and assuming that β˜ and δ˜ are slowly varying we obtain theapproximationδ˜ ≈ x˜w˜dw˜dx˜,which shows that in fact δ˜ is a function of x˜, w˜, and d ˜˜w/dx˜. Thus the terms C1 and C2appearing in (5.6) as well as C0 appearing in (5.9) represent non-linear terms. This bringsinto question the benefit of using this approximation. Motivated by the next point, in thefollowing section we will construct solutions by first assuming δ˜ is constant, solving theODE, and then obtaining corrections iteratively.34The second point we make involves the suitability of using the approximation (5.2)in the above construction. This was briefly touched on in the 0 < n ≤ 1 where we usedthat w˜(s˜) was either in the viscous or the leak-off regimes and was therefore suitablyapproximated by a function of the form w˜ = β˜x˜δ˜, where β˜ and δ˜ did not vary significantlyfrom their counterparts in (4.10) and (4.11) respectively. However, this does not addressthe question of the suitability of this approximation outside of these spatial regions. Forexample, we may question in what sense (5.2) is suitable in (5.4) when x˜ is not in theviscous regime. A similar question can be raised for the n = 0 case where we simply statedour assumption that (s˜/x˜)δ˜w˜ does not vary significantly.To address this question we turn to our discussion of the asymptotic regimes foundin §3 and in particular the spatial ordering of the parameter regimes. This discussionrevealed that provided χ 6= 0 and n 6= 0 there will be three-distinct parameter regionscorresponding to toughness, leak-off, and viscosity. Within each such region the solutiontakes the form of a monomial. The same is true for χ = 0 or n = 0 where there are nowonly two distinct parameter regions. The upshot of this is that within these regions theassumption that w˜(s˜)(s˜/x˜)δ˜is approximately constant is appropriate. Moreover, this discussionalso reveals that this approximation is not appropriate in the transition regions betweenparameter regimes where β˜ and δ˜ may undergo sufficiently rapid changes. Therefore, weshould expect the errors of the approximating ODE to be concentrated in the transitionregions.5.1 The Constant δ˜ ApproximationThe concluding comments in the section above revealed two key points worth pursuingfurther. First, in their current form, the approximating ODEs (5.6) and (5.9) provideno clear advantage due to the non-linearities arising from δ˜ ≈ x˜w˜′/w˜. Secondly, the350 0.5 1~/105C0Figure 5.2: Variation in C0 as δ˜ varies between 0 ≤ δ˜ ≤ 1.approximations used to obtain these ODEs are only suitable in those regions where thesolution takes the form of a monomial w˜ = β˜x˜δ˜ (more specifically in the toughness, leak-off, and viscosity dominating regions). Motivated by these two points, we proceed byfixing the value of δ˜ within each of the integrals in (5.7), (5.8), and (5.10) to a functiononly of x˜. Then using the result (see Appendix A.2)∫ ∞0G′(t)tκdt =2piκ(1 + κ)(1− κ) tan(pi2(κ+ 1))(−2 < κ ≤ 0), (5.11)we obtain the following closed form expressions for (5.10), (5.7), and (5.8)C0(δ˜) =16(2− δ˜)(1− δ˜)(3− δ˜) tan(pi2(1− δ˜)), (5.12)C1(δ˜, n) =16(2− n− (1 + n)δ˜)(1− n− (1 + n)δ˜)(3− n− (1 + n)δ˜) tan(pi2(1− n− (1 + n)δ˜)), (5.13)C2(δ˜, n) =16(2− n− (1 + 2n)δ˜)(1− n− (1 + 2n)δ˜)(3− n− (1 + 2n)δ˜) tan(pi2(1− n− (1 + 2n)δ˜)), (5.14)where we emphasize that δ˜ = δ˜(x˜). The use of (5.11) is readily seen to be valid for0 < n ≤ 1 and 0 < δ˜ ≤ 2−n2+n since in this case −2 < κ ≤ 0 is satisfied for each integralappearing in C0, C1, and C2. The same can be said for n = 0 when 0 < δ˜ ≤ 1, but when360 0.5 1~/203040506070C1(a)0 0.5 1~/203040506070C2n=0.1n=0.2n=0.4n=0.6n=0.8n=1.0(b)Figure 5.3: Variation in coefficients C1 and C2 as δ˜ varies between 0 ≤ δ˜ ≤ 2−n2+n for several0 < n ≤ 1 values.δ˜ = 0 we find that the integral appearing in C0 has κ = −2 and the identity (5.11) nolonger applies. Fortunately, in this case we have x˜ = 0 so the coefficient C0 disappearsfrom the ODE (5.9) altogether.In Figures 5.2 and 5.3 we see how each of the coefficients C0, C1, and C2 vary withδ˜. These figures indicate that as n decreases the range of values taken by each coefficientincreases, and in particular C0 diverges as δ˜ → 0 as expected. When n ≈ 1 we observethat there is a relatively small change in both coefficients C1 and C2 as δ˜ varies. In thepaper by Dontsov and Peirce [6] this relatively small variation motivated setting C1 andC2 to two distinct constants that accurately capture individual propagation regimes. Wewill take the same approach for 0 ≤ n ≤ 1, and although, in light of Figures 5.2 and 5.3,this is a rather extreme assumption, comparisons to full numerical results in Chapter 6will reveal that this is an acceptable approximation.Having decided to fix the values of the coefficients C0, C1 and C2, or equivalentlythe values of δ˜ ∈ [0, 2−n2+n ], the question naturally arises as to what values of δ˜ should be37chosen. In addressing this question we first focus on the 0 < n ≤ 1 case. We observe thatthere are two distinguished limits we can deduce from the 0 < n ≤ 1 ODE (5.6). One suchlimit is when w˜  χ and therefore the first term dominates (which isolates C1), whereasfor the other limit w˜  χ and therefore the second term dominates (which isolates C2).By comparing the resulting behaviour of w˜ implied by (5.6) in these two cases with thatimplied by the non-linear integral equation (4.8) we can therefore determine what valuesC1 and C2 should take, which in turn will motivate our decision of δ˜ values.We focus now on the distinguished limit w˜  χ. Then the ODE (5.6) and itsresulting solution becomedw˜dx˜∼ C1(δ˜, n) x˜1−nw˜1+n−→ w˜ =(1 +2 + n2− nC1(δ˜, n)x˜2−n) 12+n.On the other hand, as seen in §4.1 the solution to the non-singular integral equation (4.8)has the form w˜ = β˜x˜α˜ where β˜ and α˜ are given by (4.10). Comparing this to the solutionobtained from the ODE and noting that the viscous regime dominates for x˜ 1 we deducethe requirement thatC1(δ˜, n) =2− n2 + nβ2+nm = −2(4− n2)ntan(2pi2 + n). (5.15)As is easily checked using (5.13) this requirement is met by setting δ˜ = 2−n2+n . Note that inthis case the first term in the ODE accurately captures the viscous solution (w˜  χ)) byexhibiting both the correct coefficient β˜ and power α˜ in (4.10).Proceeding as above, we now consider the limit w˜  χ in which case the ODE(5.6) and it’s solution becomedw˜dx˜∼ C2(δ˜, n)χn x˜1−nw˜1+2n−→ w˜ =(1 + 21 + n2− nC2(δ˜, n)χnx˜2−n) 12(1+n).In this regime the solution of the non-singular integral equation (4.8) takes the formw˜ = β˜x˜α˜ where the coefficient and power are now given by (4.11). Comparing this38solution to that obtained from the ODE we deduce the requirement thatC2(δ˜, n) =2− n2(1 + n)β2(1+n)cm = −32(2− n)(1 + n)3n(4 + n)tan(pi(4 + n)4(1 + n)). (5.16)Using (5.14) we deduce that this requirement is met upon setting δ˜ = 2−n2(1+n) , whichcorresponds to the correct value of α˜ from (4.11).To complete this section we now focus on the n = 0 ODE (5.9). In this case withδ˜ fixed the solution is immediately seen to bew˜(x˜) = (1 + C0(δ˜)x˜2)12 .Motivated by the previous discussion we match this up with the viscous solution which inthe limit n→ 0 takes on the value w˜(x˜) = 2√pix˜ and hence we requireC0(δ˜) = 4pi, (5.17)which we see is obtained by taking the limit δ˜ → 1 in (5.12).In summary, when 0 < n ≤ 1 we obtain the ODEdw˜dx˜=γ1(1 + χ/w˜)1−nx˜1−nw˜1+n+γ2χ(1 + χ/w˜)1−nx˜1−nw˜2+n, w˜(0) = 1, (5.18)whereγ1 = −2(4− n2)ntan(2pi2 + n), (5.19)γ2 = −32(2− n)(1 + n)3n(4 + n)tan(pi(4 + n)4(1 + n)). (5.20)Similarly, for n = 0 we have the ODEdw˜dx˜= 4pix˜w˜, w˜(0) = 1. (5.21)395.2 Solutions of the Constant δ˜ ODEThe simplifications to the original ODEs (5.6) and (5.9) granted by fixing δ˜ are tremen-dous. However, even with these simplifications, apart from the special cases of n = 0 orχ = 0 the resulting ODE (5.18) does not admit an explicit solution of the form w˜ = w˜(x˜).We can overcome this difficulty in two different ways. The first is to use a numerical solversuch as Matlab’s ode45 to solve the IVP (5.18) numerically for x˜ in some finite range (say0 ≤ x˜ ≤ 1020). For example, upon specifying the grid points xt we may determine thesolution wta through the commands% calcualte ODE coefficientsg1 = -2*(4-n^2)*tan(2*pi/(2+n))/n;g2 = -32*(2-n)*(1+n)*tan(pi*(4+n)/(4*(1+n)))/(3*n*(4+n));% solve the constant delta approximating ODE using ode45odefun = @(x,w) ((1+chi./w).^(n-1)).*(g1./w.^(1+n) + chi*g2./w.^(2+n)).*x.^(1-n);[~,wta] = ode45(@(X,W)odefun(X,W),xt,1);This approach has the benefits of being sufficiently fast, and moreover it allows the userto determine w˜ values at specified x˜ values.In the second approach we discard the requirement that the solution be of the formw˜ = w˜(x˜) and instead use the fact that the ODE is separable to obtain a solution of theform x˜ = x˜(w˜) given byx˜ =[(2− n)∫ w˜1(v + χ)1−nv1+2nγ1v + γ2χdv] 12−n. (5.22)This expression has the benefit that by a suitable change of variables we can map the inter-val of integration from [1, w˜] to [−1, 1] thereby making it suitable for standard quadraturerules. The downside to this approach is that to obtain w˜ = w˜(x˜) we would need to solvean implicit equation. Alternative we could calculate x˜ for a large sample of points w˜and from this generate an interpolating polynomial solution using, for example, Matlab’spchip function. In some cases, such as the n = 1 case discussed in [6], the integral ap-40pearing in (5.22) can be evaluated explicitly. However even in this case there is littleadvantage in evaluating the integral explicitly when a standard quadrature scheme can bedone at the same or faster speed and with sufficiently high accuracy.When n = 0 or χ = 0 we are in the more favourable situation wherein we can solvethe constant δ˜ ODE explicitly. In the case of arbitrary 0 ≤ n ≤ 1 and χ = 0 we havew˜(x˜) =(1 +2 + n2− nγ1x˜2−n) 12+n(5.23)whereas when n = 0 we havew˜(x˜) = (1 + 4pix˜2)1/2 (5.24)5.3 The δ˜-Corrected SolutionIn this section we explore the so-called “δ˜-correction” procedure set out in [6]. Thisintroduces a correction to the constant δ˜ approximation by iteratively solving for anupdated δ˜ via the approximation δ˜ ≈ x˜w˜ dw˜dx˜ . We begin by focusing on the 0 < n ≤ 1 case.As the “zero-order” solution we take the constant δ˜ solution (5.22) which we repeat herefor conveniencex˜0(w˜) =[(2− n)∫ w˜1(v + χ)1−nv1+2nv + b0χdv] 12−n(5.25)where we have defined b0 ≡ γ2/γ1. Substituting this into the constant coefficient ODE(5.6) we obtain the first δ˜ correctionδ˜1 =x˜0(w˜)w˜dw˜dx˜∣∣∣∣x˜=x˜0=γ1(1 + χ/w˜)x˜0(w˜)2−nw˜2+n(1 +b0χw˜). (5.26)Since we have δ˜1 = δ˜1(w˜) it follows that (5.6), like the constant δ˜ ODE (5.18), is separableand hence we can solve for x˜ to obtain the first δ˜-corrected solutionx˜1(w˜) =[(2− n)∫ w˜1(v + χ)1−nv1+2nC1(δ˜1(v), n)v + C1(δ˜2(v), n)χdv] 12−n. (5.27)41At this stage we remark that in [6], the authors replace the v dependence of δ˜1 withinthe integral, with a dependence on w˜. This renders the C1 and C2 terms constant withinthe integral thereby simplifying the calculation. The reasons for this decision are notclear, but by comparing this to the solution obtained by solving the integral equation(4.8) numerically we find that it gives a better approximation. For this reason, we willbe making use of the same approximation, though we emphasize that the reasons for itsimproved accuracy are not completely understood. Thus, as the first δ˜-corrected solutionwe takex˜1(w˜) =[2− nC1(δ˜1(w˜), n)∫ w˜1(v + χ)1−nv1+2nv + b1(w˜)χdv] 12−n. (5.28)where we have defined b1(w˜) ≡ C2(δ˜1(w˜), n)/C1(δ˜1(w˜), n). By following the same proce-dure we can find the second δ˜ correction and so on, though numerical comparisons revealthat beyond the first δ˜ correction there is a minimal improvement in accuracy. In then = 0 case we obtain analogous results but with the absence of χ terms, and furthermorereplacing C1 with C0.In Appendix B.2 and B.3 we have included sample Matlab code for evaluating thezero and first δ˜-corrected solutions. In this code we have mapped the integration intervalto [−1, 1] and have used a 64-point Gauss-Legendre scheme to evaluate the integral.The next section involves the numerical solution of the full non-singular integralequation (4.8). This is accomplished by first discretizing the integral using Simpson’srule and then solving the resulting non-linear system in the unknowns Wj ≈ w˜(Xj) using(vector) Newton’s method. In this context, solutions to the constant δ˜ will play the roleof the initial guess in Newton’s method. Since we will be interested in finding w˜ values atspecific x˜ values, we will opt to solve the constant δ˜ ODE numerically using ode45. Wewill then compare solutions obtained in this way to those obtained using the zero- andfirst-δ˜ corrected solutions.42Chapter 6Numerical Solution of Non-linearIntegral EquationIn this section we will be solving the non-linear integral equation (4.8) numerically. Webegin by discretizing the positive real line using logarithmically distributed points0 < X1 < · · · < X2N+1 <∞where we assume X1  1 so that w˜(x˜) ≈ 1 for all x˜ ≤ X1 and X2N+1  1 so thatw˜(x˜) ≈ βmx˜2−n2+n for all x˜ ≥ X2N+1. The choice of logarithmatically distributed points isto accommodate the discretization of this large interval. In addition we have chosen touse 2N+1 points in the discretization to facilitate the use of Simpson’s rule in subsequentsteps.Now we evaluate the non-linear integral equation at x˜ = Xi and rewrite it asw˜(Xi) = 1 +8pi(I0(Xi, w˜) + I(Xi, w˜) + I∞(Xi, w˜)), (6.1)43whereI0(Xi, w˜) =∫ X10G(s˜/Xi)s˜1−nw˜(s˜)1+n(1 +χw˜(s˜))nds˜, (6.2)I(Xi, w˜) =∫ X2N+1X1G(s˜/Xi)s˜1−nw˜(s˜)1+n(1 +χw˜(s˜))nds˜, (6.3)I∞(Xi, w˜) =∫ ∞X2N+1G(s˜/Xi)s˜1−nw˜(s˜)1+n(1 +χw˜(s˜))nds˜. (6.4)In what follows we will discretize I(Xi, w˜), and find leading order approximations for I0and I∞ as functions of X1, Xi, and X2N+1 respectively.We introduce the vector W = (W1, ...,W2N+1)T where each Wi approximatesw˜(Xi) by solving the resulting non-linear system. In addition we define Gij ≡ G(Xj/Xi)so that for each i = 1, ..., 2N + 1 the Simpson’s rule discretization of I(w˜,Xi) becomesI(Xi,W) ≈ 16(2− n)N∑j=1(X2−n2j+1 −X2−n2j−1)[Gi,2j−1Wn+12j−1(1 +χW2j−1)n+ 4Gi,2jWn+12j(1 +χW2j)n+Gi,2j+1Wn+12j+1(1 +χW2j+1)n],(6.5)where we remark we have made the substitution s˜1−nds˜ = 12−nd(s˜2−n) prior to usingSimpson’s rule. In particular this substitution allows us to make the terms in squarebrackets depend solely on W.Before approximating I0(Xi,W) and I∞(Xi,W) we make note of the followingseriesG(t) = 4− 4t2∞∑k=0t2k(2k + 1)(2k + 3), (0 ≤ t ≤ 1), (6.6)G(t) =4t2∞∑k=0t−2k(2k + 1)(2k + 3), (1 ≤ t). (6.7)The convergence of the series in the designated intervals is easily obtained by using theratio test and by comparing it to the series∑∞k=1 k−2. To approximate I0(Xi,W) we use44(6.6) and w˜(s˜) ≈ 1 for all s˜ ≤ X1 to obtainI0(Xi) ≈ (1 + χ)n2− n X2−n1[4− 4∞∑k=02− n(2k + 4− n)(2k + 1)(2k + 3)(X1Xi)2k+2]. (6.8)By combining (2 − n)/(2k + 4 − n) ≤ 1 with (6.6) and 0 ≤ G(t) ≤ 4 (see Figure 5.1) weobtain the bounds0 ≤ 4∞∑k=0(2− n2k + 4− n)1(2k + 1)(2k + 3)(X1Xi)2k+2≤ 4−G(X1/Xi) ≤ 4,whence we deduce that the term in square brackets in (6.8) is at most O(1). HenceI0 = O(X2−n1 ) and since X1  1 it follows that I0 is negligibly small and will thereforebe ignored.The approximation of I∞(Xi,W) is more involved than that of I0(Xi,W) andfurthermore results in a non-negligible contribution whenever i = 2N + 1. This time weuse w˜(s˜) ≈ βms˜2−n2+n for all s˜ ≥ X2N+1 and (6.7) to obtainI∞(Xi) ≈∞∑p=04cpβn+1mX1− 2n2+n−p 2−n2+n2N+1∞∑k=0(X2N+1/Xi)−2k−2(2k + 1 + 2n2+n + p2−n2+n)(2k + 1)(2k + 3), (6.9)where each cp is the coefficient of s˜p(2−n)/(2+n) in the series expansion of (1+ χβms˜(2−n)/(2+n))n.Now using (6.7) we have0 ≤ 4∞∑k=0(X2N+1/Xi)−2k−2(2k + 1 + 2n2+n + p2−n2+n)(2k + 1)(2k + 3)≤ G(X2N+1/Xi) ≤ 4,whence the coefficient of each X1− 2n2+n−p 2−n2+n2N+1 in (6.9) is at most O(1) and thereforeI∞(Xi) =4βn+1mX2−n2+n2N+1∞∑k=0(X2N+1/Xi)−2k−2(2k + 1 + 2n2+n)(2k + 1)(2k + 3)+O(X− 2n2+n2N+1). (6.10)It remains to approximate the infinite series making up the coefficient of X2−n2+n2N+1. Wheni < 2N + 1 the terms in this infinite series decay exponentially like (Xi/X2N+1)2k+2 andwe can therefore truncate the sum after one term obtaining∞∑k=0(X2N+1/Xi)−2k−2(2k + 1 + 2n2+n)(2k + 1)(2k + 3)=2 + n3(2 + 3n)(XiX2N+1)2+O((XiX2N+1)4), (6.11)450 10 20 30 40M0. Sum0 10 20 30 40M10-510-410-310-210-1M-Partial Sum Absolute ErrorFigure 6.1: Convergence of partial sums making up coefficient of X2−n2+n2N+1 in (6.10). Theleft pane shows values of the partial sums (solid lines) and infinite sums approximatedwith M = 104 (dashed lines) for values of n = 0 at the top and n = 1 at the bottom in0.10 steps. The right pane shows the absolute error of partial sums (note that this plot isactually showing a small band that includes all the n values plotted in the left pane).which has a small error due to the logarithmic distribution of the sample points. Wheni = 2N+1 the terms in the series decay algebraically like k−3 which implies that truncatingthe sum after the first M terms will yield an error that is O(M−2). This means that toget an acceptable error in the coefficient, say of the order of 10−2, we will need to take atleast 10 terms before truncating the sum (see Figure 6.1). Since I∞(Xi) is independentof W we can compute it prior to solving the resulting non-linear system, thereby savingcomputation time.With I(Xi,W) given by (6.5), I0(Xi) ≈ 0, and I∞(Xi) given by appropriatelytruncating the sum in (6.10) we therefore obtain the non-linear system of equationsAi ≡ 1 + 8pi(I(Xi,W) + I∞(Xi))−Wi = 0, i = 1, ..., 2N + 1, (6.12)which we can solve using the vector Newton’s method. The initial guess for Newton’s46method is taken to be the solution of the constant δ˜ ODE discussed in §5.1. Due to thelarge values we expect the later entries of W to take, we will be using a relative changein the solution between iterations for the stopping criterion. The details of the numericalcomputation of w˜ using this method can be found in the code in Appendix B.1. In thiscode we see that we have chosen M = 20 when calculating the coefficient from (6.10),which according to Figure 6.1 has an associated error of approximately 10−4.In the following section we will explore the performance of this numerical solverthrough a variety of numerical experiments. This will be followed by a brief overviewof the qualitative characteristics of the solutions obtained by this method. In the finalsection we will use the solutions of the full non-singular integral equation to comparethe accuracy of the solutions obtained by solving the approximating ODE using each ofthe three methods found in §5: the solution found using ode45, and the zero and first δ˜corrected solutions.6.1 Performance of Non-Singular Integral Equation SolverWe will gauge the performance of the numerical solver discussed above by exploringhow three quantities – the total computation time, maximum number of iterations, andmaximum residual – vary with the numerical parameter N , and the problem parametersn and χ. This was done by solving the non-linear system with Newton’s method where wehave chosen X1, ..., X2N+1 to consist of 2N + 1 logarithmically distributed points between10−10 and 1015. In each simulation we have fixed the tolerance (i.e. the maximum relativedifference between iterations at which Newton’s method stops) to 10−3, and the maximumnumber of iterations to 500. Moreover, for each value of n = 0, 0.1, 0.2, ..., 1 we havecalculated the maximum of each of the three performance parameters over the range of10−5 ≤ χ ≤ 105.470 0.2 0.4 0.6 0.8 1n10-410-310-210-1100101max10!5 5@5105(TotalTime)Maximum Total Computation Time(a)0 0.2 0.4 0.6 0.8 1n10-410-310-210-1100101max10!5 5@5105(ODETime)Maximum ODE Computation TimeN=100N=300N=500N=800N=1000(b)Figure 6.2: Maximum total computation time (a) and ODE computation time (b) mea-sured in seconds, where the maximum is taken over 10−5 ≤ χ ≤ 105.0 0.2 0.4 0.6 0.8 1n01234max10!5 5@5105(Iterations)Maximum IterationsN=100N=300N=500N=800N=1000(a)0 0.2 0.4 0.6 0.8 1n10-710-610-510-410-3max10!5 5@5105(Residual)Maximum ResidualN=100N=300N=500N=800N=1000(b)Figure 6.3: Maximum number of iterations (a) and final residual (b) for Newton’s method,where the maximum is taken over 10−5 ≤ χ ≤ 105.48In Figure 6.2 we have plotted the maximum computation time on the left-handside, while on the right side we have included the component of that time dedicated tosolving the constant δ˜ ODE using ode45. This figure reveals that the total computationtime appears to remain constant as n is varied but, as expected, increases with increasingvalues of N . The same trend is present for the ODE computation time, except when n = 0for which we have an exact solution to the constant δ˜ ODE given by (5.24).Next, Figure 6.3 shows how the maximum number of iterations and maximumresidual vary with n. The left pane of the figure reveals that the program typically haltsafter two iterations regardless of the values of N and n being used. Moreover, as indicatedby the right pane of the figure, the desired tolerance is typically surpassed. These twoobservations suggest that the initial guess, i.e. that obtained by solving the approximatingODE with ode45, is already very close to the actual solution. We will be exploring thissuggestion further in the final section where we compare solutions to the approximatingODE with that obtained by solving the full non-singular integral equation.6.2 Qualitative Properties of the Numerical SolutionThis brief sections showcases some of the quantitative properties of solutions to the non-singular integral equation (4.8). In the left panes of Figures 6.4 to 6.6 we have collected acolor map showing how the fracture opening w˜ varies with both x˜ and χ. The right paneof each figure shows some sample solutions o w˜ versus x˜ for values of χ = 10−5, 103, 105.We have omitted the n = 0 case since then there is no dependence on χ which wouldotherwise trivialize the figure. In each simulation we have used the following numericalparameters (see Appendix B.1 for definitions)% numeric solution parametersp0 = -10; p1 = 15; N = 500; ittmax = 500; tol = 1e-3;49100~x100102104106~w@ = 10!5@ = 103@ = 105Figure 6.4: Numerical results for n = 1.In the left pane of each of Figures 6.4 to 6.6 we have also include dashed yellow(left), green (center), and red (right) lines which indicate the boundaries of the toughness,leak-off, and viscosity regimes respectively. Within the region left of the dashed yellowline the relative difference between the numerical solution w˜ and the toughness solutionw˜k = 1 is less than 10−2 (i.e. |w˜ − w˜k|/w˜k < 10−2). Inside the region bounded by thedashed green line the relative difference between the numerical solution and the w˜ and theleak-off solution w˜c = βcmχn2(1+n) x˜2−n2(1+n) is less than 10−2 (i.e. |w˜ − w˜c|/w˜c < 10−2). Onthe other hand in the region right of the dashed red line the relative difference between thenumerical solution w˜ and the viscous solution w˜m = βmx˜(2−n)/(2+n) is less than 10−2 (i.e.|w˜− w˜m|/w˜m < 10−2). The right pane of each figure plots sample solutions (i.e. w˜ versusx˜) for three distinct χ values. From these figures we can already draw some conclusions asto the behaviour of the fracture opening as χ and n are varied. For fixed n we see that asχ increases the (leak-off) region situated between the toughness and viscosity boundariesgrows. However the growth of this intermediate region is stunted as n decreases whichcan be seen by comparing the left panes of each figure. This is also supported by the50100~x1001051010~w@ = 10!5@ = 103@ = 105Figure 6.5: Numerical results for n = 0.60.100~x10010510101015~w@ = 10!5@ = 103@ = 105Figure 6.6: Numerical results for n = 0.20.5110-10 100 1010~x10010510101015~w@ = 10-5n=0n=0.2n=0.4n=0.6n=0.8n=110-10 100 1010~x10010510101015~w@ = 10310-10 100 1010~x10010510101015~w@ = 105Figure 6.7: Plots of solutions w˜(x˜) for various n values (see legend) and χ = 10−5, 103, 105.right pane of each figure which shows a decreasing dependence on χ for lower values ofn. This effect culminates when n = 0, which, as already mentioned, we have omitted inthe present plots since this case exhibits no χ dependence (as is clear from the governingequation (4.8)).Using the same numerical parameters as above, we now explore how the behaviourof the fracture opening changes as n is varied. For this comparison it is important torecall the scaling (4.7) in which the length scale l depends on the parameter n. Figure6.7 shows the solution w˜ for various values of n and for χ = 10−5, 103, 105. As expectedwe see that the n = 0 solution is independent of χ and the dependence on χ becomesmore pronounced as n increases. From these plots we also see that for x˜ approximatelysmaller than 1 the fracture opening is larger with greater values of n, but that the oppositeis true for values of x˜ approximately greater than 1. This behaviour is best understoodby recalling the asymptotic propagation regimes introduced in §3 and in particular the52change of the power α with n depicted in Figure 3.1. There we immediately see that,apart from αk, higher powers of n result in smaller values of α and therefore a smallerslope on the plot of log w˜ versus log x˜.6.3 Comparison to Approximating ODEWe conclude this chapter by comparing the accuracy between the solution obtainedusing ode45, the zero δ˜-corrected solution, and the first δ˜-corrected solution. In each casewe will be making comparisons to the solution obtained by solving the full non-singularintegral equation (4.8) numerically as discussed in the first section of this chapter. To re-duce errors due to the discretization of the integral equation, we will be using a large meshpoint count of N = 5, 000. Notice that since we are using Simpson’s rule, and referringto the code found in Appendix B.1, this corresponds to 2N + 1 = 10, 001 mesh points.In addition, to reduce the errors from truncating the infinite integral appearing in (4.8)we will be taking mesh points in the range 10−10 to 1020, but only drawing comparisonsin the truncated region of 10−10 to 1015. The other parameters in the numerical solverthat we are using include a maximum number of iterations of ittmax = 500 and haltingtolerance of tol = 1e-10.For our first comparison we consider the case of n = 0. In Figure 6.8 we show howthe relative error between the full numerical solution and each of the three approximationsvaries with x˜. In this figure we remark that the ode45 and the zero δ˜-corrected solution arenearly identical, as is expected since they are both solutions to the same ODE. Moreover,the maximum relative error is 0.062 for the ode45 and zero δ˜-corrected solutions, whereas itis 0.048 for the first δ˜-corrected solution. Next, in Figure 6.9 we have included the relativedifference between these solutions with the regime solutions w˜k = 1 and w˜m =√4pix˜.These figures reveal that in the appropriate spatial regimes each of the approximating5310-10 10-5 100 105 1010 101510-810-610-410-2100ode45zero-deltafirst-deltaFigure 6.8: Relative error between full numerical solution and solutions to the approxi-mating ODE for n = 0. For ode45 and zero δ˜-corrected solution the maximum relativeerror is ≈ 0.062, while for the first δ˜-corrected solution it is 0.048.ODE solutions fully capture the asymptotic solutions w˜k and w˜m.Next for each value of n = 0.2, 0.4, 0.6, 0.8, 1 we consider the relative error forχ = 10−5, 100, 105. These results are shown in Figures 6.10 – 6.14 for n = 0.2, ..., 1.0respectively. In each such figure we observe that the first δ˜-corrected solution providesa significant advantage over the ode45 and zero δ˜-corrected solutions. Perhaps one ofthe most significant conclusions we can draw from these figures is that for n ≥ 0.4 themaximum relative error of the first δ˜-corrected solution is at most 10−2. Computationally,this observation is very important because it means we can effectively use the significantlyless computationally intensive first δ˜-corrected solution rather than the full numericalsolution of the non-singular integral equation. For values of n < 0.4 we are no longerin this desirable position and it is likely better not to fully rely on the first δ˜-correctedsolution unless some other means of improving its accuracy are found.Our final commentary on the approximating ODE solutions revolves around the“bump” features present in Figure 6.8 and Figures 6.10 – 6.14. Motivated by the derivationof the approximating ODE, we claim that these bumps correspond to transition regions5410-10 10-5 100 105 1010 101510-1010-5100Rel. Diff. with Toughness Solutionode45zero-deltafirst-delta10-10 10-5 100 105 1010 101510-1010-5100Rel. Diff. with Viscous Solutionode45zero-deltafirst-deltaFigure 6.9: Relative difference of approximating ODE solutions with regime solutions(toughness and viscosity).55100~x10-810-610-410-2100@ = 10-5100~x10-810-610-410-2100@ = 100100~x10-810-610-410-2100@ = 105ode45zero-deltafirst-deltaFigure 6.10: Relative error between full numerical solution and solutions to the approxi-mating ODE for n = 0.2.between each of the propagation regimes (toughness, leak-off, and viscosity) as this iswhere δ˜ undergoes a transition. Indeed, for n = 0 by comparing Figures 6.8 and 6.9 wesee that the bump in the former occurs in the region where the latter indicates a deviationfrom the toughness solution but the beginning of convergence towards the viscous solution.To support this claim for n > 0 we have included Figures 6.15 – 6.17 which show therelative error for values of χ = 105, 100, 10−5 respectively. Each figure includes three panescorresponding to values of n = 0.2, 0.5, 0.8. We remark that we have ignored the solutionobtained using ode45 so as not to clutter the plots. In addition, each figure includestwo or three dashed vertical lines which indicate the boundaries of the toughness (black),leak-off (green), and viscous (red) regimes defined using the same criteria as was used inFigures 6.4 – 6.6. In particular, the region to the left of the black dashed line correspondsto the toughness regime, the region in between the yellow dashed lines to the leak-off56regime, and the region to the right of the red dashed line to the viscous regime. Thesefigures support our claim that peaks occur at the transitions between different propagationregimes. Moreover, the relative severity of bumps can be explained by referring to Figure3.1, since in general the transition from the toughness to the leak-off region correspondsto a larger jump in the power α than does the transition from the leak-off to the viscosityregime (i.e. αcm − αk > αm − αcm).57100~x10-810-610-410-2100@ = 10-5100~x10-810-610-410-2100@ = 100100~x10-810-610-410-2100@ = 105ode45zero-deltafirst-deltaFigure 6.11: Relative error between full numerical solution and solutions to the approxi-mating ODE for n = 0.4.100~x10-810-610-410-2100@ = 10-5100~x10-810-610-410-2100@ = 100100~x10-810-610-410-2100@ = 105ode45zero-deltafirst-deltaFigure 6.12: Relative error between full numerical solution and solutions to the approxi-mating ODE for n = 0.6.58100~x10-810-610-410-2100@ = 10-5100~x10-810-610-410-2100@ = 100100~x10-810-610-410-2100@ = 105ode45zero-deltafirst-deltaFigure 6.13: Relative error between full numerical solution and solutions to the approxi-mating ODE for n = 0.8.100~x10-810-610-410-2100@ = 10-5100~x10-810-610-410-2100@ = 100100~x10-810-610-410-2100@ = 105ode45zero-deltafirst-deltaFigure 6.14: Relative error between full numerical solution and solutions to the approxi-mating ODE for n = 1.59100~x10-1010-810-610-410-2100n=0.2100~x10-1010-810-610-410-2100n=0.5100~x10-1010-810-610-410-2100n=0.8zero-deltafirst-deltaFigure 6.15: Relative error of zero- and first-delta corrected solution for χ = 105 withdashed vertical lines showing regime boundaries.100~x10-1010-810-610-410-2100n=0.2100~x10-1010-810-610-410-2100n=0.5100~x10-1010-810-610-410-2100n=0.8zero-deltafirst-deltaFigure 6.16: Relative error of zero- and first-delta corrected solution for χ = 100 withdashed vertical lines showing regime boundaries.60100~x10-1010-810-610-410-2100n=0.2100~x10-1010-810-610-410-2100n=0.5100~x10-1010-810-610-410-2100n=0.8zero-deltafirst-deltaFigure 6.17: Relative error of zero- and first-delta corrected solution for χ = 10−5 withdashed vertical lines showing regime boundaries.61Chapter 7ConclusionIn this thesis we have considered the problem of a semi-infinite hydraulic fracture underplane strain conditions driven by a shear-thinning power-law fluid. The motivation forstudying the semi-infinite fracture arises from its close connection to the tip behaviourof a finite fracture [19] and its relevance to their numerical simulation. Furthermore, wehave focused on shear-thinning power-law fluids as this is expected to be the dominantfracturing fluid rheology in oil reservoir stimulation operations [8, 3]. Our model usesPoiseuille’s law for a power-law fluid to describe the fracturing fluid’s motion, but otherwiseneglects the power-law rheology when describing fluid leak-off. This is based on theunderlying assumption that the power-law properties of the fracturing fluid arise from theaddition of filter-cake-building polymers, which themselves do not experience leak-off.The main contribution of this thesis has been to expand on the work of Dontsovand Peirce [6] to accommodate shear-thinning power-law fluids. In particular, in Chapter3 we have included a detailed derivation of the asymptotic solutions that dominate in eachof the three propagation regimes. This was followed in Chapters 4 and 5 by reformulatingthe governing equations to a non-singular integral equation and then reducing it to anapproximating ODE. In solving the approximating ODE we considered three methods,62the last of which involved the use of the “delta-correction” method outlined in [6]. Thismethod involves evaluating two integrals, which is easily accomplished using standardquadrature rules, and therefore potentially represents a computationally fast alternativefor determining the fracture’s tip behaviour.In Chapter 6 we considered a numerical scheme for solving the non-singular integralequation in which the integral was first discretized using Simpson’s rule and the resultingnon-linear algebraic system was solved using Newton’s method. As a first guess in theNewton’s method we used the solution obtained by solving the approximating ODE usingMatlab’s ode45 function. The latter part of Chapter 6 focused on comparing the solutionsto the approximating ODE to those obtained by solving the non-singular integral equation.These comparisons revealed that the solutions to the approximating ODE are in goodagreement with the full numerical solution. In particular, for values of 0 ≤ n < 0.4 themaximum relative error was found to be about 6% whereas for 0.4 ≤ n ≤ 1 the firstdelta-corrected solution yielded a maximum relative error that is less than 1%. Thus,at least for n ≥ 0.4 the first delta-corrected solution could effectively replace the morecomputationally intensive solution obtained by solving the non-singular integral equation.We will conclude by listing some directions for future study. This list will be limitedto those directions most relevant to the work presented above. For a more general overviewand future problems we direct the reader to the recent review article by Detournay [4]. Inno particular order, the proposed directions for future study include:1. A rigorous derivation and error analysis of the approximating ODE, as well as ofthe first delta-corrected solution.2. Finding an improved method of approximating solutions to the non-singular integralequation for values of n < 0.4, capable of obtaining relative errors that are ≤ 1%.633. Exploring the effect of the power-law rheology on fluid lag, and in particular de-termining if, as is the case for Newtonian fluids, the lag region is negligible undercommon field conditions.4. Implementing the tip solutions obtained using methods from this thesis in the caseof a penny-shaped fracture driven by a power-law fluid.64Bibliography[1] J. Adachi, E. Siebrits, A. Peirce, and J. Desroches. Computer simulation of hydraulicfractures. International Journal of Rock Mechanics and Mining Sciences, 44(5):739–757, 2007.[2] J. I. Adachi. Fluid-driven fracture in permeable rock, 2001. Copyright - Databasecopyright ProQuest LLC; ProQuest does not claim copyright in the individual un-derlying works; Last updated - 2016-05-05.[3] J. I. Adachi and E. Detournay. Self-similar solution of a plane-strain fracture drivenby a power-law fluid. International Journal for Numerical and Analytical Methods inGeomechanics, 26(6):579–604, 2002.[4] E. Detournay. Mechanics of hydraulic fractures. Annual Review of Fluid Mechanics,48:311–339, 2016.[5] E. Donaldson, W. Alam, and N. Begum. Hydraulic Fracturing Explained: Evaluation,Implementation, and Challenges. Gulf Drilling. Elsevier Science, 2014.[6] E. V. Dontsov and A. P. Peirce. A non-singular integral equation formulation toanalyse multiscale behaviour in semi-infinite hydraulic fractures. Journal of FluidMechanics, 781:R1 (11 pages), 10 2015.65[7] J. Dundurs. Elastic interaction of dislocations with inhomogeneities. In T. Mura,editor, Mathematical Theory of Dislocations, chapter 4, pages 70–115. The AmericanSociety of Mechanical Engineers, New York, N.Y. 10017, 1969.[8] M. J. Economides, K. G. Nolte, U. Ahmed, and D. Schlumberger. Reservoir stimu-lation, volume 18. Wiley Chichester, 2000.[9] D. Garagash and E. Detournay. The tip region of a fluid-driven fracture in an elasticmedium. Journal of Applied Mechanics, 67:1 (10 pages), 6 1999.[10] D. Garagash, E. Detournay, and J. Adachi. Multiscale tip asymptotics in hydraulicfracture with leak-off. Journal of Fluid Mechanics, 669:260–297, 2 2011.[11] A. Griffith. The phenomena of rupture and flow in solids. Philosophical Transactionsof the Royal Society of London A: Mathematical, Physical and Engineering Sciences,221(582-593):163–198, 1921.[12] D. A. Hills, P. A. Kelly, D. N. Dai, and A. M. Korsunsky. Solution of crack problems,volume 44 of Solid Mechanics and its Applications. Kluwer Academic PublishersGroup, Dordrecht, 1996. The distributed dislocation technique.[13] J. Hirth. Introduction to the mathematical theory of dislocations. In T. Mura, editor,Mathematical Theory of Dislocations, chapter 1, pages 1–24. The American Societyof Mechanical Engineers, New York, N.Y. 10017, 1969.[14] M. F. Kanninen and C. L. Popelar. Advanced fracture mechanics. 1985.[15] A. M. Linkov. Analytical solution of hydraulic fracture problem for a non-newtonianfluid. Journal of Mining Science, 49(1):8–18, 2013.66[16] A. M. Linkov. The particle velocity, speed equation and universal asymptotics for theefficient modelling of hydraulic fractures. J. Appl. Math. Mech., 79(1):54–63, 2015.[17] J. A. Mechanics, F. E. C. N. U. (1969.), .-e. Mura, Toshio, and A. S. of MechanicalEngineers. Applied Mechanics Division. Mathematical theory of dislocations, 1969.“Sponsored by Applied Mechanics Division of the American Society of MechanicalEngineers.”.[18] T. Mura. Method of continuously distributed dislocations. In T. Mura, editor, Math-ematical Theory of Dislocations, chapter 2, pages 25–48. The American Society ofMechanical Engineers, New York, N.Y. 10017, 1969.[19] A. Peirce and E. Detournay. An implicit level set method for modeling hydrauli-cally driven fractures. Computer Methods in Applied Mechanics and Engineering,197(33):2858–2885, 2008.[20] M. Williams. Stress singularities resulting from various boundary conditions. Journalof Applied Mechanics, 19:526–528 (2 pages), 4 1952.67Appendix AClosed Formulae for two IntegralsA.1 Proof of Closed Form Formula for∫∞0G(t)tα dtWe will prove the identityJ(α) ≡∫ ∞0G(t)tαdt =2piα(2− α) tan(piα2), (−1 < α < 1).Let 0 < ε 1 and defineJε(α) ≡(∫ 1−ε0+∫ ∞1+ε)G(t)tαdt = J(α)−∫ 1+ε1−εG(t)tαdt.As ε −→ 0 we find that the rightmost term is O(ε) and thereforeJε(α) −→ J(α) (ε −→ 0).Now we writeJε(α) =(∫ 1−ε0+∫ ∞1+ε)1− t2t1+α[log∣∣∣∣1 + t1− t∣∣∣∣+ 2t1− t2]dt,68ΓRε10Figure A.1: Schematic of the integration contour.and integrate by parts to getJε(α) =− (2− α)t−α + αt2−αα(2− α)[log∣∣∣∣1 + t1− t∣∣∣∣+ 2t1− t2]∣∣∣∣∞0+(2− α)t−α + αt2−αα(2− α)[log∣∣∣∣1 + t1− t∣∣∣∣+ 2t1− t2]∣∣∣∣1+ε1−ε+4α(2− α)(∫ 1−ε0+∫ ∞1+ε)(2− α)t−α + αt2−α(1− t2)2 dtSince 1 < α < 1 andlog∣∣∣∣1 + t1− t∣∣∣∣+ 2t1− t2 =O(t), as t→ 0,O(t−3), as t→∞,the boundary terms at 0 and ∞ vanish. Carrying out the calculation for the boundaryterms at 1− ε and 1 + ε we obtain a single O(ε−1) term leading toJε(α) = − 4α(2− α)1ε+4α(2− α)(∫ 1−ε0+∫ ∞1+ε)(2− α)t−α + αt2−α(1− t2)2 dt. (A.1)To evaluate the remaining integral we calculate the contour integral∫Γ(2− α)z−α + αz2−α(1− z2)2 dz,69where the contour Γ consists of one large circular contour of radius R centred at z = 0and two small circular contours centred at z = 0 and z = 1 as shown in Figure A.1.Furthermore we introduce a branch cut along the positive real axis to account for thepower zα. In the limit of ε −→ 0 the contribution from the small circle at z = 0 vanishesand the remaining contributions equalO(R−α−3)+O(R−1−α) + (1− e−i2piα)(∫ 1−ε0+∫ ∞1+ε)(2− α)t−α + αt2−α(1− t2)2 dt−∫ pi0(2− α)(1 + εeiθ)−α + α(1 + εeiθ)2−α(−2εeiθ − ε2ei2θ)2 εieiθdθ−e−i2piα∫ pi0(2− α)(1 + εe−iθ)−α + α(1 + εe−iθ)2−α(−2εe−iθ − ε2e−i2θ)2 εie−iθdθ,where the O(R−α−3) and O(R−1−α) terms are coming from integrating along the largeouter circle. Meanwhile the latter two terms equal ipi2 − 1ε and ipi2 e−i2piα + 1εe−i2piα respec-tively with an O(ε) correction. On the other hand, the contour integral can be evaluatedusing Cauchy’s residue theorem which yields ipie−ipiα. Combining these results we obtain(∫ 1−ε0+∫ ∞1+ε)(2− α)t−α + αt2−α(1− t2)2 dt =pi2tan(piα2)+1ε+O(ε),whenceJε(α) =2piα(2− α) tan(piα2)+O(ε),and the desired result follows by taking the limit ε→ 0.A.2 Proof of Closed Form Formula for∫∞0G′(t)tκ dtWe will prove the identity∫ ∞0G′(t)tκdt =2piκ(1 + κ)(1− κ) tan(pi2(κ+ 1))(−2 < κ ≤ 0).70First, we consider the case −2 < κ < 0. Integrating by parts and using G(t) ∼ O(1) ast→ 0 and G(t) ∼ O(t−2) as t→∞ to cancel the boundary terms we obtain∫ ∞0G′(t)tκdt = κ∫ ∞0G(t)tκ+1dt.Since −1 < κ+ 1 < 1 we can use the identity (4.9) from which the desired result follows.When κ = 0 we obtain∫ ∞0G′(t)dt = G(∞)−G(0) = −4,which agrees with the desired identity upon setting κ = 0.71Appendix BMatlab CodeB.1 Newton’s Method Solver for Non-Singular Integral Equa-tionfunction [xt,wt,wta,itt,Res,totTime,odeTime] =MKC_calc_w(n,chi,p0,p1,N,ittmax,tol)% ----------------------- Preamble ----------------------------------------%% DESCRIPTION:% This function solves for $\tilde{w}(\tilde{x})$ by first discretizing% the non-linear integral equation using Simpson’s rule and then solving% the resulting non-linear system of equations using Newton’s method. The% initial guess used for Newton’s method comes from solving the linear% approximating ODE using ode45.%% INPUTS:% n - flow behaviour index% chi - leakoff parameter% p0,p1 - powers determining range of xt via 10^p0 < xt < 10^p1% N - number of intervals for Simpsons rule (note length(xt)=2*N+1)% ittmax - maximum number of iterations in Newton’s method% tol - desired relative tolerance in Newton’s method%% OUTPUTS:% xt - \tilde{x} values (logarithmically distributed)% wt - the solution $\tilde{w}$% wta - solution of the approximating ODE% itt - number of iterations before Newton’s method termianted72% Res - maximum relative change during last iteration% totTime - total computation time% ------------------------ Initialization ---------------------------------% initiate the Matlab timertic% generate logarithmically distributed nodes for $\tilde{xt}$xt = logspace(p0,p1,2*N+1);% beta_m for the viscous regimeif n == 0beta_m = 2*sqrt(pi);elsebeta_m = (-2*((n+2)^2)/(n*cot(2*pi/(n+2))))^(1/(n+2));end% -------------------- Infinity Residual Computation ----------------------% initialize infinity residual vectorI_inf = zeros(length(xt),1);% infinity residual for X_i < X_{2N+1}I_inf(1:1:end-1) = (8/pi)*(4*(2+n)/(3*(beta_m^(n+1))*(2+3*n)))...*(xt(1:1:end-1)’.^2)./(xt(end)^((2+3*n)/(2+n)));% coefficient calculation for X_i = X_{2N+1}M = 20;n_par = 2*n/(2+n);Coeff_Sum = 1/(3*(1+n_par));for k=1:1:MCoeff_Sum = Coeff_Sum + 1/((2*k+1*n_par)*(2*k+1)*(2*k+3));end% infinity residual for X_i = X_{2N+1}I_inf(end) = (8/pi)*(4*Coeff_Sum/(beta_m^(n+1)))*(xt(end)^((2-n)/(2+n)));% ---------- Initial Approximation for Newton’s Method -------------------% solve approximating linear ODEif (chi == 0) || (n == 0)% exact solution available for n = 0 and chi = 0wta = (1 + (beta_m^(2+n))*xt.^(2-n)).^(1./(2+n))’;else% solve numerically using ode45% calculate ODE coefficients73g1 = -2*(4-n^2)*tan(2*pi/(2+n))/n;g2 = -32*(2-n)*(1+n)*tan(pi*(4+n)/(4*(1+n)))/(3*n*(4+n));% solve the constant delta approximating ODE using ode45odefun = @(x,w) ((1+chi./w).^(n-1)).*(g1./w.^(1+n) +chi*g2./w.^(2+n)).*x.^(1-n);[~,wta] = ode45(@(X,W)odefun(X,W),[0,xt],1);wta = wta(2:end);endodeTime = toc;% ---------- Simpson’s Rule Set Up ---------------------------------------% ds^(2-n)/(2-n) for Simpson’s ruleds = zeros(2*N+1,1);ds(1) = xt(3)^(2-n) - xt(1)^(2-n);ds(2:2:2*N) = 4*(xt(3:2:2*N+1).^(2-n)-xt(1:2:2*N-1).^(2-n));ds(3:2:2*N-1) = xt(5:2:2*N+1).^(2-n)-xt(1:2:2*N-3).^(2-n);ds(2*N+1) = xt(2*N+1)^(2-n) - xt(2*N-1)^(2-n);ds = ds/(6*(2-n));% matrix with G values[ST,XT] = meshgrid(xt,xt);G = func_G(ST./XT);% matrix representing the discretized kernel as a linear operator[dS,~] = meshgrid(ds,ds);B = (8/pi)*G.*dS;% ---------- Newton’s Method Solver --------------------------------------% initialize Newton’s methodRes = 1;itt = 0;wt = wta;% main loop for Newton’s methodwhile (Res>tol)&&(itt<ittmax)% incremenet iteration counteritt=itt+1;% check if maximum allowable iterations have been reachedif itt==ittmax-1disp(’No convergence, MKC2’);disp(Res);end74% non-linear algebraic equation (solving F = 0)F = 1 - wt + B*((1./wt.^(n+1)).*((1+chi./wt).^n)) + I_inf;% calculate Jacobian of F for Newton’s methodJ = -1*(n+1)*B*diag((1./wt.^(n+2)).*((1+chi./wt).^n)) ...- chi*n*B*diag((1./wt.^(n+3)).*((1+chi./wt).^(n-1))) ...- eye(length(wt));% calculate next iterated solutionwt_temp = wt - J\F;Res = max(abs(wt_temp-wt)./wt_temp);wt = wt_temp;end% total elapsed timetotTime = toc;endB.2 Zero Delta-Corrected Solutionfunction [ xt0 ] = fcnV_xt0(wt,chi,n)% vectorized xt0 function% make sure wt is a column vectorN = length(wt);wt = reshape(wt,N,1);if n ~= 0% n != 0 case% load Gaussian Quadrature points and weightsgauss_quad = get_gauss();quad_weights = gauss_quad(:,1);quad_points = gauss_quad(:,2);M = length(quad_weights);% constants for zero-order solutiong1 = -2*(4-n^2)*tan(2*pi/(2+n))/n;g2 = -32*(2-n)*(1+n)*tan(pi*(4+n)/(4*(1+n)))/(3*n*(4+n));b0 = g2/g1;% change of variables coefficients75A = 0.5*(wt-1); B = 0.5*(wt+1);% change of variables evaluated at the quadrature pointsV = A*quad_points’+B*ones(1,M);% integrand in matrix formJ0 = (V+chi).^(1-n).*V.^(1+2*n)./(V+b0*chi);% zero-order solution (n ~= 0)xt0 = (((2-n)/g1)*A.*(J0*quad_weights)).^(1./(2-n));else% zero-order solution (n == 0)xt0 = 0.5*((wt.^2-1)./pi).^0.5;endendB.3 First Delta-Corrected Solutionfunction [ xt1 ] = fcnV_xt1(wt,chi,n)% function [ xt2, C1, C2 ] = fcnV_xt1(wt,chi,n)% vectorized xt1 function% make sure wt is a column vectorN = length(wt);wt = reshape(wt,N,1);if n ~= 0% n != 0 case% load Gaussian Quadrature points and weightsgauss_quad = get_gauss();quad_weights = gauss_quad(:,1);quad_points = gauss_quad(:,2);M = length(quad_weights);% constants for zero-order solutiong1 = -2*(4-n^2)*tan(2*pi/(2+n))/n;g2 = -32*(2-n)*(1+n)*tan(pi*(4+n)/(4*(1+n)))/(3*n*(4+n));b0 = g2/g1;% first delta-correctiond1 = (g1./(1+chi./wt).^(1-n)).*(fcnV_xt0(wt,chi,n).^(2-n)./...76wt.^(2+n)).*(1+chi.*b0./wt);% constants for first-delta iterationC1 = 16*(2-n-(n+1).*d1).*tan(0.5*pi*(1-n-(n+1).*d1))./((1-n-(n+1).*d1)....*(3-n-(n+1).*d1));C2 = 16*(2-n-(1+2*n).*d1).*tan(0.5*pi*(1-n-(1+2*n).*d1))./((1-n-(1+2*n)....*d1).*(3-n-(1+2*n).*d1));% change of variables coefficientsA = 0.5*(wt-1); B = 0.5*(wt+1);% change of variables evaluated at the quadrature pointsV = A*quad_points’+B*ones(1,M);% integrand in matrix formJ1 = (V+chi).^(1-n).*V.^(1+2*n)./...(V.*(C1*ones(1,M))+chi*C2*ones(1,M));% first-delta-corrected solution (n ~= 0)xt1 = ((2-n)*A.*(J1*quad_weights)).^(1./(2-n));else% n = 0 case% first delta-correctiond1 = (wt.^2-1)./wt.^2;% constant for first-delta correctionC0 = 16.*(2-d1).*tan(0.5*pi*(1-d1))./((1-d1).*(3-d1));C0(abs(1-d1)<1e-6) = 4*pi*ones(size(C0(abs(1-d1)<1e-6)));% first-delta-corrected solution (n == 0)xt1 = ((wt.^2-1)./C0).^0.5;end77


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items