Skip to main content

Hydraulic fracturing and its peculiarities

Abstract

Background

Simulation of pressure-induced fracture in two-dimensional (2D) and three-dimensional (3D) fully saturated porous media is presented together with some peculiar features.

Methods

A cohesive fracture model is adopted together with a discrete crack and without predetermined fracture path. The fracture is filled with interface elements which in the 2D case are quadrangular and triangular elements and in the 3D case are either tetrahedral or wedge elements. The Rankine criterion is used for fracture nucleation and advancement. In a 2D setting the fracture follows directly the direction normal to the maximum principal stress while in the 3D case the fracture follows the face of the element around the fracture tip closest to the normal direction of the maximum principal stress at the tip. The procedure requires continuous updating of the mesh around the crack tip to take into account the evolving geometry. The updated mesh is obtained by means of an efficient mesh generator based on Delaunay tessellation. The governing equations are written in the framework of porous media mechanics and are solved numerically in a fully coupled manner.

Results

Numerical examples dealing with well injection (constant inflow) in a geological setting and hydraulic fracture in 2D and 3D concrete dams (increasing pressure) conclude the paper. A counter-example involving thermomecanically driven fracture, also a coupled problem, is included as well.

Conclusions

The examples highlight some peculiar features of hydraulic fracture propagation. In particular the adopted method is able to capture the hints of Self-Organized Criticality featured by hydraulic fracturing.

Background

Fluid-driven fracture propagating in porous media is widely used in geomechanics to improve the permeability of reservoirs in oil and gas recovery or of geothermal wells. Another application of importance is related to the overtopping stability analysis of dams. In the case of reservoir engineering, water is forced under high pressure deep into the ground by injection into a well. The fluid, usually mixed with sand and some chemicals, penetrates in the reservoir rock, opening long cracks (fracking). Horizontal drilling together with hydraulic fracturing makes the extraction of tightly bound natural gas from shale formations economically feasible [1]. In the field, it is unfortunately rather difficult to obtain direct information about the evolution of the crack in the ground, and very little data are known or accessible. Two types of measurements are mainly performed: monitoring of pressure fluctuations at the injection pump and registration of acoustic emissions at the surface [2]. Fracking can also induce small earthquakes [3]. Pressure-induced fracture propagation presents some peculiar features such as pressure peaks and stepwise advancement, which have been discovered only recently and need further investigation. It is recalled that differently from tensile experiments where the crack surfaces are stress free, in hydraulic fracturing, these surfaces are loaded by a pressure distribution resulting from the invading fluid or gas [2]. Simulation is an extremely useful tool to obtain more insight into the problem. The paper addresses this issue.

Contributions to the mathematical modelling of fluid-driven fractures have been made continuously since the 1960s, beginning with Perkins and Kern [4]. These authors made various simplifying assumptions, for instance, regarding fluid flow, fracture shape and velocity leakage from the fracture. For other analytical solutions in the frame of linear fracture mechanics, assuming the problem to be stationary, see [59]. They suffer the limits of an analytical approach, in particular the inability to represent an evolutionary problem in a domain with a real complexity. An analysis of solid and fluid behaviour near the crack tip can be found in [10, 11]. Boone and Ingraffea [12] present a numerical model in the context of linear fracture mechanics which allows for fluid leakage in the medium surrounding the fracture and assumes a moving crack depending on the applied loads and material properties. Tzschichholz and Herrmann [2] used a two-dimensional (2D) lattice model for constant injection rate and homogeneous and heterogeneous material which only breaks under tension. Carter et al. [13] show a fully three-dimensional (3D) hydraulic fracture model which neglects the fluid continuity equation in the medium surrounding the fracture. A discrete fracture approach with remeshing in an unstructured mesh and automatic mesh refinement is used by Schrefler et al. [14]. An element threshold number (number of elements over the cohesive zone) was identified to obtain mesh-independent results. This approach has been extended to 3D situation in [15]. Extended finite elements (XFEM) have been applied to hydraulic fracturing in a partially saturated porous medium by Réthoré et al. [16] in a 2D setting. In this case, a two-scale model has been developed for the fluid flow: in the cohesive crack, Darcy's equation is used for flow in a porous medium, and identities are derived that couple the local momentum and mass balances to the governing equations for the unsaturated medium at macroscopic level. As an example, the rupture of a saturated square plate (0.25 × 0.25 m) in plane strain conditions is investigated under a prescribed fixed vertical velocity v = 2.35 × 10-2 μm/s in the opposite direction at the top and bottom of the plate (tensile loading). The mesh used consists of 20 × 20 quadrilateral elements (12.5 × 12.5 mm each) with bilinear shape functions, and the time step size is 1 s. In the cracked region, the elements are further divided in four triangles. Mohammadnejad and Khoei [17] solve the same problem also with XFEM, using full two-phase flow throughout the region. Darcy flow is assumed within the crack. Finer meshes are used as above (smallest element size 4.5 × 4.5 mm) and much lower time steps (0.25 to 0.125 s). Cavitation is found in both papers, also due to the impervious boundary conditions chosen. Partition of unity finite elements (PUFEM) are used for 2D mode I crack propagation in saturated ionized porous media by Kraaijeveldt et al. [18]. A pull test, a delamination test and an osmopolarity test are simulated with rather fine regular meshes (quadrangular elements with side length of the order of 2 mm and lower) and time step size down to 0.1 s. The time and space discretizations, including the element threshold number used for the solutions, are extremely important for catching the phenomena described next.

We address now the peculiar behaviour of hydraulic fracture propagation which has been observed only by a minority of the above-mentioned authors, but has been confirmed experimentally. Tzschichholz and Herrmann [2] have evidenced with their lattice model and constant injection rate a drop in pressure in time and oscillations on short time scales. These authors explain this by the fact that at the beginning high pressures are needed to push the fluid into the crack. The crack is enlarged and the pressure drops because the enlarged crack can now be opened much more easily than before. The pressure goes down although additional fluid has been added to the crack in the time step. If the pressure drops too much, the stresses at the crack tip fall below their cohesion value and the crack cannot grow at the next time step. By injecting more fluid into the crack, the pressure increases linearly in time until the cohesion forces can be overcome again. Using arguments from continuum mechanics, the authors show that the obtained value for pressure decline in the long term agrees acceptably with their numerical results. The short-term deviations are due the lattice model and the ensuing pressure drops. Oscillations are also obtained for the stored lattice energy. The breaking process is discontinuous in time with time intervals of quiescence where all beams on the crack surface are stressed below their cohesion thresholds and the acting pressure increases linearly in time. Tzschichholz and Herrmann [2] also find a temporal clustering of the breaking events, calling such a sequence bursts (avalanche behaviour). The bursts are unevenly distributed in time and occur relatively often for small times and become rarer later. There is resemblance between the obtained data and magnitude records of earthquakes or of acoustic emission records from laboratory experiments. We have shown with our porous media mechanics model in a 2D setting [14] that in the case of hydraulic fracturing the fracture advances stepwise. Two types of mesh refinement in space and refinement in time were used, but the stepwise advancement did not disappear. Such steps do not appear in other coupled solutions involving cohesive fracture, as e.g. the thermo-elastic one of [19] where the crack surfaces are stress free. The stepwise advancement and flow jumps were also found by Kraaijeveld [20] with a strong and a weak discontinuity model for flow. In [18], the stepwise advancement in mode I crack propagation is difficult to see because a continuous pressure profile across the crack is used. However, continuous pressure profile only works for sufficiently fine meshes. If the mesh is sufficiently fine, then the discretization can resolve the steep pressure gradients along the crack, but the advantage of PUFEM which allows keeping the mesh pretty rough all over the continuum is lost. Hence, dealing with the stepwise progression of the crack in this mode I model is only possible with a finer mesh than the one used (JM Huyghe, personal communication). This is why the authors state that the physical phenomenon challenges the numerical scheme. In mode II, as shown in [20, 21], this problem does not appear because a discontinuous pressure across the crack is accounted for. There it is not attempted to resolve the steep pressure gradient, but this gradient is reconstructed afterwards, using the Terzaghi analytical solution for pressure diffusion. This two-step procedure allows using a rough mesh and still handling a realistic pressure gradient. Stepwise crack advancement can clearly be observed in the crack length histories of Figure four of [17], while it does not appear in the solution for the same problem in [16]. The cohesive fracture length for this problem is estimated with Barenblatt's expression (see Equation 22) [22] to be about 136 mm. Hence, there are about 10 elements over the cohesive zone in [16] and 30 elements in [17]. The first value is probably below the element threshold number for this type of problem, even with XFEM (see also the large time steps used), while the second one is sufficient even for standard elements. While both use XFEM, the two-step procedure and the large time step size and coarse mesh in [16] hide the problem.

Finally, stepwise advancement and flow are also mentioned in [23], where PUFEM is used for 2D poroelastic media. Their method still suffers from mesh dependence because the crack propagates through one element each time step. Hence, their conclusions are not definite. However, Pizzocolo et al. [24] confirmed stepwise advancement experimentally with a test on a small hydrogel disk. The duration of the pause Δt between steps is found to be inversely related to the hydraulic permeability K according to Δt = Δx2 / KE with E Young's modulus and Δx length of the step. A possible explanation for the stepwise behaviour observed in [20, 21, 24] put forward in [24] is that an incompressible fluid consolidation comes into play which prevents tip advancement until the overpressures due to the last advancement have been dissipated, and the stress has been transferred again to the solid phase. This implies the existence of pressure peaks after each advancement stage. During the period of quiescence, the effective stress is below the breaking threshold. Consolidation as a possible explanation for the stepwise advancement needs further investigation in the case of fluid injection, because for some permeability values the tip pressure goes down to zero as shown below on an example. The existence of periods of quiescence is in line with the findings of [2]. We will show that this phenomenon is not only relevant for small structures, where it has been observed experimentally, but also for large structures such as underground soil masses and dams. In that case, the fracture length is much larger, but the phenomenon is still there and the bursts can be felt at great distances compared to the crack length.

Methods

The first subsection presents the fracture model, the second subsection summarizes the governing equations and their numerical solution by means of the finite element method and the third subsection explains the adopted fracture advancement procedure and the required refinements necessary to obtain mesh-independent results.

The fracture model

We use a discrete crack model for a situation depicted in Figure 1: Ω is the domain, Γe is the boundary of the fully saturated porous material surrounding the crack, Γ′ is the crack boundary and Ω ¯ is the domain inside the crack filled with fluid only. There is fluid exchange between the crack and the surrounding porous medium. The mechanical behaviour of the solid phase at a distance from the process zone is assumed to obey a Green elastic or hyperelastic material behaviour [14].

Figure 1
figure 1

Hydraulic fracture domain and cohesive crack geometry. Definition of the hydraulic fracture domain, reprinted from [15], Copyright (2012), with permission from Springer, and the cohesive crack geometry, reprinted from [14], Copyright (2006), with permission from Elsevier.

For the fracture itself, we use a cohesive fracture model. Between the real fracture apex which appears at macroscopic level and the apex of a fictitious fracture, there is a process zone where cohesive forces act (see Figure 1). Following [22, 25, 26], the cohesive law for mode I crack opening with monotonically increasing opening is

σ = σ 0 1 - δ σ δ σ cr
(1)

σ0 being the maximum cohesive traction (closed crack), δ σ the current relative displacement normal to the crack, δσ cr the maximum opening with exchange of cohesive tractions and Gc = σ0 × δσ cr / 2 the fracture energy. If after some opening δσ 1 < δσ cr, the crack begins to close and tractions obey a linear unloading as

σ = σ 0 1 - δ σ 1 δ σ cr δ σ δ σ 1
(2)

When the crack reopens, Equation 2 is reversed until the opening δσ 1 is recovered; then, tractions obey again Equation 1.

When tangential relative displacements of the sides of a fracture in the process zone cannot be disregarded, mixed mode crack opening takes place. This is often the case of a crack moving along an interface separating two solid components. In fact, whereas the crack path in a homogeneous medium is governed by the principal stress direction, the interface has an orientation that is usually different from the principal stress direction. The mixed cohesive mechanical model involves the simultaneous activation of normal and tangential displacement discontinuity and corresponding tractions. For the pure mode II, the relationship between tangential tractions and displacements is

τ = τ 0 δ σ δ σ cr δ τ δ τ
(3)

τ0 being the maximum tangential stress (closed crack), δ τ the relative displacement parallel to the crack and δσ cr the limiting value opening for stress transmission. The unloading/loading from/to some opening δσ 1 < δσ cr follows the same behaviour as for mode I.

For the mixed mode crack propagation, the interaction between the two cohesive mechanisms is treated as in [27]. By defining an equivalent or effective opening displacement δ and the scalar effective traction t as

δ = β 2 δ τ 2 δ σ 2 , t = β - 2 τ 2 + σ 2
(4)

the resulting cohesive law is

t = t δ β 2 δ τ + δ σ
(5)

β being a suitable material parameter that defines the ratio between the shear and the normal critical components. For more details, see [14].

Governing equations and their discretization in space and time

Taking into account the cohesive forces and the symbols of Figure 1, the linear momentum balance of the mixture, discretized in space with finite elements according to the standard Galerkin procedure [28] is written as

M v ˙ + Ω B T σ - Qp - f 1 - Γ N u T c d Γ = 0
(6)

where c is the cohesive traction on the process zone as defined above.

The fully saturated medium surrounding the fracture has constant absolute permeability, while for the permeability within the crack, the Poiseuille or cubic law is assumed. This permeability does not depend on the rock type or stress history but is defined by crack aperture only. Deviation from the ideal parallel surface conditions causes only an apparent reduction in flow and can be incorporated into the cubic law, which reads as [29]

k ij = 1 f w 2 12
(7)

w being the fracture aperture and f a coefficient in the range 1.04 to 1.65 depending on the solid material. In the following, this parameter will be assumed as constant and equal to 1.0. Incorporating the Poiseuille law into the weak form of water mass balance equation within the crack and discretizing in space by means of the finite element method results in

H ˜ p + S ˜ p ˙ + Γ N p T q w d Γ = 0
(8)

With

H ˜ = Ω ¯ N p T w 2 12 μ w N p d Ω ¯
(9)
S ˜ = Ω ¯ N p T 1 K f N p d Ω ¯
(10)

μ w is the dynamic viscosity and K f the bulk modulus of the fluid. The last term of (8) represents the leakage flux into the surrounding porous medium across the fracture borders and is of paramount importance in hydraulic fracturing techniques. This term can be represented by means of Darcy's law using the medium permeability and pressure gradient generated by the application of water pressure on the fracture lips. No particular simplifying hypotheses are hence necessary for this term. This equation can be directly assembled at the same stage as the mass balance Equation 11 for the saturated medium surrounding the crack, because both have the same structure: only the parameters have to be changed in the appropriate elements depending whether they belong to the fracture or to the surrounding medium.

The discretized mass balance equation for the porous medium surrounding the fracture is

Q T u ˙ + Hp + S p ˙ - f 2 - Γ N p T q w d Γ = 0
(11)

where qw represents the water leakage flux along the fracture toward the surrounding medium of Equation 7. This term is defined along the entire fracture, i.e. the open part and the process zone. It is worth mentioning that the topology of the domains Ω and Ω ¯ changes with the evolution of the fracture. In particular, the fracture path, the position of the process zone and the cohesive forces are unknown and must be regarded as products of the mechanical analysis.

Discretization in time is then performed with time discontinuous Galerkin approximation following [30, 31]. Denoting with I n = t n - , t n + 1 + a typical incremental time step of size Δt = tn + 1 - t n , the weighted residual forms are

I n δ v T M v ˙ + Ku - Qp - f 1 dt + I n δ u T K u ˙ - v dt + + δ u T t n K u n + - u n - dt + δ v T t n M v n + - v n - = 0
(12)
I n δ p T Q T v + Ss + Hp - f 2 dt + I n δ p T S p ˙ - s dt + + δ p T t n S p n + - p n - dt = 0
(13)

with the constraint conditions

u ˙ - v = 0 p ˙ - s = 0
(14)

Subscripts -/+ indicate quantities immediately before and after the generic time station. Field variables and their first time derivatives at time  t  [t n , tn + 1 ] are interpolated by linear time shape functions, and the following discretized equations are obtained

u n = u n - + Δt 2 v n + - v n + 1 - u n + 1 = u n - + Δt 2 v n + + v n + 1 - s n = 1 Δt p n + 1 + 3 p n - 4 p n - s n + 1 = 1 Δt p n + 1 + 3 p n + 2 p n -
(15)
1 2 M - 5 36 Δ t 2 K v n + 1 2 M + 1 36 Δ t 2 K v n + 1 + Δt 3 Q p n + + Δt 6 Q p n + 1 = - Δt 2 K u n - + M v n - + I n N 1 t f 1 dt - 1 2 M - 7 36 Δ t 2 K v n + 1 2 M + 5 36 Δ t 2 K v n + 1 + Δt 3 Q p n + + Δt 3 Q p n + 1 = - Δt 2 K u n - + I n N 2 t f 1 dt Δt 3 Q T v n Δt 6 Q T v n + 1 + 1 2 S + Δt 3 H p n + 1 2 S + Δt 6 H p n + 1 = = S p n - + I n N 1 t f 2 dt Δt 6 Q T v n Δt 3 Q T v n + 1 + - 1 2 S + Δ t H p n + 1 2 S + Δt 3 H p n + 1 = I n N 1 t f 2 dt
(16)

The nodal displacement, velocity and pressure, u n - , v n - and p n - , respectively, for the current step coincide with the unknowns at the end of the previous one, hence are known in the time marching scheme and coincide with the initial condition for the first time step. The system of algebraic equations is solved with a monolithic approach using an optimized non-symmetric sparse matrix algorithm. The number of unknowns is doubled with respect to the traditional trapezoidal method.

In a quasistatic situation, adopted for the examples, the submatrices of the above equations are the usual ones of soil consolidation [28], except for

f ˙ 1 = Ω N u T ρ b ˙ + Γ t N u T t ˙ + Γ N u T c ˙ d Γ
(17)

where c ˙ is the cohesive traction rate and is different from zero only if the element has a side on the lips of the fracture Γ ′. Given that the liquid phase is continuous over the whole domain, leakage flux along the opened fracture lips is accounted for through the H matrix together with the flux along the crack. Finite elements are in fact present along the crack (not shown in Figure 1), which account only for the pressure field and have no mechanical stiffness. In the present formulation, non-linear terms arise through cohesive forces in the process zone and permeability along the fracture.

Fracture advancement and refinement strategy

Because of the continuous variation of the domain as a consequence of the propagation of the cracks, also the boundary Γ′ and the related mechanical conditions change. Along the formed crack edges and in the process zone, boundary conditions are the direct result of the field equations, while the mechanical parameters have to be updated. The following remeshing techniques account for all these changes [15, 32].

For the fracture nucleation and advancement, the Rankine criterion is used. More than one crack can open and fractures can also branch. Fracture forms and advances if the maximum principal stress exceeds in a point the fixed threshold. The fracture advancement procedure differs for 2D and 3D situations: in 2D, the fracture follows directly the direction normal to the maximum principal stress, while in 3D, the fracture follows the face of the element around the fracture tip which is closest to the normal direction of the maximum principal stress. In this last case, the fracture tip becomes a curve in space (front). The advancement of a fracture creates new nodes: in 2D, the resulting new elements for the filler at the front are triangles, while in 3D situation, they are tetrahedral. If an internal node along the process zone advances in a 3D setting, a new wedge element results in the filler [15].

At each time station t n , j successive tip (front) advancements are possible until the Rankine criterion is satisfied (Figure 2). Their number in general depends on the chosen time step increment Δt, the adopted crack length increment Δs and the variation of the applied loads. This requires continuous remeshing with a consequent transfer of nodal vectors from the old to the continuously updated mesh by a suitable operator v m (Ωm + 1) = (v m (Ω m )). For momentum and energy conservation, the solution is repeated with the quantities of mesh m but re-calculated on the new mesh m + 1 before advancing the crack tip [33].

Figure 2
figure 2

Multiple advancing fracture step at the same time station. Reprinted from [15], Copyright (2012), with permission from Springer.

Three types of refinement are needed to obtain satisfactory results: the refinement in space in general, the satisfaction of an element threshold number over the process zone and a refinement in time. For refinement and de-refinement in space, the Zienkiewicz-Zhu error estimator is used [34]. Fluid lag, i.e. negative fluid pressures at the crack tip, may arise in the case of injection if the speed at which the crack tip advances is sufficiently high so that for a given permeability water cannot flow in fast enough to fill the created space. This as well as mesh-independent results can be obtained numerically only if an element threshold number is satisfied over the process zone. This number is given by the ratio of elements over the process zone and its length and can be estimated in advance from the problem at hand and the expected process zone. The number of elements over the process zone is of paramount importance and has not received sufficient attention by many authors. It is a sort of object-oriented refinement and is extensively dealt with in [14]. Adaptivity in time is obtained by means of the adopted discontinous Galerkin method in the time domain (DGT) [33]. The error of the time-integration procedure can be defined through the jump of the solution

u n = u n - u n - v n = v n - v n - p n = p n - p n -
(18)

at each time station, i.e. the difference between the final point of time step n - 1 and the first point of time step n. By adopting the total energy norms as error measure, we define the following terms:

e u n = v n T M v n + u n T K u n 1 / 2 e u , p n = u n T Q p n 1 / 2 e p n = p n T Q T u n + p n T H p n Δt + p n T P T p n 1 / 2 e n = max e u n , e u , p n , e p n
(19)

Error measures defined in Equation 19 account at the same time for the cross effects among the different fields and the ones between space and time discretizations.

The relative error is defined as in [30]

η n = e n e max
(20)

where emax is the maximum total energy norm

e max = max e i , 0 < i < n

When η > ηtoll, the time step Δt n is modified and a new Δ t n < Δ t n is obtained according to

Δ t n = θ η toll η 1 / 3 Δ t n
(21)

where θ < 1.0 is a safety factor. If the error is smaller than a defined value ηtoll,min, the step is increased using a rule similar to Equation 21.

As it stands, the refinements in space and time are carried out sequentially, starting with the space refinement, followed by the element threshold number and then the refinement in time. An eye is kept on the satisfaction of the discrete maximum principle [35] which states that it is not possible to refine in time below a certain limit depending on the material properties without also refining in space. A proper functional would be needed to link all the three refinements. A flow chart of the numerical procedure is given in the ‘Appendix’.

Results and discussion

First, we show for comparison purposes the results of cohesive fracture propagation in a thermo-elastic medium, a coupled problem solved with the method for fracture advancement outlined above [19]. The method itself of the crack tip advancement does not introduce steps, once mesh-independent results are achieved. This was also found in static problems (not shown here) such as the three-point bending test, the four-point shear test and the case of a plate with a circular hole. The problem shown here evidences also the importance of the element threshold number, i.e. the number of elements over the process zone. A three-point bending test is performed on a bimaterial specimen subjected to a thermo-mechanical loading [36]. One part of the sample is made of aluminium 6061 and the other of polymethylmethacrylate (PMMA), bonded with methacrylate adhesive. The geometry is presented in Figure 3: the sample has a notch with a sharp tip of 1-mm width and 30-mm height shifted 3 mm from the interface in the PMMA zone. The two materials present very different Young's moduli and thermal expansion coefficients, so that, when the system is subjected to heat, stresses arise near the interface as a result of the mismatch in thermal expansion.

Figure 3
figure 3

Geometry of the three-point bending test for a bimaterial specimen. Reprinted from [19], Copyright (2004), with permission from Elsevier.

Two different experiments are reported in [36]. In the first, at a room temperature of 25°C, a load was applied 3 mm from the interface in the PMMA zone (Figure 3) to trigger the fracture process. The loading rate was very low and the resulting speed of crack propagation at the initial stages was also quite slow, so that quasistatic conditions can be assumed. The crack path was individuated, and stresses near the crack tip in the PMMA were measured using a shearing interferometer.

In the second experiment, the same operations were performed when the temperature of the aluminium was 60°C in steady state conditions. To reach these conditions, a cartridge heater (Q in Figure 3) was inserted into the aluminium part near the external vertical side. The variation in time of the PMMA temperature was checked before the fracture test, which was performed when steady state conditions were reached. The temperature of PMMA was recorded at the crack tip location, at 5 and 7 mm from the interface. Also in this case, the crack path was spotted. From the differences between the two situations, the authors gathered the thermal effects, which were independent of the magnitude of the applied mechanical load.

In the two experiments, the crack propagation trajectories differ as shown in Figure 4a,b where a zoom of the fractured specimens in correspondence of the notch is presented. In particular, the crack path is closer to the interface when the temperature is higher. The numerical results are shown in Figure 4c. The agreement is remarkable, see also [19].

Figure 4
figure 4

Zoom of the notch of the specimen with crack path trajectories. (a, b) Experimental results (reproduced from [36]). (c) Numerical results: case A, uniform temperature (25°C); case B, thermal load with E, σ0, δσ cr varying with temperature; and case C, thermal load with E = E(25°C), σ0 = σ0(25°C), δσ cr = δσ cr (25°C). Reprinted from [19], Copyright (2004), with permission from Elsevier.

Application of Barenblatt's theory [22] for calculation of characteristic cohesive zone size l yields for PMMA

= πE G c 8 1 - ν 2 σ 0 2 0.75 mm
(22)

Our numerical results (0.8 mm) are in good accordance with this value. From this value, the choice of the crack tip advancement length can be estimated. It should be such that the heuristically determined element threshold number is satisfied (five elements in the thermo-mechanical case). Using linear elements of decreasing size, the value of the force F (Figure 3), corresponding to an applied vertical displacement on the same point, is calculated. Results are summarized in Figure 5. The peak of the external load and the softening branch are mesh independent once the process zone is subdivided into at least five elements with edges of 0.15 mm or smaller. This situation is handled by the mesh generator simply by locating an element source [32] at the crack tip. Its weight may be a priori stated and/or can be a posteriori updated during the adaptive remeshing procedure once the length of the process zone has been determined. What is important here is that the diagrams in this coupled problem are smooth reasonably once mesh-independent results are obtained.

Figure 5
figure 5

External force vs. vertical displacement and mesh size. Reprinted from [19], Copyright (2004), with permission from Elsevier.

The next application deals with a hydraulically driven fracture due to fluid pumped at constant flow rate Q into a well in 2D conditions (plane strain) [14]. Figure 6 shows the geometry of the problem together with the initial finite element discretization. A notch with a sharp tip is present along the symmetry axis of the analyzed area.

Figure 6
figure 6

Problem geometry for water injection benchmark and overall discretization. Reprinted from [14], Copyright (2006), with permission from Elsevier.

The effects of combined spatial/temporal discretizations are clearly seen in Figure 7, where the crack length is drawn versus time for different tip advancements, Δs, and time step increments, Δt. The correct time history (case E) is obtained by simultaneously reducing these two parameters, whereas the reduction of only one discretization parameter leads to errors (about ±20%) even using small tip advancement, if compared to the crack length. Again, the importance of the element threshold number is evident for the choice of Δs (the length of the process zone according to Equation 22 is 0.8 m, and about 30 elements are needed over it). It clearly appears that the crack tip velocity is very mesh sensitive. Hence, the element threshold number must be satisfied to obtain mesh-independent results.

Figure 7
figure 7

Crack length time history for μ w = 1 × 10-9MPa s and Q= 0.0001 m3/s. Reprinted from [14], Copyright (2006), with permission from Elsevier.

A lower number of elements results in wrong crack tip velocity, and the possible development of fluid lag may be missed [14]. Fluid lag corresponds to negative pressure in the process zone and determines hence different body forces (see Equation 6). The distribution of the pressure over the fracture length at time station 10 min is shown in Figure 8 for the following three combinations of dynamic viscosity and injection rate: μ w  = 1 × 10-9 MPa s, Q = 0.0001 m3/s; μ w  = 1 × 10-11 MPa s, Q = 0.0001 m3/s; and μ w  = 1 × 10-9 MPa s, Q = 0.0002 m3/s. The fracture length clearly varies with the chosen data. For the first combination, the pressure at the fracture tip goes almost to zero, while for lower values of μ w , the pressure is almost constant. For high μ w and doubled injection rate, cavitation occurs.

Figure 8
figure 8

Distribution of the fluid pressure over the fracture length. At time station 10 min for the combinations of dynamic viscosity and injection rate: μ w  = 1 × 10-9 MPa s, Q = 0.0001 m3/s (red circles); μ w  = 1 × 10-11 MPa s, Q = 0.0001 m3/s (green diamonds); and μ w  = 1 × 10-9 MPa s, Q = 0.0002 m3/s (white squares).

The third case deals with the benchmark exercise A2 proposed by ICOLD [37]. The benchmark consists in the evaluation of failure conditions as a consequence of overtopping wave acting on a concrete gravity dam. Contrarily to the previous example here, we have increasing pressure. The geometry of the dam is shown in Figure 9 together with boundary conditions and an intermediate discretization. Differently from the original benchmark, the dam concrete foundation is also considered, which has been assumed homogeneous with the dam body. In such a situation, the crack path is unknown. On the contrary, when a rock foundation is present, the crack naturally develops at the interface between the dam and foundation. In Figure 9 also, the influence of the viscosity on the crack direction is evidenced (circle).

Figure 9
figure 9

Problem geometry for ICOLD benchmark and calculated crack positions. Reprinted from [14], Copyright (2006), with permission from Elsevier.

The initial condition is obtained under self-weight and the hydrostatic pressure due to water in the reservoir up to a level of 52 m. From this point, the water level increases until the overtopping level is reached (higher than the dam crest [14]). The increase of water level in the reservoir is specified in days according to the benchmark.

For an intermediate situation, the principal stress contours and the cohesive forces are shown on Figure 10. Also, fluid lag has been obtained for this situation, not shown in the picture (see [14]). The crack mouth opening displacement versus days is depicted in Figure 11 for different values of the crack tip advancement. The smallest value corresponds to the proper element threshold number. Clearly, stepwise advancement can be observed with some clustering of the steps.

Figure 10
figure 10

Zoom near the fracture on maximum principal stress contour and cohesive forces. Reprinted from [14], Copyright (2006), with permission from Elsevier.

Figure 11
figure 11

Crack mouth opening displacement versus time (days) for different values of the crack tip advancements (mm).

The effects of the stepwise advancement can also be felt at great distance from the actual crack: the horizontal displacements on the dam crest are effected, as can be seen from Figure 12. Only the diagram for the purely elastic solution (no crack) is smooth. Note that here the vertical scale is logarithmic and in the abscissa appear the time steps, not the actual time. This is the reason why the diagram for the elastic case is above the others.

Figure 12
figure 12

Horizontal displacements versus time step of the dam crest. For different values of the crack tip advancements (mm) and without fracture (elastic).

For a similar problem, a 3D solution has been obtained in [15]. In Figure 13, the mesh, the fracture, the process zone and the stress contours are shown when the fracture length is about 15 m corresponding to an intermediate step of the analysis when the water level is 80 m. The horizontal displacement of the dam crest is drawn versus time in Figure 14. The following situations are considered: no fracture at all (elastic); dry fracture (fracture), i.e. water pressure acts only on the dam, not on the crack lips; hydrostatic water pressure in the crack, constant over the crack length (hydraulic fracture); and fully coupled solution with water pressure varying over the crack length (u- p). The last one has fluid exchange between the crack and surroundings. The results for the last case correspond to an intermediate value between the others because the pressure is diminishing towards the crack tip, reaching even negative values there (cavitation).

Figure 13
figure 13

Mesh, fracture, process zone and stress contours. With a fracture length of about 15 m corresponding to an intermediate step of the analysis when the water level is 80 m. Reprinted from [15], Copyright (2012), with permission from Springer.

Figure 14
figure 14

Horizontal displacement of the dam crest versus time. No fracture at all (elastic); dry fracture (fracture), i.e. water pressure acts only on the dam, not on the crack lips; hydrostatic water pressure in the crack, constant over the crack length (hydraulic fracture); and fully coupled solution with water pressure varying over the crack length and fluid exchange between the crack and surroundings (u- p).

The relative variations of the horizontal crest displacements according to

u = u i u el - 1 100
(23)

with u i referring to the studied cases and uel to the elastic solution, are drawn in Figure 15. The largest steps correspond to the situations where fluid is present in the crack and may have pressure exchange (consolidation) with the material surrounding the process zone. Note that these variations are felt on the dam crest, while the pressure-induced fracturing happens on the bottom of the dam. The 3D results have however only qualitative value because the element threshold number would require finer meshes over the cohesive zone which makes the solution very expensive (elements of about 300 mm minimum side length and time steps of a mean value of 4 days were used).

Figure 15
figure 15

Relative displacements versus time at the dam crest.

In all three examples of hydraulic fracturing, the fracture site is not easily accessible. However, the fact that the effects of the stepwise advancement can be felt also at distance as shown in two examples would make them possible to be monitored remotely. Field data and experimental evidence on reservoir rocks and large bodies are still missing. Data could possibly come from fracking sites or from some fracking-induced earthquakes [3].

Conclusions

A fully coupled model for pressure-induced cohesive fracture in a saturated porous medium and its solution by the finite element method has been shown. The model is of the discrete crack type and requires continuous updating of the mesh as the crack tip advances. This is achieved with powerful mesh generators. Three types of refinement are necessary to obtain mesh-independent results: a refinement over the domain of the Zienkiewicz-Zhu type, an element threshold number over the process zone and a refinement in time, here with DGT. The results show that in case of pressure induced fracture with pressure exchange and flow between the fracture and the surrounding medium the crack tip advances stepwise. This was found also by few other authors. Smooth diagrams are found on the contrary in a thermo-elastic fracture which is a coupled problem but with stress-free fracture surfaces. From a comparison of results obtained with different methods by other authors, it appears that in some situations a particular adopted method hides the problems discussed in this paper because the required refinements clash with the raison d'être of such methods like XFEM or PUFEM (adoption of rough meshes). This is the reason why ‘the physical phenomenon challenges the numerical scheme’ [18] and why several authors dealing with hydraulic fracturing have not noticed the peculiar behaviour shown here. Also, two-step procedures may introduce some bias in the solution. Two different explanations are found in the literature for the discussed phenomena: one invokes pressure drop [2] and the other pressure peak [24] after crack advancement. The respective loading conditions are different, but the question deserves further scrutiny. Finally, the stepwise advancement may be relevant for earthquake engineering, see e.g. the resemblance between the obtained data of [2] and magnitude records of earthquakes. In many earthquake-prone regions, there is plenty of water available at the level where the rupture takes place [38, 39]. The problem solved in [17] has been solved again with XFEM and finer mesh in [40] and the steps in the fracture advancement featured in [17] disappeared. This implies that XFEM yields a smooth solution for a phenomenon which in nature is not smooth: as shown in [2] hydraulic fracturing exhibits avalanche behaviour and hints of Self-Organized Criticality.

Appendix

The flow chart of Figure 16 shows the numerical procedure adopted with particular emphasis on the refinement part.

Figure 16
figure 16

Flow chart.

References

  1. Vidic RD, Brantley SE, Vandenbossche JM, Joxtheimer D, Abad JD: Impact of shale gas development on regional water quality. Science 2013, 340: 826–836.

    Article  Google Scholar 

  2. Tzschichholz F, Herrmann HJ: Simulations of pressure fluctuations and acoustic emission in hydraulic fracturing. Phys Rev E 1995, 5: 1961–1970.

    Article  Google Scholar 

  3. Ellsworth WL: Injection induced earthquakes. Science 2013, 341: 142–150.

    Article  Google Scholar 

  4. Perkins TK, Kern LR: Widths of hydraulic fractures. SPE J 1961, 222: 937–949.

    Google Scholar 

  5. Rice JR, Cleary MP: Some basic stress diffusion solutions for fluid saturated elastic porous media with compressible constituents. Rev Geophs Space Phys 1976, 14: 227–241. 10.1029/RG014i002p00227

    Article  Google Scholar 

  6. Cleary MP: Moving singularities in elasto-diffusive solids with applications to fracture propagation. Int J Solids Struct 1978, 14: 81–97. 10.1016/0020-7683(78)90045-8

    Article  Google Scholar 

  7. Huang NC, Russel SG: Hydraulic fracturing of a saturated porous medium—I: general theory. Theor Appl Fract Mech 1985, 4: 201–213. 10.1016/0167-8442(85)90005-9

    Article  Google Scholar 

  8. Huang NC, Russel SG: Hydraulic fracturing of a saturated porous medium—II: special cases. Theor Appl Fract Mech 1985, 4: 215–222. 10.1016/0167-8442(85)90006-0

    Article  Google Scholar 

  9. Detournay E, Cheng AH: Plane strain analysis of a stationary hydraulic fracture in a poroelastic medium. Int J Solids Struct 1991, 27: 1645–1662. 10.1016/0020-7683(91)90067-P

    Article  Google Scholar 

  10. Advani SH, Lee TS, Dean RH, Pak CK, Avasthi JM: Consequences of fluid lag in three-dimensional hydraulic fracture. Int J Num Anal Methods Geomech 1997, 21: 229–240. 10.1002/(SICI)1096-9853(199704)21:4<229::AID-NAG862>3.0.CO;2-V

    Article  Google Scholar 

  11. Garagash D, Detournay E: The tip region of a fluid-driven fracture in an elastic medium. J Appl Mech 2000, 67: 183–192. 10.1115/1.321162

    Article  Google Scholar 

  12. Boone TJ, Ingraffea AR: A numerical procedure for simulation of hydraulically driven fracture propagation in poroelastic media. Int J Num Ana Methods Geomech 1990, 14: 27–47. 10.1002/nag.1610140103

    Article  Google Scholar 

  13. Carter BJ, Desroches J, Ingraffea AR, Wawrzynek PA: Simulating fully 3-D hydraulic fracturing. In Modeling in geomechanics. Edited by: Zaman M, Booker JR, Gioda G. Wiley, Chichester; 2000:525–567.

    Google Scholar 

  14. Schrefler BA, Secchi S, Simoni L: On adaptive refinement techniques in multifield problems including cohesive fracture. Comp Methods Appl Mech Engrg 2006, 195: 444–461. 10.1016/j.cma.2004.10.014

    Article  Google Scholar 

  15. Secchi S, Schrefler BA: A method for 3-D hydraulic fracturing simulation. Int J Fracture 2012, 178: 245–258. 10.1007/s10704-012-9742-y

    Article  Google Scholar 

  16. Réthoré J, de Borst R, Abellan MA: A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks. Comput Mech 2008, 42: 227–238. 10.1007/s00466-007-0178-6

    Article  Google Scholar 

  17. Mohammadnejad T, Khoei AR: Hydromechanical modelling of cohesive crack propagation in multiphase porous media using extended finite element method. Int J Numer Anal Meth Geomech 2013, 37: 1247–1279. 10.1002/nag.2079

    Article  Google Scholar 

  18. Kraaijeveldt F, Huyghe JM, Remmers JJC, de Borst R: 2-D mode one crack propagation in saturated ionized porous media using partition of unity finite elements. J Appl Mech 2013, 80: 020907–1-12.

    Google Scholar 

  19. Secchi S, Simoni L, Schrefler BA: Cohesive fracture growth in a thermoelastic bi-material medium. Comput Struct 2004, 82: 1875–1887. 10.1016/j.compstruc.2004.03.059

    Article  Google Scholar 

  20. Kraaijeveld F: Propagating discontinuities in ionized porous media, Dissertation. Eindhoven University of Technology; 2009.

    Google Scholar 

  21. Kraaijeveldt F, Huyghe JM, Remmers JJC, de Borst R, Baaijens FPT: Shearing in osmoelastic fully saturated media: a mesh-independent model. Engineering Fracture Mechanics; 2014. in press

    Google Scholar 

  22. Barenblatt GI: The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses: axially-symmetric cracks. J Appl Math Mech 1959, 23: 622–636. 10.1016/0021-8928(59)90157-1

    Article  MathSciNet  Google Scholar 

  23. Remij EW, Pizzoccolo F, Remmers JJ, Smeulders D, Huyghe JM: Nucleation and mixed-mode crack propagation in porous material. ASCE, Poromechanics V; 2013:2260–2269. doi:10.1061/9780784412992.247

    Google Scholar 

  24. Pizzocolo F, Hyughe JM, Ito K: Mode I crack propagation in hydrogels is stepwise. Eng Fract Mech 2013, 97: 72–79.

    Article  Google Scholar 

  25. Dugdale DS: Yielding of steel sheets containing slits. J Mech Phys Solids 1960, 8: 100–104. 10.1016/0022-5096(60)90013-2

    Article  Google Scholar 

  26. Hilleborg A, Modeer M, Petersson PE: Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976, 6: 773–782. 10.1016/0008-8846(76)90007-7

    Article  Google Scholar 

  27. Camacho GT, Ortiz M: Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996, 33: 2899–2938. 10.1016/0020-7683(95)00255-3

    Article  Google Scholar 

  28. Lewis RW, Schrefler BA: The finite element method in the static and dynamic deformation and consolidation of porous media. Wiley, Chichester; 1998.

    Google Scholar 

  29. Witherspoon PA, Wang JSY, Iwai KJE, Gale JE: Validity of cubic law for fluid flow in a deformable rock fracture. Water Resour Res 1980, 16: 1016–1024. 10.1029/WR016i006p01016

    Article  Google Scholar 

  30. Li XD, Wiberg NE: Implementation and adaptivity of a space-time finite element method for structural dynamics. Comp Methods Appl Mech Engrg 1998, 156: 211–229. 10.1016/S0045-7825(97)00207-7

    Article  MathSciNet  Google Scholar 

  31. Secchi S, Simoni L, Schrefler BA: Numerical difficulties and computational procedures for thermo-hydro-mechanical coupled problems of saturated porous media. Comput Mech 2008, 43: 179–189. 10.1007/s00466-008-0302-2

    Article  Google Scholar 

  32. Secchi S, Simoni L: An improved procedure for 2-D unstructured Delaunay mesh generation. Adv Eng Softw 2003, 34: 217–234. 10.1016/S0965-9978(02)00131-X

    Article  Google Scholar 

  33. Secchi S, Simoni L, Schrefler BA: Numerical procedure for discrete fracture propagation in porous materials. Int J Num Anal Methods Geomech 2007, 31: 331–345. 10.1002/nag.581

    Article  Google Scholar 

  34. Zhu JZ, Zienkiewic OC: Adaptive techniques in the finite element method. Com Appl Num Methods 1988, 4: 197–204. 10.1002/cnm.1630040210

    Article  Google Scholar 

  35. Rank E, Katz C, Werner H: On the importance of the discrete maximum principle in transient analysis using finite element methods. Int J Num Methods Engng 1983, 19: 1771–1782. 10.1002/nme.1620191205

    Article  MathSciNet  Google Scholar 

  36. Bae JS, Krishnaswamy S: Subinterfacial cracks in bimaterial systems subjected to mechanical and thermal loading. Eng Fract Mech 2001, 68: 1081–1094. 10.1016/S0013-7944(01)00005-4

    Article  Google Scholar 

  37. ICOLD: Fifth international benchmark workshop on numerical analysis of dams. Theme A2, Denver, Colorado. 1999.

    Google Scholar 

  38. Doglioni C, Barba S, Carminati E, Riguzzi F: Fault on-off fluids response. Geoscience Frontiers; 2013:1–14. doi.org/10.1016/j.gsf.2013.08.004

    Google Scholar 

  39. Kelbert A, Schultz A, Egbert G: Global electromagnetic induction constraints on transition-zone water content variations. Nature 2009, 460: 1003–1006. doi:10.1038/nature08257 10.1038/nature08257

    Article  Google Scholar 

  40. Mohammadnejad T, Khoei AR: An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model. Finite Elements in Analysis and Design 2013, 73: 77–95.

    Article  MathSciNet  Google Scholar 

Download references

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Bernhard A Schrefler.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contribution

SS developed the code, devised the crack tip advancement procedure and carried out the simulations BS drafted the manuscript and explained the results in light of new findings in literature. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Secchi, S., Schrefler, B.A. Hydraulic fracturing and its peculiarities. Asia Pac. J. Comput. Engin. 1, 8 (2014). https://doi.org/10.1186/2196-1166-1-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/2196-1166-1-8

Keywords