Access the full text.
Sign up today, get DeepDyve free for 14 days.
Loading next page...
/lp/springer-journals/hydraulic-fracturing-and-its-peculiarities-eTNZvudxnQ
- Publisher
- Springer Journals
- Copyright
- Copyright © 2014 by Secchi and Schrefler; licensee Springer.
- Subject
- Engineering; Theoretical and Applied Mechanics; Computational Science and Engineering; Classical Continuum Physics; Mathematical Applications in the Physical Sciences
- eISSN
- 2196-1166
- DOI
- 10.1186/2196-1166-1-8
- Publisher site
- See Article on Publisher Site
Abstract
bernhard.schrefler@dicea.unipd.it Department of Civil, Environmental Background: Simulation of pressure-induced fracture in two-dimensional (2D) and and Architectural Engineering, three-dimensional (3D) fully saturated porous media is presented together with some University of Padua, 9 via Marzolo, Padua 35123, Italy peculiar features. Full list of author information is Methods: A cohesive fracture model is adopted together with a discrete crack and available at the end of the article 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. Keywords: Hydraulic fracturing; Cohesive model; Fluid lag; Finite elements 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 nat- ural 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 © 2014 Secchi and Schrefler; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 2 of 21 http://www.apjcen.com/content/1/1/8 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 re- cently and need further investigation. It is recalled that differently from tensile experi- ments 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]. Simu- lation 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 au- thors made various simplifying assumptions, for instance, regarding fluid flow, frac- ture 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 [5-9]. 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 homo- geneous 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 investi- −2 μm/s in the opposite gated under a prescribed fixed vertical velocity v =2.35×10 direction at the top and bottom of the plate(tensileloading). The meshusedcon- sists 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 fur- ther divided in four triangles. Mohammadnejad and Khoei [17] solve the same prob- lem 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 satu- rated 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 (quad- rangular elements with side length of the order of 2 mm and lower) and time step Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 3 of 21 http://www.apjcen.com/content/1/1/8 size down to 0.1 s. The time and space discretizations, including the element thresh- old number used for the solutions, are extremely important for catching the phe- nomena 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 con- firmed experimentally. Tzschichholz and Herrmann [2] have evidenced with their lat- tice 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 thetimestep. If thepressuredrops toomuch, thestressesatthe cracktip fall below their cohesion value and the crack cannot grow at the next time step. By inject- ing more fluid into the crack, the pressure increases linearly in time until the cohe- sion 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 ac- ceptably 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 resem- blance 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 step- wise 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 discontinu- ous 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 Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 4 of 21 http://www.apjcen.com/content/1/1/8 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 num- ber 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 con- clusions are not definite. However, Pizzocolo et al. [24] confirmed stepwise advance- ment 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 ac- cording to Δt = Δx / KE with E Young's modulus and Δx length of the step. A pos- sible 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 dissi- pated, 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 quies- cence, 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 struc- tures 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 dis- tances 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, Γ 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 mechan- ical behaviour of the solid phase at a distance from the process zone is assumed to obey a Green elastic or hyperelastic material behaviour [14]. 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 Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 5 of 21 http://www.apjcen.com/content/1/1/8 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. process zone where cohesive forces act (see Figure 1). Following [22,25,26], the cohe- sive law for mode I crack opening with monotonically increasing opening is σ ¼ σ 1− ð1Þ σcr σ being the maximum cohesive traction (closed crack), δ the current relative displace- 0 σ ment normal to the crack, δ the maximum opening with exchange of cohesive trac- σcr tions and G = σ × δ / 2 the fracture energy. If after some opening δ < δ , the c 0 σcr σ1 σcr crack begins to close and tractions obey a linear unloading as δ δ σ1 σ σ ¼ σ 1− ð2Þ δ δ σcr σ1 When the crack reopens, Equation 2 is reversed until the opening δ is recovered; σ1 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 direc- tion. 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 δ δ σ τ τ ¼ τ ð3Þ δ jj δ σcr τ τ being the maximum tangential stress (closed crack), δ the relative displacement par- 0 τ allel to the crack and δ the limiting value opening for stress transmission. The σcr unloading/loading from/to some opening δ < δ follows the same behaviour as for σ1 σcr mode I. Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 6 of 21 http://www.apjcen.com/content/1/1/8 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 dis- placement δ and the scalar effective traction t as qffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 −2 2 2 δ ¼ β δ δ ; t ¼ β τ þ σ ð4Þ τ σ the resulting cohesive law is t ¼ β δ þ δ ð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 momen- tum balance of the mixture, discretized in space with finite elements according to the standard Galerkin procedure [28] is written as Z Z T ″ ðÞ 1 u ′ Mv_ þ B σ dΩ−Qp−f − ðÞ N 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] 1 w k ¼ ð7Þ ij f 12 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 p T w ′ ~ ~ Hp þ S p_ þ ðÞ N q dΓ ¼ 0 ð8Þ With p T p H ¼ ðÞ ∇N ∇N dΩ ð9Þ 12μ p p S ¼ ðÞ N N dΩ ð10Þ μ is the dynamic viscosity and K the bulk modulus of the fluid. The last term of (8) w f represents the leakage flux into the surrounding porous medium across the fracture Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 7 of 21 http://www.apjcen.com/content/1/1/8 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 pres- sure 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 satu- rated 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 T ðÞ 2 p w ′ Q u_ þ Hp þ Sp_ −f − ðÞ N q dΓ ¼ 0 ð11Þ where q 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 pos- ition 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 approxi- − þ mation following [30,31]. Denoting with I ¼ t ; t a typical incremental time step n nþ1 of size Δt = t − t , the weighted residual forms are n +1 n Z Z T ðÞ 1 T δv Mv_ þ Ku−Qp−f dt þ δu KðÞ u_ −v dtþ I I ð12Þ n n T þ − T þ − þδu Ku −u dt þ δv Mv −v ¼ 0 n n n n Z Z T T ðÞ 2 T δp Q v þ Ss þ Hp−f dt þ δp SðÞ p_ −s dtþ I I n n ð13Þ T þ − þδp Sp −p dt ¼ 0 n n n with the constraint conditions u_ −v ¼ 0 ð14Þ p_ −s ¼ 0 Subscripts −/+ indicate quantities immediately before and after the generic time station. Field variables and their first time derivatives at time t ∈ [t , t ] are in- n n +1 terpolated by linear time shape functions, and the following discretized equations are obtained Δt þ − u ¼ u þ v −v n n nþ1 Δt þ − u ¼ u þ v þ v nþ1 n n nþ1 ð15Þ s ¼ p þ 3p −4p n nþ1 n Δt s ¼ p þ 3p þ 2p nþ1 nþ1 n Δt Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 8 of 21 http://www.apjcen.com/content/1/1/8 1 5 1 1 Δt 2 2 M− Δt K v þ M þ Δt K v þ Qp þ n nþ1 2 36 2 36 3 Δt Δt − − ðÞ 1 þ Qp ¼ − Ku þ Mv þ N ðÞ t f dt nþ1 n n 6 2 1 7 1 5 Δt 2 2 − M− Δt K v þ M þ Δt K v þ Qp þ n nþ1 2 36 2 36 3 Δt Δt − ðÞ 1 þ Qp ¼ − Ku þ N ðÞ t f dt nþ1 2 3 2 ð16Þ Δt Δt 1 Δt 1 Δt T T Q v Q v þ S þ H p þ S þ H p ¼ n nþ1 n nþ1 3 6 2 3 2 6 ðÞ 2 ¼ Sp þ N ðÞ t f dt Δt Δt 1 1 Δt T T ðÞ 2 Q v Q v þ S þ ΔtH p þ S þ H p ¼ N ðÞ t f dt n nþ1 1 n nþ1 6 3 2 2 3 − − − The nodal displacement, velocity and pressure, u , v and p , respectively, for the n n n 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 Z Z Z T T T ðÞ 1 u u u ′ _ _ f ¼ ðÞ N ρb dΩ þ ðÞ N t dΓ þ ðÞ N 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 propaga- tion 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 up- dated. 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 ad- vancement procedure differs for 2D and 3D situations: in 2D, the fracture follows dir- ectly 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 Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 9 of 21 http://www.apjcen.com/content/1/1/8 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 set- ting, a new wedge element results in the filler [15]. At each time station t , j successive tip (front) advancements are possible until the Ran- kine 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 (Ω )= ℵ(v (Ω )). m m +1 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]. 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 hi u ¼ u −u n n hi v ¼ v −v ð18Þ n n hi p ¼ p −p n n n Figure 2 Multiple advancing fracture step at the same time station. Reprinted from [15], Copyright (2012), with permission from Springer. Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 10 of 21 http://www.apjcen.com/content/1/1/8 at each time station, i.e. the difference between the final point of time step n −1and the first point of time step n. By adopting the total energy norms as error measure, we define the following terms: 1=2 T T kk e ¼ hi v Mvhi þhi u Kuhi u n n n n 1=2 e ¼ hi u Qphi u;p n ð19Þ 1=2 T T T T T e ¼ hi p Q hi u þhi p HphiΔt þhi p P hi p p n n n n n n no kk e ¼ max kk e ; e ; e u u;p p n n n 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] kk e η ¼ ð20Þ kk e max where ‖e‖ is the maximum total energy norm kk e ¼ max kk e ; 0 < i < n max max i When η > η , the time step Δt is modified and a new Δt < Δt is obtained accord- toll n n ing to 1=3 θη ′ toll Δt ¼ Δt ð21Þ where θ < 1.0 is a safety factor. If the error is smaller than a defined value η , the toll,min 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 re- finement 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 ad- vancement 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 ele- ments 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 inter- face in the PMMA zone. The two materials present very different Young's moduli and Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 11 of 21 http://www.apjcen.com/content/1/1/8 Figure 3 Geometry of the three-point bending test for a bimaterial specimen. Reprinted from [19], Copyright (2004), with permission from Elsevier. 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. 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 trig- ger 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 con- ditions, 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 differ- ences 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]. Application of Barenblatt's theory [22] for calculation of characteristic cohesive zone size l yields for PMMA πEG ℓ ¼ ≅0:75 mm ð22Þ 81ðÞ −ν σ 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 ele- ments 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 Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 12 of 21 http://www.apjcen.com/content/1/1/8 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, σ , δ varying with temperature; and case C, thermal load with E = E(25°C), σ = σ (25°C), δ = δ 0 σcr 0 0 σcr σcr (25°C). Reprinted from [19], Copyright (2004), with permission from Elsevier. the same point, is calculated. Results are summarized in Figure 5. The peak of the ex- ternal 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 deter- mined. What is important here is that the diagrams in this coupled problem are smooth reasonably once mesh-independent results are obtained. 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. 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 simultan- eously 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. A lower number of elements results in wrong crack tip velocity, and the possible de- velopment of fluid lag may be missed [14]. Fluid lag corresponds to negative pressure Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 13 of 21 http://www.apjcen.com/content/1/1/8 Figure 5 External force vs. vertical displacement and mesh size. Reprinted from [19], Copyright (2004), with permission from Elsevier. 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: −9 3 −11 3 μ =1×10 MPa s, Q = 0.0001 m /s; μ =1×10 MPa s, Q = 0.0001 m /s; and μ = w w w −9 3 1×10 MPa s, Q = 0.0002 m /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 μ , the pressure is almost constant. For high μ and doubled w w injection rate, cavitation occurs. 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 overtop- ping 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 Figure 6 Problem geometry for water injection benchmark and overall discretization. Reprinted from [14], Copyright (2006), with permission from Elsevier. Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 14 of 21 http://www.apjcen.com/content/1/1/8 −9 3 Figure 7 Crack length time history for μ =1 ×10 MPa s and Q = 0.0001 m /s. Reprinted from [14], Copyright (2006), with permission from Elsevier. with boundary conditions and an intermediate discretization. Differently from the ori- ginal benchmark, the dam concrete foundation is also considered, which has been as- sumed 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 de- velops at the interface between the dam and foundation. In Figure 9 also, the influence of the viscosity on the crack direction is evidenced (circle). 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 in- creases until the overtopping level is reached (higher than the dam crest [14]). The in- crease 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 corre- sponds to the proper element threshold number. Clearly, stepwise advancement can be observed with some clustering of the steps. 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. 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 frac- ture 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); Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 15 of 21 http://www.apjcen.com/content/1/1/8 Figure 8 Distribution of the fluid pressure over the fracture length. At time station 10 min for the −9 3 combinations of dynamic viscosity and injection rate: μ =1× 10 MPa s, Q = 0.0001 m /s (red circles); −11 3 −9 3 μ =1× 10 MPa s, Q = 0.0001 m /s (green diamonds); and μ =1 × 10 MPa s, Q = 0.0002 m /s w w (white squares). Figure 9 Problem geometry for ICOLD benchmark and calculated crack positions. Reprinted from [14], Copyright (2006), with permission from Elsevier. Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 16 of 21 http://www.apjcen.com/content/1/1/8 Figure 10 Zoom near the fracture on maximum principal stress contour and cohesive forces. Reprinted from [14], Copyright (2006), with permission from Elsevier. 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 re- sults 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 11 Crack mouth opening displacement versus time (days) for different values of the crack tip advancements (mm). Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 17 of 21 http://www.apjcen.com/content/1/1/8 Figure 12 Horizontal displacements versus time step of the dam crest. For different values of the crack tip advancements (mm) and without fracture (elastic). The relative variations of the horizontal crest displacements according to kk u ¼ −1 ⋅100 ð23Þ el with u referring to the studied cases and u to the elastic solution, are drawn in i el 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 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. Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 18 of 21 http://www.apjcen.com/content/1/1/8 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 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 (ele- ments of about 300 mm minimum side length and time steps of a mean value of 4 days were used). In all three examples of hydraulic fracturing, the fracture site is not easily access- ible. 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 re- motely. 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]. Figure 15 Relative displacements versus time at the dam crest. Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 19 of 21 http://www.apjcen.com/content/1/1/8 Figure 16 Flow chart. 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 refine- ment 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 re- quired 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 chal- lenges the numerical scheme’ [18] and why several authors dealing with hydraulic Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 20 of 21 http://www.apjcen.com/content/1/1/8 fracturing have not noticed the peculiar behaviour shown here. Also, two-step proce- dures 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 advance- ment 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 im- plies 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. 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. Author details 1 2 Institute of system science ISIB – CNR, corso Stati Uniti 4, Padua 35127, Italy. Department of Civil, Environmental and Architectural Engineering, University of Padua, 9 via Marzolo, Padua 35123, Italy. Received: 15 November 2013 Accepted: 21 February 2014 Published: 16 May 2014 References 1. Vidic RD, Brantley SE, Vandenbossche JM, Joxtheimer D, Abad JD (2013) Impact of shale gas development on regional water quality. Science 340:826–836 2. Tzschichholz F, Herrmann HJ (1995) Simulations of pressure fluctuations and acoustic emission in hydraulic fracturing. Phys Rev E 5:1961–1970 3. Ellsworth WL (2013) Injection induced earthquakes. Science 341:142–150 4. Perkins TK, Kern LR (1961) Widths of hydraulic fractures. SPE J 222:937–949 5. Rice JR, Cleary MP (1976) Some basic stress diffusion solutions for fluid saturated elastic porous media with compressible constituents. Rev Geophs Space Phys 14:227–241 6. Cleary MP (1978) Moving singularities in elasto-diffusive solids with applications to fracture propagation. Int J Solids Struct 14:81–97 7. Huang NC, Russel SG (1985) Hydraulic fracturing of a saturated porous medium—I: general theory. Theor Appl Fract Mech 4:201–213 8. Huang NC, Russel SG (1985) Hydraulic fracturing of a saturated porous medium—II: special cases. Theor Appl Fract Mech 4:215–222 9. Detournay E, Cheng AH (1991) Plane strain analysis of a stationary hydraulic fracture in a poroelastic medium. Int J Solids Struct 27:1645–1662 10. Advani SH, Lee TS, Dean RH, Pak CK, Avasthi JM (1997) Consequences of fluid lag in three-dimensional hydraulic fracture. Int J Num Anal Methods Geomech 21:229–240 11. Garagash D, Detournay E (2000) The tip region of a fluid-driven fracture in an elastic medium. J Appl Mech 67:183–192 12. Boone TJ, Ingraffea AR (1990) A numerical procedure for simulation of hydraulically driven fracture propagation in poroelastic media. Int J Num Ana Methods Geomech 14:27–47 13. Carter BJ, Desroches J, Ingraffea AR, Wawrzynek PA (2000) Simulating fully 3-D hydraulic fracturing. In: Zaman M, Booker JR, Gioda G (ed) Modeling in geomechanics. Wiley, Chichester, pp 525–567 14. Schrefler BA, Secchi S, Simoni L (2006) On adaptive refinement techniques in multifield problems including cohesive fracture. Comp Methods Appl Mech Engrg 195:444–461 15. Secchi S, Schrefler BA (2012) A method for 3-D hydraulic fracturing simulation. Int J Fracture 178:245–258 16. Réthoré J, de Borst R, Abellan MA (2008) A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks. Comput Mech 42:227–238 Secchi and Schrefler Asia Pacific Journal on Computational Engineering 2014, 1:8 Page 21 of 21 http://www.apjcen.com/content/1/1/8 17. Mohammadnejad T, Khoei AR (2013) Hydromechanical modelling of cohesive crack propagation in multiphase porous media using extended finite element method. Int J Numer Anal Meth Geomech 37:1247–1279 18. Kraaijeveldt F, Huyghe JM, Remmers JJC, de Borst R (2013) 2-D mode one crack propagation in saturated ionized porous media using partition of unity finite elements. J Appl Mech 80:020907-1-12 19. Secchi S, Simoni L, Schrefler BA (2004) Cohesive fracture growth in a thermoelastic bi-material medium. Comput Struct 82:1875–1887 20. Kraaijeveld F (2009) Propagating discontinuities in ionized porous media, Dissertation. Eindhoven University of Technology 21. Kraaijeveldt F, Huyghe JM, Remmers JJC, de Borst R, Baaijens FPT (2014) Shearing in osmoelastic fully saturated media: a mesh-independent model. Engineering Fracture Mechanics. in press 22. Barenblatt GI (1959) The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses: axially-symmetric cracks. J Appl Math Mech 23:622–636 23. Remij EW, Pizzoccolo F, Remmers JJ, Smeulders D, Huyghe JM (2013) Nucleation and mixed-mode crack propagation in porous material. ASCE, Poromechanics V, pp 2260–2269. doi:10.1061/9780784412992.247 24. Pizzocolo F, Hyughe JM, Ito K (2013) Mode I crack propagation in hydrogels is stepwise. Eng Fract Mech 97:72–79 25. Dugdale DS (1960) Yielding of steel sheets containing slits. J Mech Phys Solids 8:100–104 26. Hilleborg A, Modeer M, Petersson PE (1976) Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 6:773–782 27. Camacho GT, Ortiz M (1996) Computational modelling of impact damage in brittle materials. Int J Solids Struct 33:2899–2938 28. Lewis RW, Schrefler BA (1998) The finite element method in the static and dynamic deformation and consolidation of porous media. Wiley, Chichester 29. Witherspoon PA, Wang JSY, Iwai KJE, Gale JE (1980) Validity of cubic law for fluid flow in a deformable rock fracture. Water Resour Res 16:1016–1024 30. Li XD, Wiberg NE (1998) Implementation and adaptivity of a space-time finite element method for structural dynamics. Comp Methods Appl Mech Engrg 156:211–229 31. Secchi S, Simoni L, Schrefler BA (2008) Numerical difficulties and computational procedures for thermo-hydro- mechanical coupled problems of saturated porous media. Comput Mech 43:179–189 32. Secchi S, Simoni L (2003) An improved procedure for 2-D unstructured Delaunay mesh generation. Adv Eng Softw 34:217–234 33. Secchi S, Simoni L, Schrefler BA (2007) Numerical procedure for discrete fracture propagation in porous materials. Int J Num Anal Methods Geomech 31:331–345 34. Zhu JZ, Zienkiewic OC (1988) Adaptive techniques in the finite element method. Com Appl Num Methods 4:197–204 35. Rank E, Katz C, Werner H (1983) On the importance of the discrete maximum principle in transient analysis using finite element methods. Int J Num Methods Engng 19:1771–1782 36. Bae JS, Krishnaswamy S (2001) Subinterfacial cracks in bimaterial systems subjected to mechanical and thermal loading. Eng Fract Mech 68:1081–1094 37. ICOLD (1999) Fifth international benchmark workshop on numerical analysis of dams. Theme A2, Denver, Colorado 38. Doglioni C, Barba S, Carminati E, Riguzzi F (2013) Fault on-off fluids response. Geoscience Frontiers, pp 1–14. doi.org/10.1016/j.gsf.2013.08.004 39. Kelbert A, Schultz A, Egbert G (2009) Global electromagnetic induction constraints on transition-zone water content variations. Nature 460:1003–1006. doi:10.1038/nature08257 40. Mohammadnejad T, Khoei AR (2013) An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model. Finite Elements in Analysis and Design 73:77–95 doi:10.1186/2196-1166-1-8 Cite this article as: Secchi and Schrefler: Hydraulic fracturing and its peculiarities. Asia Pacific Journal on Computational Engineering 2014 1:8. Submit your manuscript to a journal and benefi t from: 7 Convenient online submission 7 Rigorous peer review 7 Immediate publication on acceptance 7 Open access: articles freely available online 7 High visibility within the fi eld 7 Retaining the copyright to your article Submit your next manuscript at 7 springeropen.com
Journal
Asia Pacific Journal on Computational Engineering – Springer Journals
Published: May 16, 2014