Modelling fracture formation and propagation in geo-materials

The modelling of fracture formation and propagation in geo-materials is a subject of interest in several problems and applications. In this work, finite elements with high aspect ratio are used to study the hydraulic fracture process in the context of conventional continuum constitutive (stress-strain) relationships based on the damage theory. The deformable porous material is solved in a coupled manner by considering a fully hydro-mechanical approach. The adopted approach is explained and validated against analytical and numerical solutions.


Introduction
The formation and propagation of discontinuities in porous media (i.e., in the form of fractures, cracks or fissures) correspond to a relevant and complex phenomena of interest for different engineering applications, some of them discussed as follows. For example, the potential migration of gas and/or pollutants flow through discontinuities can seriously jeopardize the performance of repositories for high-level nuclear waste and other type of barrier systems. The mechanical integrity of the geomaterials (i.e. rocks and soils) comprising a geological system is a key component in several problems involving energy production. For instance, the injection of fluids (e.g. CO2, H2O) at very high pressure is generally necessary to enhance oil production from reservoirs. This technique may trigger the formation of fractures and can also reactivate preexisting geological faults [1].
Furthermore, in unconventional fields oil and/or gas are trapped in very low porosity/permeability rocks (e.g., gas/oil shales) and special recovery operations (outside the conventional techniques) are required to produce hydrocarbons for this type of system. Hydraulic fracturing is the most common enhanced recovery technique used in shale gas fields and consists in injecting engineered fluids at very-high pressures into a reservoir via injection wells [e.g.2,3]. The high pressurized fluid breaks the rocks, forming the fractures and allowing the release of the trapped oil and gas. Similar techniques are required to assist oil production from recently found very-deep oil reservoirs. The formation and propagation of drying cracks in soils is also a complex problem involving the presence of evolving discontinuities in geomaterials [4].
The numerical modelling of developing discontinuities, as the ones in the engineering problems discussed above, it is a quite challenging problem that involves the coupling of several physical processes. One of the first theoretical/analytical approaches dealing with the formation of fractures correspond to the so-called PKN [2] and KGD models [3]. Both of them consider simplifying assumptions, and therefore, they are applicable only when specific conditions and hypotheses are considered. As for the numerical modelling for dealing with this type of problem, several simulations techniques have been developed to investigate the formation of fractures. As for example, models based on the boundary element method [5,6]; extended finite elements (X-FEM) [7,8,9]; and also zero-thickness interface elements [10].
In this work, a finite-element (FE) based framework incorporating high aspect ratio (HAR) elements is used to simulate the formation and propagation of fractures in geo-materials. A key component of the proposed approach is a tensile damage constitutive model that incorporates the characteristic length of the problem, allowing the regularization of the numerical solution of geomaterials with softening behaviour. The proposed technique is general and can be potentially used to tackle any of the engineering problems discussed above. In this paper is verified for the near-K solution of the leak offtoughness dominated regime [3,4] and validated against the numerical results obtained by Carrier and Granet [10].
This paper is organized as follows, first the main aspects of the adopted approach is presented; then the main equations of the framework are introduced; afterwards the application case is studied; and finally the paper closes with main conclusions of this research.

Adopted approach
The proposed framework is based on the finite element program CODE_BRIGHT [11], which is a powerful platform capable of dealing with coupled Thermo-Hydro-Mechanical (THM) processes and interactions typically observed in problems in geological media involving simultaneous multiphysics actions. A multiphase/multispecies mathematical formulation consisting on three main set of equations is adopted, namely: i) balance equations (e.g., mass balance of species, balance of internal energy, and momentum balance); ii) constitutive equations (e.g., Darcy's and Fourier's laws, stress-strain relationship), and iii) equilibrium restrictions (e.g., Henry's law). More details can be found elsewhere [11]. In this work, the following main phenomena are considered: water flow (via liquid phase advection) and geomaterials deformation controlled by effective stresses (including the influence of load history). This formulation has been widely validated and applied to solve different coupled THM problems in geological media [1,11,12].
CODE_BRIGHT was developed to deal with continuous porous media and therefore was not able to tackle problems involving the formation and propagation of discontinues. In this work, CODE_BRIGHT was upgraded to simulate evolving fractures by implementing in it the approach suggested by Manzoli et al. [13,14,15] that is based on the inclusion of HAR finite elements with tailored mechanical and hydraulic constitutive models to represent the behaviour of discontinuities. The mechanical behaviour of the HAR elements is controlled by a tensile damage model, while the hydraulic response depends on the fracture aperture through the well know cubic law. As for the regular finite elements of the mesh (i.e. the bulk elements), the mechanical behaviour can be described by means of traditional constitutive models (e.g. a linear/non-linear elastic, or elasto-plastic models), while the hydraulic behaviour is modelled by the Darcy's generalized law.
This numerical technique allows to explicitly include the formation and propagation of fractures in the numerical analysis. CODE_BRIGHT uses GiD [16] or ParaView [17] to generate the mesh and to prepare all the information required for a finite element simulation. The following section presents the governing equations used to model the hydraulic fracturing problem.

Governing equations
The governing equations comprise the balance equations, constitutive relationship, and phase laws, together with boundary and initial conditions of the problem. The next sections show the main expressions used to model the fracturing phenomenon.

Balance equations
The equation of mass balance solid that controls porosity changes can be expressed as: where  denotes porosity, s  is the solid density and u is the velocity vector. The momentum balance equation for the porous material neglecting inertial terms is given by: where σ is the total stress vector, g is the gravity vector (i.e. g T ={0,0,-g}) and  is the density of the porous medium.
The water mass balance equation can be written as: where l  is the water density and l q is the liquid Darcy's velocity vector.

Constitutive equations
The mechanical constitutive models are introduced first and the hydraulic law afterwards.

Mechanical constitutive equation
The mechanical behaviour of the material is ruled by the effective stresses ( ' σ ) expressed in terms of total stresses and water pressure by: where   T 1,1,1, 0, 0, 0 m = is an auxiliary vector, l p is the liquid pressure and b is the Biot's coefficient, defined as with K and s K denoting the bulk moduli of the porous medium and the solid phase, respectively. A scalar tension damage model is adopted to tackle the dissipative process related to the formation of fractures. Tab. 1 contains the main equations associated with the damage model, which, as mentioned before, controls the behaviour of the HAR elements. In this work it is assumed that the mechanical behaviour of the bulk elements is controlled by a linear elastic model.

Hydraulic constitutive equations
It is assumed that Darcy's law governs the water flow in both, the solid materials and the discontinuities: where k is the intrinsic permeability tensor of the porous medium and l  is the water dynamic viscosity.
The permeability of the HAR-FE is computed in terms of the crack aperture (or width) by means of the well-known cubic law and the solid elements through the Kozeny-Carman equation in terms of the porosity.

Applications
The storage-toughness dominated regime of the hydraulic fracturing process is investigated in this section by using the proposed numerical technique.  Two cases were considered in order to study the effect of the isotropic (kx = ky = 5.0x10 -15 m 2 ) and orthotropic (kx =1.0x10 -16 m 2 , ky = 5.0x10 -15 m 2 ) rock permeability on fracture aperture, length and contour of pressure. The parameters of the near-K solutions (analytical results) were reported in Carrier and Granet [10]. Fig. 1 presents the results in terms of the aperture for the 1D and 2D diffusion cases associated with the orthotropic and isotropic rock permeability conditions, respectively, together with the analytic solution (available for the '1D-diffusion' case only). Fig. 2 presents the comparisons in terms of the predicted fracture length. In general terms, a very satisfactory performance of the proposed method incorporating HAR finite elements is observed when compared against the numerical and analytic solutions reported in Carrier and Granet [10].  The contour of liquid pressure at 100 t  s is shown in Figures 3 and 4. The unidimensional flow hypothesis is valid for the 1D-diffusion case because the fluid flow does not goes further than the fracture tip [10]. On the other hand, in the 2D-diffusion case, the flow is clearly bi-dimensional, as show in these figures.

Conclusions
In this work, a numerical approach capable of modelling the coupled displacement and pressure field problems in porous media with evolving discontinuities is suggested. The main equations and assumptions of the proposed framework are presented and discussed. The predictions obtained with the new approach are compared against already published analytic and numerical solutions involving the formation and propagation of fractures in geomaterials. The comparisons were very satisfactory, in general terms, showing the ability of the model to replicate the main trends and results involved in this type of problem. This is very promising and simple numerical technique that allows extending classical finite elements formulations and codes to deal with evolving discontinuities in geomaterials.
The financial support from NEUP-DOE, USA, through Award DE-NE0008762 is acknowledged. The first author would like to acknowledge CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico, Brazil) for funding his doctoral studies.