Coupled hydromechanical‐fracture simulations of nonplanar three‐dimensional hydraulic fracture propagation

This paper presents an algorithm and a fully coupled hydromechanical‐fracture formulation for the simulation of three‐dimensional nonplanar hydraulic fracture propagation. The propagation algorithm automatically estimates the magnitude of time steps such that a regularized form of Irwin's criterion is satisfied along the predicted 3‐D fracture front at every fracture propagation step. A generalized finite element method is used for the discretization of elasticity equations governing the deformation of the rock, and a finite element method is adopted for the solution of the fluid flow equation on the basis of Poiseuille's cubic law. Adaptive mesh refinement is used for discretization error control, leading to significantly fewer degrees of freedom than available nonadaptive methods. An efficient computational scheme to handle nonlinear time‐dependent problems with adaptive mesh refinement is presented. Explicit fracture surface representations are used to avoid mapping of 3‐D solutions between generalized finite element method meshes. Examples demonstrating the accuracy, robustness, and computational efficiency of the proposed formulation, regularized Irwin's criterion, and propagation algorithm are presented.     

Developing constitutive models from EPR‐based self‐learning finite element analysis

A constitutive model that captures the material behavior under a wide range of loading conditions is essential for simulating complex boundary value problems. In recent years, some attempts have been made to develop constitutive models for finite element analysis using self‐learning simulation (SelfSim). Self‐learning simulation is an inverse analysis technique that extracts material behavior from some boundary measurements (eg, load and displacement). In the heart of the self‐learning framework is a neural network which is used to train and develop a constitutive model that represents the material behavior. It is generally known that neural networks suffer from a number of drawbacks. This paper utilizes evolutionary polynomial regression (EPR) in the framework of SelfSim within an automation process which is coded in Matlab environment. EPR is a hybrid data mining technique that uses a combination of a genetic algorithm and the least square method to search for mathematical equations to represent the behavior of a system. Two strategies of material modeling have been considered in the SelfSim‐based finite element analysis. These include a total stress‐strain strategy applied to analysis of a truss structure using synthetic measurement data and an incremental stress‐strain strategy applied to simulation of triaxial tests using experimental data. The results show that effective and accurate constitutive models can be developed from the proposed EPR‐based self‐learning finite element method. The EPR‐based self‐learning FEM can provide accurate predictions to engineering problems. The main advantages of using EPR over neural network are highlighted.     

A fully coupled element‐free Galerkin model for hydro‐mechanical analysis of advancement of fluid‐driven fractures in porous media

Hydraulic fracturing (HF) of underground formations has widely been used in different fields of engineering. Despite the technological advances in techniques of in situ HF, the industry uses semi‐analytical tools to design HF treatment. This is due to the complex interaction among various mechanisms involved in this process, so that for thorough simulations of HF operations a fully coupled numerical model is required. In this study, using element‐free Galerkin (EFG) mesh‐less method, a new formulation for numerical modeling of hydraulic fracture propagation in porous media is developed. This numerical approach, which is based on the simultaneous solution of equilibrium and continuity equations, considers the hydro‐mechanical coupling between the crack and its surrounding porous medium. Therefore, the developed EFG model is capable of simulating fluid leak‐off and fluid lag phenomena. To create the discrete equation system, the Galerkin technique is applied, and the essential boundary conditions are imposed via penalty method. Then, the resultant constrained integral equations are discretized in space using EFG shape functions. For temporal discretization, a fully implicit scheme is employed. The final set of algebraic equations that forms a non‐linear equation system is solved using the direct iterative procedure. Modeling of cracks is performed on the basis of linear elastic fracture mechanics, and for this purpose, the so‐called diffraction method is employed. For verification of the model, a number of problems are solved. According to the obtained results, the developed EFG computer program can successfully be applied for simulating the complex process of hydraulic fracture propagation in porous media. Copyright © 2016 John Wiley & Sons, Ltd.     

A practical and efficient numerical scheme for the analysis of steady state unconfined seepage flows

The scaled boundary finite‐element method (SBFEM), a novel semi‐analytical technique, is applied to the analysis of the confined and unconfined seepage flow. This method combines the advantages of the finite‐element method and the boundary element method. In this method, only the boundary of the domain is discretized; no fundamental solution is required, and singularity problems can be modeled rigorously. Anisotropic and nonhomogeneous materials satisfying similarity are modeled without additional efforts. In this paper, SBFE equations and solution procedures for the analysis of seepage flow are outlined. The accuracy of the proposed method in modeling singularity problems is demonstrated by analyzing seepage flow under a concrete dam with a cutoff at heel. As only the boundary is discretized, the variable mesh technique is advisable for modeling unconfined seepage analyses. The accuracy, effectiveness, and efficiency of the method are demonstrated by modeling several unconfined seepage flow problems. Copyright © 2011 John Wiley & Sons, Ltd.     

Simulation of non‐planar three‐dimensional hydraulic fracture propagation

Hydraulic fracturing is the method of choice to enhance reservoir permeability and well efficiency for extraction of shale gas. Multi‐stranded non‐planar hydraulic fractures are often observed in stimulation sites. Non‐planar fractures propagating from wellbores inclined from the direction of maximum horizontal stress have also been reported. The pressure required to propagate non‐planar fractures is in general higher than in the case of planar fractures. Current computational methods for the simulation of hydraulic fractures generally assume single, symmetric, and planar crack geometries. In order to better understand hydraulic fracturing in complex‐layered naturally fractured reservoirs, fully 3D models need to be developed. In this paper, we present simulations of 3D non‐planar fracture propagation using an adaptive generalized FEM. This method greatly facilitates the discretization of complex 3D fractures, as finite element faces are not required to fit the crack surfaces. A solution strategy for fully automatic propagation of arbitrary 3D cracks is presented. The fracture surface on which pressure is applied is also automatically updated at each step. An efficient technique to numerically integrate boundary conditions on crack surfaces is also proposed and implemented. Strongly graded localized refinement and analytical asymptotic expansions are used as enrichment functions in the neighborhood of fracture fronts to increase the computational accuracy and efficiency of the method. Stress intensity factors with pressure on crack faces are extracted using the contour integral method. Various non‐planar crack geometries are investigated to demonstrate the robustness and flexibility of the proposed simulation methodology. Copyright © 2014 John Wiley & Sons, Ltd.     

Thermal recovery from a fractured medium in local thermal non‐equilibrium

Thermal recovery from a hot dry rock (HDR) reservoir viewed as a deformable fractured medium is investigated with a focus on the assumption of local thermal non‐equilibrium (LTNE). Hydraulic diffusion, thermal diffusion, forced convection and deformation are considered in a two‐phase framework, the solid phase being made by impermeable solid blocks separated by saturated fractures. The finite element approximation of the constitutive and field equations is formulated and applied to obtain the response of a generic HDR reservoir to circulation tests. A change of time profile of the outlet fluid temperature is observed as the fracture spacing increases, switching from a single‐step pattern to a double‐step pattern, a feature which is viewed as characteristic of established LTNE. A dimensionless number is proposed to delineate between local thermal equilibrium (LTE) and non‐equilibrium. This number embodies local physical properties of the mixture, elements of the geometry of the reservoir and the production flow rate. All the above properties being fixed, the resulting fracture spacing threshold between LTNE and LTE is found to decrease with increasing porosity or fluid velocity. The thermally induced effective stress is tensile near the injection well, illustrating the thermal contraction of the rock, while the pressure contribution of the fracture fluid is negligible during the late period. Copyright © 2012 John Wiley & Sons, Ltd.     

Boundary Element based Discontinuous Deformation Analysis

A Boundary Element based Discontinuous Deformation Analysis (BE‐DDA) method is developed by implementing the improved dual reciprocity boundary element method into the open close iterations based DDA. This newly developed BE‐DDA is capable of simulating both the deformation and movement of blocks in a blocky system. Based on geometry updating, it adopts an incremental dynamic formulation taking into consideration initial stresses and dealing with external concentrated and contact forces conveniently. The boundaries of each block in the discrete blocky system are discretized with boundary elements while the domain of each block is divided into internal cells only for the integration of the domain integral of the initial stress term. The contact forces among blocks are treated as concentrated forces and the open–close iterations are applied to ensure the computational accuracy of block interactions. In the current method, an implicit time integration scheme is adopted for numerical stability. Three examples are used to show the effectiveness of the algorithm in simulating block movement, sliding, deformation and interaction of blocks. At last, block toppling and tunnel stability examples are conducted to demonstrate that the BE‐DDA is applicable for simulation of blocky systems. Copyright © 2016 John Wiley & Sons, Ltd.     

Large‐scale poroelastic fractured reservoirs modeling using the fast multipole displacement discontinuity method

An effective approach to modeling the geomechanical behavior of the network and its permeability variation is to use a poroelastic displacement discontinuity method (DDM). However, the approach becomes rather computationally intensive for an extensive system of cracks, particularly when considering coupled diffusion/deformation processes. This is because of additional unknowns and the need for time‐marching schemes for the numerical integration. The Fast Multipole Method (FMM) is a technique that can accelerate the solution of large fracture problems with linear complexity with the number of unknowns both in memory and CPU time. Previous works combining DDM and FMM for large‐scale problems have accounted only for elastic rocks, neglecting the fluid leak‐off from the fractures into the matrix and its influence on pore pressure and stress field. In this work we develop an efficient geomechanical model for large‐scale natural fracture networks in poroelastic reservoirs with fracture flow in response to injection and production operations. Accuracy and computational performance of the proposed method with those of conventional poroelastic DDM are compared through several case studies involving up to several tens of thousands of boundary elements. The results show the effectiveness of the FMM approach to successfully evaluate field‐scale problems for the design of exploitation strategies in unconventional geothermal and petroleum reservoirs. An example considering faults reveals the impact of reservoir compartmentalization because of sealing faults for both geomechanical and flow variables under elastic and poroelastic rocks. Copyright © 2015 John Wiley & Sons, Ltd.     

Accelerating strategies to the numerical simulation of large‐scale models for sequential excavation

In this paper, a novel combination of well‐established numerical procedures is explored in order to accelerate the simulation of sequential excavation. Usually, large‐scale models are used to represent these problems. Due to the high number of equations involved, the solver algorithm represents the critical aspect which makes the simulation very time consuming. The mutable nature of the excavation models makes this problem even more pronounced. To accomplish the representation of geometrical and mechanical aspects in an efficient and simple manner, the proposed solution employs the boundary element method with a multiple‐region strategy. Together with this representational system, a segmented storage scheme and a time‐ordered tracking of the changes form an adequate basis for the usage of fast updating methods instead of frontal solvers. The present development employs the Sherman–Morrison–Woodbury method to speed up the calculation due to sequential changes. The efficiency of the proposed framework is illustrated through the simulation of test examples of 2D and 3D models. Copyright © 2006 John Wiley & Sons, Ltd.     

Block fracturing analysis using nodal‐based discontinuous deformation analysis with the double minimization procedure

As a hybrid method, the nodal‐based discontinuous deformation analysis (NDDA) greatly improves the stress accuracy within each DDA block by coupling a well‐defined finite element mesh inside the DDA block; at the same time, the NDDA inherits the unique block kinematics of the standard DDA method. Each finite element mesh line inside the DDA block is treated as a potential crack, which enables the transformation of the block material from continuum to discontinuum through the tensile and shear fracturing mechanism. This paper introduces a double minimization procedure into the NDDA method to further improve the accuracy of the stresses evaluated at the finite element mesh lines and thus to obtain a more realistic fracture model. Three numerical examples are employed to demonstrate the improved stress accuracy by the implemented double minimization procedure and the accuracy and capability of the enhanced NDDA method in capturing brittle fracturing process. Copyright © 2013 John Wiley & Sons, Ltd.     

Semi‐analytical far field model for three‐dimensional finite‐element analysis

A challenging computational problem arises when a discrete structure (e.g. foundation) interacts with an unbounded medium (e.g. deep soil deposit), particularly if general loading conditions and non‐linear material behaviour is assumed. In this paper, a novel method for dealing with such a problem is formulated by combining conventional three‐dimensional finite‐elements with the recently developed scaled boundary finite‐element method. The scaled boundary finite‐element method is a semi‐analytical technique based on finite‐elements that obtains a symmetric stiffness matrix with respect to degrees of freedom on a discretized boundary. The method is particularly well suited to modelling unbounded domains as analytical solutions are found in a radial co‐ordinate direction, but, unlike the boundary‐element method, no complex fundamental solution is required. A technique for coupling the stiffness matrix of bounded three‐dimensional finite‐element domain with the stiffness matrix of the unbounded scaled boundary finite‐element domain, which uses a Fourier series to model the variation of displacement in the circumferential direction of the cylindrical co‐ordinate system, is described. The accuracy and computational efficiency of the new formulation is demonstrated through the linear elastic analysis of rigid circular and square footings. Copyright © 2004 John Wiley & Sons, Ltd.     

Hybrid FE‐DQ for dynamic analysis of multilayered half‐space subjected to concentrated point impulsive loading

A hybrid finite element method and differential quadrature method (DQM) is developed to estimate the dynamic response of two‐dimensional multilayered half‐spaces subjected to impulsive point loading. Nonreflecting absorbing boundary conditions consist of appropriate springs, and dampers are considered. The capabilities of the finite element method for solving boundary value problems with general domain, loading and systematic boundary treatment are combined with accurate and stable time marching capabilities of the DQM to develop an accurate and efficient numerical technique. The capability, efficiency, robustness and convergence of the DQM for solving the dynamic problem are demonstrated through numerical simulations of various half‐spaces with different time increments and layer arrangement. Also, comparison study when using Newmark's time integration scheme for the same problem is done. It can be concluded that the DQM as an unconditionally stable method is suitable for solving such a problem. Also, parametric study is performed to show the effect of the absorbing boundary conditions on the dynamic response. Copyright © 2012 John Wiley & Sons, Ltd.     

Modeling three‐dimensional hydraulic fracture propagation using virtual multidimensional internal bonds

Propagation of fractures, especially those emanating from wellbores and closed natural fractures, often involves Mode I and Mode II, and at times Mode III, posing significant challenges to its numerical simulation. When an embedded inclined fracture is subjected to compression, the fracture edge is constrained by the surrounding materials so that its true propagation pattern cannot be simulated by 2D models. In this article, a virtual multidimensional internal bond (VMIB) model is presented to simulate three‐dimensional (3D) fracture propagation. The VMIB model bridges the processes of macro fracture and micro bond rupture. The macro 3D constitutive relation in VMIB is derived from the 1D bond in the micro scale and is implemented in a 3D finite element method. To represent the contact and friction between fracture surfaces, a 3D element partition method is employed. The model is applied to simulate fracture propagation and coalescence in typical laboratory experiments and is used to analyze the propagation of an embedded fracture. Simulation results for single and multiple fractures illustrate 3D features of the tensile and compressive fracture propagation, especially the propagation of a Mode III fracture. The results match well with the experimental observation, suggesting that the presented method can capture the main features of 3D fracture propagation and coalescence. Moreover, by developing an algorithm for applying pressure on the fracture surfaces, propagation of a natural fracture is also simulated. The result illustrates an interesting and important phenomenon of Mode III fracture propagation, namely the fracture front segmentation. Copyright © 2012 John Wiley & Sons, Ltd.     

Forward Calculation of 2-D and 3-D Structures with a Cover by the Boundary Element Method while Using Electrical Methods

This paper develops the boundary element method,the authors employtwo-layered earth Green’s functions as the weighting functions of residual andderive boundary integral equations.The forward problems of point sources on2-D and 3-D structures with an influencing cover are solved by this method.The results show that this method markedly improves the original boundaryelement method.The features of the improved method are greater numericalaccuracy and much smaller systems of equations and thus considerable savingsfor the storage capacity of computers,allowing us to solve the above problemswith only ordinary microcomputers.The results in this paper extend the scopeof applying the boundary element method while using electrical methods forgeophysical prospecting.     

3‐D mesoscale simulation of crack‐permeability coupling in the Brazilian splitting test

In this paper, mesoscale hydromechanical simulations are performed to study (1) fracture features and (2) crack‐gas permeability coupling evolution in the context of the tensile splitting test. The mesostructure is based on a 2‐phase 3‐D representation of heterogeneous materials, such as concrete, where stiff aggregates are embedded into a mortar matrix. To take into account these heterogeneities without any mesh adaptation, a weak discontinuity is introduced into the strain field. In addition, a strong discontinuity is also added to take into account microcracking. This mechanical model is cast into the framework of the enhanced finite element method. Concerning the coupling with gas permeability, a double‐porosity method is used to simulate the flow through the cracks and the porosity. The apparent gas permeability is afterwards evaluated by a homogenization method. On the basis of finite element simulations, influence of aggregate size on ultimate crack opening, macroscopic ultimate tensile stress, total dissipated energy, and gas permeability evolution is numerically investigated. Furthermore, gas permeability evolution is also compared with experimental results from the literature. In addition, in the spirit of a sequential multiscale approach, macroscale gas permeability equations are identified from the hydromechanical results coming from the mesoscale computations. These equations lead to a relation between macroscale gas permeability evolution and crack opening. Besides, we show how the aggregate size influences the percolation threshold and that after this threshold, a cubic relation between macroscale gas permeability and crack opening is obtained.     

A high‐frequency open boundary for transient seepage analyses of semi‐infinite layers by extending the scaled boundary finite element method

A high‐frequency open boundary has been developed for the transient seepage analyses of semi‐infinite layers with a constant depth. The scaled boundary finite element equation of pore water pressure is formulated first in the frequency domain. With the eigenvalue problem, the equation can be decoupled into modal equations whose modal dynamic permeability equation can be determined. The continued fraction technique is adopted to formulate the continued fraction solution in the frequency domain. All constants in the solution are determined recursively at the high‐frequency limit. By introducing auxiliary variables and the continued fraction solution to the relationship between the prescribed seepage flow and the pore water pressure in the frequency domain, the open boundary condition is obtained. After transformed to the time domain, the open boundary condition is expressed as a system of fractional differential equations. No convolution integral is required. The accuracy of the analysis results increases with the increasing orders of continued fraction. Copyright © 2015 John Wiley & Sons, Ltd.     

A numerical method for the study of shear band propagation in soft rocks

This paper investigates the possibility of interpreting progressive shear failure in hard soils and soft rocks as the result of shear propagation of a pre‐existing natural defect. This is done through the application of the principles of fracture mechanics, a slip‐weakening model (SWM) being used to simulate the non‐linear zone at the tips of the discontinuity. A numerical implementation of the SWM in a computation method based on the boundary element technique of the displacement discontinuity method (DDM) is presented. The crack and the non‐linear zone at the advancing tip are represented through a set of elements, where the displacement discontinuity (DD) in the tangential direction is determined on the basis of a friction law. A residual friction angle is assumed on the crack elements. Shear resistance decreases on elements in the non‐linear zone from a peak value at the tip, which is characteristic of intact material, to the residual value. The simulation of a uniaxial compressive test in plane strain conditions is carried out to exemplify the numerical methodology. The results emphasize the role played by the critical DD on the mechanical behaviour of the specimen. A validation of the model is shown through the back analysis of some experimental observations. The results of this back analysis show that a non‐linear fracture mechanics approach seems very promising to simulate experimental results, in particular with regards to the shear band evolution pattern. Copyright © 2009 John Wiley & Sons, Ltd.     

Coupling of finite and hierarchical infinite elements: application to a non‐homogeneous cross‐anisotropic half‐space subjected to a non‐uniform circular loading

A method is presented for coupling cubic‐order quadrilateral finite elements with the finite side of a new coordinate ascent hierarchical infinite element. At a common side shared by a hierarchical infinite element and an arbitrary number of finite elements, the displacements are minimized in the least square sense with respect to the degrees‐of‐freedom of the finite elements. This leads to a set of equations that relate the degrees‐of‐freedom of the finite and hierarchical infinite elements on the shared side. The method is applied to a non‐homogeneous cross‐anisotropic half‐space subjected to a non‐uniform circular loading with Young's and shear moduli varying with depth according to the power law. A constant mesh constructed from coupled finite and hierarchical infinite elements is used and convergence is sought simply by increasing the degree of the interpolating polynomial. The displacements and stresses produced by conical and parabolic circular loads applied on the surface are obtained. The efficiency of the proposed method is demonstrated through convergence and comparison studies. New results produced by a frusto‐conical circular load applied on the surface of a half‐space made up of heavily consolidated London clay are provided. The non‐homogeneity parameter and degree of anisotropy are shown to influence the soil response. Copyright © 2012 John Wiley & Sons, Ltd.     

Transient seepage analysis in zoned anisotropic soils based on the scaled boundary finite‐element method

The scaled boundary finite‐element method, a semi‐analytical computational scheme primarily developed for dynamic stiffness of unbounded domains, is applied to the analysis of unsteady seepage flow problems. This method is based on the finite‐element technology and gains the advantages of the boundary element method as well. Only boundary of the domain is discretized, no fundamental solution is required and singularity problems can be modeled rigorously. Anisotropic and non‐homogeneous materials satisfying similarity are modeled with no additional efforts. In this study, firstly, formulation of the method for the transient seepage flow problems is derived followed by its solution procedures. The accuracy, simplicity and applicability of the method are demonstrated via four numerical examples of transient seepage flow – three of them are available in the literature. Homogenous, non‐homogenous, isotropic and anisotropic material properties are considered to show the versatility of the technique. Excellent agreement with the finite‐element method is observed. The method out‐performs the finite‐element method in modeling singularity points. Copyright © 2014 John Wiley & Sons, Ltd.     

Effects of non‐uniform traction and specimen height in the direct shear test on stress and deformation in a rock fracture

In the direct shear test (DST), an internal moment is distributed within the rock specimen by non‐coaxial shear loads applied to the specimen, which cause non‐uniform distributions of both the traction on the loading planes and the stress and deformation in the specimen. To examine the validity of the DST for a rock fracture and to clarify the effect of specimen height, both the stress and deformation in a fracture in the DST were analyzed for specimens with three different heights using a three‐dimensional finite element method with quadratic joint elements for a fracture model. The constitutive law of the fracture considers the dependence of the non‐linear behavior of closure on shear displacement and that of shear stiffness on normal stress and was implemented in simulation code to give a conceptional fracture with uniform mechanical properties to extract only the effect of non‐uniform traction on the stress and deformation in the fracture. The results showed that both normal and shear stresses are concentrated near the end edges of the fracture, and these stress concentrations decrease with a decrease in the specimen height according to the magnitude of the moment produced by the non‐coaxial shear loads. Furthermore, although closure is greater near the end edges of the fracture, where normal stress is concentrated, this concentration of closure is not so significant within the range of this study because of the non‐linear behavior of closure, that is, closure does not significantly increase with an increase in normal stress at large normal stresses. Copyright © 2012 John Wiley & Sons, Ltd.