A coupled hydro-mechanical analysis for prediction of hydraulic fracture propagation in saturated porous media using EFG mesh-less method

Computers and Geotechnics 55 (2014) 254–266

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

A coupled hydro-mechanical analysis for prediction of hydraulic fracture propagation in saturated porous media using EFG mesh-less method Mohammad Norouz Oliaei a,⇑, Ali Pak b, Kenichi Soga c a

Department of Civil Engineering, Tarbiat Modares University, Tehran 14115-143, Iran Department of Civil Engineering, Sharif University of Technology, Tehran 11365-9313, Iran c Engineering Department, University of Cambridge, Cambridge CB2 1PZ, UK b

a r t i c l e

i n f o

Article history: Received 30 October 2012 Received in revised form 15 July 2013 Accepted 3 September 2013

Keywords: Hydraulic fracture Propagation Initiation Mesh-less EFG Coupled hydro-mechanical analysis

a b s t r a c t The details of the Element Free Galerkin (EFG) method are presented with the method being applied to a study on hydraulic fracturing initiation and propagation process in a saturated porous medium using coupled hydro-mechanical numerical modelling. In this EFG method, interpolation (approximation) is based on nodes without using elements and hence an arbitrary discrete fracture path can be modelled. The numerical approach is based upon solving two governing partial differential equations of equilibrium and continuity of pore water simultaneously. Displacement increment and pore water pressure increment are discretized using the same EFG shape functions. An incremental constrained Galerkin weak form is used to create the discrete system of equations and a fully implicit scheme is used for discretization in the time domain. Implementation of essential boundary conditions is based on the penalty method. In order to model discrete fractures, the so-called diffraction method is used. Examples are presented and the results are compared to some closed-form solutions and FEM approximations in order to demonstrate the validity of the developed model and its capabilities. The model is able to take the anisotropy and inhomogeneity of the material into account. The applicability of the model is examined by simulating hydraulic fracture initiation and propagation process from a borehole by injection of fluid. The maximum tensile strength criterion and Mohr–Coulomb shear criterion are used for modelling tensile and shear fracture, respectively. The model successfully simulates the leak-off of fluid from the fracture into the surrounding material. The results indicate the importance of pore fluid pressure in the initiation and propagation pattern of fracture in saturated soils. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Hydraulic fracturing (HF) refers to the parting of geomaterials due to excessively high fluid pressure. The most common application of HF is enhance oil recovery, in which increase in functional permeability of the ground is achieved by the formation of permeable fractures. In geotechnical engineering, HF is one of the causes of cracks in earth and rockfill dams. It is also used for grouting to create impermeable barriers or to control ground settlements. In spite of the significant advances in the technique of in situ hydraulic fracturing, its design still requires a good deal of engineering judgment and practical experience. That is, the ability to determine the fracture path, shape, dimensions and fracture conductivities is still not fully developed. This is due to the complex interaction among the different mechanisms that are involved in ⇑ Corresponding author. Tel.: +98 21 82884395; fax: +98 21 82884914. E-mail addresses: (M.N. Oliaei).

[email protected],

[email protected]

0266-352X/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2013.09.001

hydraulic fracturing, namely ground deformation, fluid flow and fracturing process. As a result, a numerical model capable of analysing different aspects of hydraulic fracturing can be a valuable tool to overcome some uncertainties in the design of hydraulic fracturing in soils. Since capturing the pattern of arbitrary fracture growth is an essential ingredient in HF modelling, a numerical scheme for modelling the hydraulic fracturing path is developed in this research using a fully coupled hydro-mechanical element free Galerkin (EFG) mesh-less method. The FEM for modelling complex problems in engineering science is well established. However, it is not without shortcomings. One of them is crack propagation problems. The use of a mesh in modelling these problems creates difficulties in the treatment of discontinuities and large distortions that do not coincide with the original mesh lines. This is an inherent weakness of FEM, which is related to its element-based shape functions. One solution for such a problem is remeshing the domain using adaptive computational algorithm. However, projection of field variables between the meshes used in successive stages of an

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

analysis leads to degradation of the accuracy. For large threedimensional problems, the computational cost of remeshing at each step of an analysis becomes prohibitively expensive. Furthermore, it is almost impossible to automatically remesh finite elements for an arbitrarily growing crack. Therefore most of the FEM methods developed to date involve interactive analysis where the user guides the remeshing after each step of analysis [1]. An effective solution to avoid mesh dependency is to use a mesh-less method. This method does not require any element for shape function construction and hence provides no connectivity between elements and nodes. The principle attractions of meshless methods are their simplified adaptability and ability to cope with moving boundaries and discontinuities. There are a number of mesh-less (or mesh-free) methods that have been proposed in the past (e.g. [2]). Among these methods some of them are well known and achieved remarkable success in recent years. These methods can be classified into two major categories based on the shape function construction approach: i Finite integral representation methods, i.e. Smooth Particle Hydrodynamics or SPH [3,4]. ii Finite series representation methods, i.e. the Element Free Galerkin (EFG) method [5] based on MLS (Moving Least Squares) and the Point Interpolation Method (PIM) [6]. The MLS shape functions adopted in EFG are both consistent and compatible; the PIM shape functions, however, are consistent but not compatible. The SPH shape functions are compatible but not consistent near the boundary of the problem domain. Consistency is the capability of the field function approximation method to reproduce the fields of lowest orders of complete polynomials at any point in the problem domain. Compatibility refers to the continuity of the approximation on the boundaries between subdomains, based on which shape the functions are constructed. Both consistency and compatibility affect the accuracy and convergence of the numerical results [2]. The first attempt to apply mesh-less strategies to a soil–water coupled problem was made by Modaressi and Aubert [7] using the coupled EFG–FEM technique. In their work, displacements of solid skeleton are modelled by the standard FEM, whereas fluid pressures are determined by element–free nodes. This attempt was followed by Murakami et al. [8] with an elasto-plastic constitutive relation within the infinitesimal strain theory and by Murakami and Arimoto [9] and Arimoto et al. [10] in finite strain with the aim to capture the localization phenomenon in saturated soil. Murakami et al. [11] also used EFG in the analysis of a punch problem for a soft soil foundation which has stress singularity problems. They used the forward difference approximation (fully explicit scheme) for discretization in time domain. Wang et al. [12] proposed to use the backward difference approximation (fully implicit scheme) in time domain to avoid spurious ripple effect. Oliaei et al. [13] further examined the stability criteria for coupled soil deformation and pore fluid flow problems. Generally, two major approaches are possible for modelling of fracture: one is considering discrete fractures and the other is the use of smeared cracking. Discrete fracture modelling is the preferred technique for simulating the propagation of a small number of individual fractures over the use of smeared cracking since it does not suffer from stress locking, spurious modes and directional bias [14]. Among available mesh-free methods, the EFG method developed by Belytschko et al. [5] is adopted in this study for modelling of discrete hydraulic fracture propagation. The EFG method has been used for the simulation of growing crack problems, both static and dynamic [1,15]. Its key advantage for modelling this class of problems using the EFG method

255

is that it does not require any element connectivity data and does not suffer much degradation in accuracy when nodal arrangements are irregular [15]. Because of the gridless character of the method, a growing crack can be modelled simply by extending the free surfaces in the domain which correspond to the crack. This dramatically simplifies the modelling of moving cracks, since it eliminates the need for remeshing. Another advantage of the EFG method compared with the FEM is the dependent variable and its gradient are continuous in the entire domain, and hence post processing to obtain smooth gradient fields is unnecessary [15]. The shape functions obtained by the EFG approximation based on nodes (not elements) are both consistent and compatible. They are of a higher order compare to those used in the ordinary FEM, which effectively induces more accurate approximations compared to solutions obtained by a distorted mesh in FEM. The main objective of this paper is to investigate the validity and capability of the EFG formulation developed in this study in order to predict hydraulic fracture propagation. The coupled analysis of soil deformation and fluid flow in inhomogeneous/anisotropic saturated porous media is conducted. In Section 2, the EFG shape function construction and the MLS approximation are briefly described. In Section 3, the weak forms are developed through a global equilibrium at each time-step. Then spatial variables (i.e. displacement and pore water pressure increments) are discretized by the same EFG shape functions. A fully implicit scheme in time domain is used to avoid spurious ripple effect. An algorithm for numerical solution is proposed to solve hydraulic fracture initiation and propagation problems based on EFG. The (so-called) diffraction technique [16] is introduced for the construction of mesh-less approximation with discontinuities and nonconvex boundaries. Section 4 presents the numerical analysis of 2D coupled consolidation problem and crack propagation problem using fracture mechanics. Hydraulic fracture initiation pressure and the pattern of hydraulic fracture propagation are computed and results are compared to the closed-form and/or numerical (FEM) solutions to examine the accuracy of the new algorithm. Conclusions are given in the final section.

2. Element Free Galerkin (EFG) shape function and Moving Least Squares (MLS) approximation The Element Free Galerkin (EFG) method establishes system of equations for the whole problem domain without the use of a predefined mesh. To represent (not discretize) the problem domain and its boundaries, a set of nodes are scattered within the problem domain as well as on the boundaries. Therefore, the construction of shape functions is only based on the nodes as shown in Fig. 1. The EFG method employs the MLS approximation to construct the shape function and approximate the desired function as given by Eq. (7) in Fig. 1. As shown in Eqs. (3) and (8) through Eq. (15) of Fig. 1, the EFG shape function and its spatial derivatives within the problem domain are constructed from two_ components: (i) a weight function associated with each node-ðW ðX  XI ÞÞ; and (ii) a basis function usually consisting of a polynomial-(PT(X)). The weight function is nonzero over a small subdomain around a node, which is called its support. The support of the weight function defines a node’s influence domain, which is the area (or volume) over a particular node contributing to the approximation. The overlap of the nodal influence domains defines the nodal connectivity. P(X) – Eqs. (2) and (3) in Fig. 1 – is a vector of basis functions that contains most often of monomials of the lowest orders to ensure minimum completeness. Most mesh-less weight functions are bell-shaped and the following equations have been proposed in the past [2].

256

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

Fig. 1. EFG shape function and its derivatives construction flow chart in the presented EFG code. based on the work of Belytschko et al. [5].

257

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

(a) Cubic spline weight function (W3):

8 2 2 3 > < 3  4r þ 4r _ _ W ðX  XI Þ ¼ W ðrÞ ¼ 43  4r þ 4r 2  43 r 3 > : 0 _

r:ðqv Þ  G ¼ 

for r 6 12 for

1 2


ð16Þ

for r > 1

(b) Quadratic spline(weight function (W4): _

W ðX  XI Þ ¼ W ðrÞ ¼

_

(2) Continuity equation of fluid flow

1  6r 2 þ 8r 3  3r 4

for r 6 1

0

for r > 1

ð17Þ

e (c) Exponential weight ( function (W ):

_

W ðX  XI Þ ¼ W ðrÞ ¼

eðr=aÞ 0

2

for r 6 1 for r > 1

ð18Þ

where a is constant (a = 0.3 has been recommended by Liu [2]). In Eqs. (16)–(18), r = kX  XIk/(Rid)I where kX  XIk is the distance between the evaluation point and the node I and (Rid)I is the radius of the influence domain of the Ith node. These weight functions have been tested and have proved to work very well for many applications and have basically been adopted for various kinds of problems [2]. MLS, which is originated by mathematicians for data fitting and surface construction, is often termed local regression and loss [17]. Nayroles et al. [18] were the first to use MLS procedure to construct shape functions for the Diffuse Element Method (DEM). DEM was modified by Belytschko et al. [5], as the EFG method, where the MLS approximation is also employed. The MLS approximation – with a bell-shaped weight function – has two major features that make it popular: (1) the approximated field function is continuous and smooth in the entire problem domain; and (2) it is capable of producing an approximation with the desired order of consistency. One attractive property of the MLS approximants is that their continuity is related to the continuity of the weight function; therefore, a low order polynomial basis, e.g., a linear basis – which has used in this study – can be used to generate highly continuous approximations by choosing an appropriate weight function [5]. The shape functions of this method are global and can be used all over the domain. Although the MLS approximation is both consistent and compatible, the use of the MLS approximation produces shape function that do not possess the Kronecker delta function property, which implies that one cannot impose essential boundary conditions in the same way as in conventional FEM. In order to solve this problem, the penalty method proposed by Liu and Yang [19] is used in this study to create constrained Galerkin weak form for imposing essential boundary conditions. The use of the penalty method produces equation systems of the same dimensions such that the conventional FEM produces for the same number of nodes, and the modified stiffness matrix is still positive definite, banded, and symmetric. The treatment of boundary conditions is as simple as it is in conventional FEM. However, the disadvantage of penalty method is the way of selection a suitable large positive number for each set of the equations. 3. Formulation 3.1. Strong and weak forms For soil deformation–water flow coupled problems, there are two sets of governing equations [20,21]: (1) Equilibrium equation in incremental form

Drij;j þ Dbi ¼ 0

ð19Þ

where rij is the total stress tensor, and bi is the body force vector.

@ðnqÞ @t

ð20Þ

where q is the density of fluid, v is the Darcy’s velocity vector of fluid flow, G is the fluid mass flux from sink or source, n is the porosity of soil mass and t is time. In this derivation, the pore fluid is considered to be water with constant density. This assumption is the same as the incompressibility equation of solid–water mixture in Biot’s consolidation theory [22]. Considering that sink or source term can be considered later as a boundary condition, Eq. (20) is written in the following form.

r:ðv Þ ¼ 

@n @t

ð21Þ

The soil deformation is governed by Terzaghi’s effective stress principle.

Dr0ij ¼ Drij  aDpdij

ð22Þ

where r is the effective stress tensor, a is Biot’s coefficient (1 for soils), Dp is the pore water pressure increment, and dij is kronecker delta. The constitutive law for soil skeleton in incremental form is written as 0 ij

1 Dr0ij ¼ Dijkl ðDekl  cs dkl DpÞ 3

ð23Þ

where Dijkl is the soil skeleton stiffness matrix, Dekl is the total strain increment tensor, and cs is the compressibility of solid particles of soil. The constitutive relationship for water flow in porous media is Darcy’s law.

v i ¼ K ij

  p zþ

c

ð24Þ

;j

where Kij is the permeability tensor of soil skeleton, z is the elevation head, p is the pore water pressure, and c is the unit weight of water. Boundary conditions for soil skeleton boundary and fluid boundary can be written as follows



on Cu  ½0; 1Þ

i ui ¼ u

rij nj ¼ ti on Ct  ½0; 1Þ 

 p¼p v i ¼ v i

on  ½0; 1Þ on  ½0; 1Þ

for soil skeleton boundary

for fluid boundary

ð25Þ

ð26Þ

 i is the boundary value of diswhere ui is the displacement vector, u placement, Cu is the displacement boundary, nj is the unit normal vector at the boundary, t i is the boundary value of traction, Ct is  is the boundary vathe traction boundary, p is the pore pressure, p lue of pore pressure, Cp is the pore pressure boundary, vi is Darcy’s velocity vector, v i is the boundary value of Darcy’s velocity, and Cv is the velocity (flux) boundary. The initial conditions can be written as



ui ¼ ui jt¼0

on X 0

p ¼ pjt¼0

on X 0

ð27Þ

where X is the domain. By applying the weighted residual method on Eq. (19) and inserting Eqs. (22), (23), and (25) into Eq. (19), the following constrained Galerkin weak form of the equilibrium equation is obtained.

258

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

R

R

R

R dðLDuÞT Dijkl Dekl dX  X dðDuÞT Dbi dX  Ct dðDuÞT Dt i dC þ X dðLDuÞT adij DpdX R R  ÞT apu ðDu  Du  ÞdC  X dðLDuÞT ð1=3Þcs Dijkl dkl DpdX  d Cu ð1=2ÞðDu  Du R T  þ  Du   Þ bpu ðDu  þ  Du   ÞdC ¼ 0 d Cs ð1=2ÞðDu X

ð28Þ

where d(Du) is the test function, L is the differential operator, Du is  is the prescribed incremental the incremental displacements, Du displacements on the essential boundary, apu is the penalty factor for the weak form of equilibrium equation to penalize the difference between the incremental displacements of the MLS approximation and the prescribed incremental displacements on the essential boundary, Cs is the interface of the different materials in an inho þ and Du   are the incremental displacements mogeneous body, Du in the different materials but on the two sides of interfaces and bpu is the penalty factor for the weak form of equilibrium equation  þ and Du  . to penalize the difference between Du According to the time domain discretization methods, the following relation can be used for a field function f in the time interval of [t, t + Dt]:

f ¼ ð1  hÞft þ hftþDt ¼ ft þ hDf

ð29Þ

where h can vary from zero (fully explicit scheme) to 1.0 (fully implicit scheme). The approximation is unconditionally stable when h P 0:5, but for any value of h – 1 the numerical solution can exhibit a spurious rippling effect [12]. This time integration is applied to Eq. (24). By using the weighted residual method and inserting Eqs. (24) and (26) into Eq. (21), the weak form for the continuity equation of pore water is expressed as R

dðDpÞT ðv i ni ÞdC þ

R

T

R

T

ðK ij =cÞpt;j dX R R R T T ÞT app ðDp  Dp ÞdC þ X dðLp DpÞ hðK ij =cÞDp;j dX þ X dðDpÞ ð@n=@tÞdX  d Cp ð1=2ÞðDp  Dp R  þ  Dp  ÞT bpp ðDp  þ  Dp  ÞdC ¼ 0 d Cs ð1=2ÞðDp Cv

X dðLp DpÞ

K i3 dX þ

X dðLp DpÞ

ð30Þ

where d(Dp) is the test function, Lp is the differential operator, Dp is  is the prescribed pore water the pore water pressure increment, Dp pressure increment on the essential boundary, app is the penalty factor for the weak form of continuity equation to penalize the difference between the pore water pressure increment of MLS approximation and the prescribed pore water pressure increment on the essential boundary, Cs is the interface of the different materials in þ and Dp  are the pore water pressure an inhomogeneous body, Dp increments in the different materials but on the two sides of interfaces and bpp is the penalty factor for the weak form of continuity þ and Dp  . equation to penalize the difference between Dp

the Ith node. The procedure of constructing uI and its spatial derivatives is given in Fig. 1. Differential operator matrices L and Lp are given by

2

0

0

@=@y

@=@y

0

@=@x

0

@=@z

0

@=@x

6 LT ¼ 4 0 0

Displacement increments Du and pore water pressure increment Dp at any time and at any point are approximated using Eq. (7) in Fig. 1.

8 9h 2 uI > n < Dux > = X 6 h Du ¼ Duy ¼ 4 0 > > : ; I 0 Duz

Dph ¼

n X

/I DpI

0

uI 0

9 38 > n < DuxI > = X 7 0 5 DuyI ¼ UI DuI > > I uI : DuzI ; 0

ð31Þ

ð32Þ

I

where Duh is the approximated displacement increments, Dph is the approximated pore water pressure increment, DuI is the displacement increments of the Ith node, DpI is the pore water pressure increment of the Ith node, n is the number of nodes in the neighbourhood of the point of interest (Gauss point) which the point is in the influence domain of node, uI is the MLS shape function of

@=@z

3

7 0 5 @=@y @=@x

ð33Þ

@=@z

LTp ¼ ½ @=@x @=@y @=@z 

ð34Þ

By using Eqs. (31)–(34), the products of LDuh and LoDph become

2

uI;x

6 0 6 n 6 X 6 0 h 6 LDu ¼ 6u I 6 I;y 6 4 0

uI;z

0

3

0

9 0 7 78 DuxI > 7> n < = X 7 uI;z 7 DuyI ¼ BI DuI 7 > 0 7> : ; I 7 DuzI uI;y 5

uI;y 0

uI;x uI;z

ð35Þ

uI;x

0

2 3 uI;x n n X X 6u 7 Lp Dph ¼ BpI DpI 4 I;y 5DpI ¼

uI;z

I

ð36Þ

I

where uI,x, uI,y, uI,z are the MLS shape function spatial derivatives of the Ith node. As Eq. (30) is discretized, app is a scalar, where as apu in Eq. (28) is a diagonal matrix of penalty factors.

2

apux

apu ¼ 6 4 0

0

0 0

3

0

7 0 5

apuy

ð37Þ

apuz

The penalty factors can be functions of coordinates and they can be different from each other, but in practice an identical constant of a large positive number is assigned for each set of the equations [19]. bpp and bpu have the same form of app and apu, but their values may be different. The appropriate selection of the penalty factors is given in Oliaei et al. [13] for coupled soil deformation and pore fluid flow problems. Substituting Eqs. (31), (32), (35), and (37) into Eq. (28), the first system equation can be obtained for the equilibrium equation. Similarly, by substituting Eqs. (32) and (36) into Eq. (30), and with consideration of

Dev ¼

n X 

uI;x uI;y

I

3.2. Numerical discretization

0

8 9 > DuxI > n = X < uI;z DuyI ¼ C I DuI > > : ; I DuzI

ð38Þ

the second system equation can be obtained for the continuity equation. These two equations constitute the final system of discrete equations for the entire problem domain in the EFG method. These equations should be solved simultaneously for a fully coupled analysis. The matrix equation in coupled form is therefore written as

"

K 11 þ K au þ K bu

K 12

K 21

K 22 þ K ap þ K bp

#

DU DP

(

 ¼

DF u þ DF au DF p þ DF ap

) ð39Þ

where DU and DP are the global displacement increments parameter vector and the global pore pressure increment parameter vector, respectively. K11, K12, K21 and K22 are the parts of global coefficient matrix. The non-diagonal terms in the matrix [K] of Eq. (39) represent the coupling terms in the analysis; K12 represents force induced by pore pressure and K21 represents fluid flow caused by the ground deformation. DFu and DFp are the global force increments vector and the global fluid flux increment vector,

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

respectively. They are the standard terms in the conventional soil– water coupled problem derivations. In EFG method the nodal matrices and vectors can, in fact, be viewed as the basic components for assembling the global matrices and vectors. It has been clearly shown that EFG method operates on nodes, in contrast to the operation on elements in the FEM, where the element matrices and vectors are the basic components for assembling. The nodal matrices and vectors for the above mentioned terms can be written as

K 11ij ¼

K 12ij ¼

K 21ij ¼

Z X

Z X

Z

BTi DBj dX Z

BTi ð1=3Þcs Dm/j dX 

X

BTi am/j dX

ð41Þ

/i C j dX

ð42Þ

X

K 22ij ¼ hDt

DF u i ¼

ð40Þ

Z X

Z X

BTpi ðK w =cÞBpj dX

UTi DbdX þ

DF pi ¼ Dt

Z X

Z Ct

ð43Þ

UTi DtdC

BTpi K w3 dX  Dt

Z X

BTpi ðK w =cÞBpi pi dX  Dt

ð44Þ Z Cv

/i v Ti ni dC ð45Þ

where ‘m’ represents dij in a vector form i.e. m ¼ h 1 1 1 0 0 0 iT , Kw represents the permeability tensor and K Tw3 ¼ ½ 0 0 K z  introduces the datum effect. In DF pi the first and second terms are the flows due to changes in velocity, while the third term indicates the effect of specified flux on the boundaries. In Eq. (39), K au and K ap are the global penalty matrices assembled using the nodal penalty matrices, whereas DF au and DF ap are the global penalty increments vectors assembled using the nodal penalty increments vectors. The nodal matrices and vectors for these terms are defined as follows.

K auij ¼ K apij ¼

Z Cu

Z

UTi apu Uj dC

ð46Þ

/i app /j dC

ð47Þ

Cp

DF aui ¼ DF api ¼

Z Cu

Z

UTi apu Du dC

ð48Þ

 dC /i app Dp

ð49Þ

Cp

Furthermore, in Eq. (39), K bu and K bp are the global penalty matrices assembled using the nodal penalty matrices in other to connect the different materials in an inhomogeneous body. The nodal matrices for these terms are defined as follows. b

Z

b

Z

K uij ¼

T

Cs

K pij ¼

Cs

½Uþi  Ui  bpu ½Uþj  Uj dC

ð50Þ

½/þi  /i bpp ½/þj  /j dC

ð51Þ

 where /þ i and /i are the shape functions created using the nodes in different materials that have influence on the quadrature points along the interfaces by a nonpenetration rule [23].

259

The procedure sequence for the numerical algorithm of the developed EFG code is illustrated in Fig. 2. 3.3. Discontinuities in approximations With the use of node-based interpolation techniques, mesh-less methods offer great opportunities to handle problems of complex geometries. However, since there is no connectivity between nodes, difficulties arise in determining the domain of influence of a node. The influence domain usually takes a simple shape, for example circular, for simple geometries of problems. The influence of a node on any other point is directly computed by the distance between the node and the point via the use of the weight functions such as Eq. (16) through Eq. (18). When a domain involves discontinuities such as cracks, an influence domain may contain irregular boundary fragments. Hence the computation of nodal weights simply based on physical distance can be erroneous as the discontinuity of the field variables caused by the complex boundary is not accounted for (Fig. 3). A number of techniques have been reported on the construction of mesh-less approximation with discontinuities and nonconvex boundaries. There are three typical methods: the visibility, diffraction and transparency methods [16]. These methods have been used for defining the profile of the influence domain and computing the weight of influence in mesh-less approximations with irregular boundaries. The visibility method is simple, but the weight function by visibility criterion is discontinuous within the influence domain around the tip of a discontinuity line. As a consequence the shape functions constructed are also discontinuous, which is undesirable in a mesh-less approximation. Moreover, the visibility criterion may not be very easy to adapt to concave boundaries such as holes. The transparency method was developed to smooth a mesh-less approximation around the tip of a discontinuity line with a varying degree of transparency, from completely transparent at the tip of discontinuity to completely opaque at a distance away from the tip. As a consequence, the shape functions constructed are depended to degree of transparency, which is undesirable in a mesh-less approximation. The diffraction method was used in the EFG code developed for this study. The development of this method has been motivated by the way light diffracts around a sharp corner, such as crack tip, but the equations used in constructing the domain of influence and the weight function bear almost no relationship to the equations of diffraction. This technique applies only to polar-type weight functions, such as Eq. (16) through Eq. (18), where the weights are defined as a function of a single distance-related weight parameter. The mechanism of the diffraction method is illustrated in Fig. 4, where the straight line represents a crack. By the diffraction method, the crack line is treated as opaque; therefore, the length of the ray is evaluated by a path that passes around the crack tip. According to Belytschko et al. [16], the weight parameter s associated with point P is computed by

sðxÞ ¼

 k s1 þ s2 ðxÞ s0 ðxÞ s0 ðxÞ

ð52Þ

where s0, s1, and s2 are shown in Fig. 4, and k is a user-defined parameter. The domain of influence is determined by the condition that the weight function vanishes. As a result, the domain of influence, as given by the criterion w(s) > 0, contracts around the crack tip (see Fig. 4). The diffraction method has been applied to nonconvex boundaries such as hole. Fig. 5 illustrates such an example, where the parameter s is constructed from the lengths of two line segments that just graze the boundary.

260

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

Fig. 2. The procedure sequence for the numerical algorithm of the developed EFG code.

The major benefit of the diffraction method is that the weight function and shape function are continuous within the influence domain but are discontinuous across the crack line. The derivatives of the shape function are multiple valued at the crack tip but this should cause no difficulties. However, it results in the cost of lower efficiency in computing the derivatives of shape functions. Care should also be exercised because when node and point are very close, i.e., s0(x) is close to zero, the weight parameter s becomes infinite, causing difficulty in computation. Therefore, it is required that all nodes have a minimum distance to the crack surface.

4. Numerical results Two examples were used to verify the developed EFG code. The first example considered the numerical validity of the EFG meshless method by solving fully coupled hydro-mechanical problem [24] and the second example evaluated the numerical capability of the method to model discrete fracture growth in rock [25]. The last one was considered as an example of interaction between ground deformation and fracturing process. To investigate the numerical stability and the desired level of accuracy of the EFG solution for coupled hydro-mechanical problems, a sensitivity analysis was conducted to examine the numerical parameters that influence the solutions [13]. Based on the findings, the parameters are chosen to guarantee the accuracy of

the EFG solution for soil deformation–fluid flow coupled problems. A summary of recommended parameters is listed in Table 1. Then, the application of the developed EFG code to determine hydraulic fracture initiation pressure and to predict hydraulic fracture propagation path in homogeneous and inhomogeneous saturated soils will be conducted and complex interactions among ground deformation, fluid flow and fracturing process will be examined. 4.1. Prediction of hydraulic fracture propagation in saturated soils Hydraulic fracture propagation from a hole inside a saturated cemented soil sample for low permeability (undrained) and high permeability (drained) conditions using tensile and shear criteria is simulated. The maximum tensile strength criterion and Mohr– Coulomb shear criterion are used for modelling tensile and shear fracture, respectively. This example is designed to compare the EFG results with analytical, experimental and field results and to investigate the capability of the developed EFG code. However, the shear criterion in the developed EFG code can predict very good the hydraulic fracture initiation pressure (similar to tension criterion), using shear fracture criterion for modelling of hydraulic fracture propagation in saturated soils by the developed EFG code has been shown that the rate of extension in plastic shear failure zone in the sample with respect to the rate of crack propagation is extremely high; which means the shear band cannot

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

261

Fig. 2 (continued)

Fig. 3. Highly irregular influence domain [2].

created properly. It coincides to the field and laboratory test results for boreholes in cemented soils which the propagation of tensile fracture is almost dominant with respect to shear fracture. Therefore, tensile fracture criterion has been presented in this study. (i) Homogeneous soil Fig. 6 shows the geometry of the problem simulating laboratory testing of hydraulic fracturing. A plane strain analysis is

considered. The characteristic parameters for the problem are as follow: The initial total outside and inside confining pressures of the sample r = 200 kPa, The initial pore pressure in the sample P = 100 kPa, Young’s modulus E = 10 MPa and Poisson’s ratio v = 0.45. The permeability is K = 1010 m/s for the low permeability case and K = 105 m/s for the high permeability case. The borehole radius is a = 1 cm and the dimension of the sample is W = 10 cm. The tensile strength of cemented soil ðr0t  rt Þ is 100 kPa. It is assumed that the shear strength of the cemented soil is very large compared to the tensile strength and hence the soil always fracture in tension. When the initial stresses on the homogeneous sample are isotropic, increasing the pressure of the injected fluid in the borehole up to hydraulic fracture initiation pressure, causes tensile crack to propagate from the borehole in any radial directions. However, in reality, deformation around the borehole is not uniform and the fracture initiates at a localized point or zone. Therefore, in this example, after fracturing, a microcrack (1 mm length) is considered at the borehole wall in positive x direction to initiate hydraulic fracture propagation from this point due to hydraulic pressure and leak-off from crack and also borehole. The nodal arrangement of the EFG model is shown in Fig. 7(a). It consists of (i) a regular arrangement in the domain, (ii) a radial– tangential arrangement on the borehole wall, (iii) a radial arrangement at the crack tip which is progressive with the crack tip [1] and (iv) an arrangement on crack surfaces after propagation. The arrangement of the background integration mesh is shown in

262

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266 Table 1 Guidelines to guarantee the accuracy of the EFG solution for coupled hydromechanical problems [13]. 2

2

h =ðcv DtÞ 6 2 for 1D, h =ðcv DtÞ 6 3 for 2D 612

app = 10

(K/c)

apu = 1048E Dimension (radius) of influence domain = (1.25–1.75) average local diagonal nodal spacing Number of Gauss points = at least about 5 times the number of nodes (It reflects the number of integration cells and the order of integration. These two have to be balanced) Type of weight function Cubic spline weight function Quadratic spline weight function Exponential weight function (with constant parameter P 0:3Þ Type of weight function has not much effect on results

Fig. 4. The diffraction method to determine the influence domain [2] (a) Definition of parameters (b) Domain of influence.

Fig. 6. Geometry of the model for the problem of hydraulic fracture initiation pressure determination and propagation path prediction.

Fig. 5. The diffraction method applied to concave boundary [2].

Fig. 7(b). The Gauss points arrangement is 3  3 for each cell, except in three rows and three columns in the middle of the model, which is 5  5 for each cell. After applying initial stresses and pore pressures on the sample, a constant fluid flux is applied on the borehole wall such that the flow rate to the borehole is Q = 1  109 m3/s. Nodal displacements and nodal pore pressures are calculated in each time step; then the pore pressure calculated on the borehole wall is applied as fluid traction on the borehole wall in the next time step together with the constant flux on the borehole wall. For the outer boundary, the flow boundary conditions are zero flux for the low permeability (undrained) case and constant pressure (P = 100 kPa) for the high permeability (drained) case. After fracturing, constant flux on the borehole wall and crack surfaces is applied in each time step. Then,

nodal displacements and nodal pore pressures are calculated at the end of each time step and the procedure will be repeated as mentioned above. Although a constant flux is assumed around the crack for simplicity, more investigation in this area is required to study the effect of bulk fluid flow inside the crack. The applied flux on the borehole wall acts as traction on the wall until the domain around the borehole reaches to tensile failure as shown in Fig. 8(a) for the low permeability soil case. The excess pore pressure does not generate inside the domain and hence the pore pressure in the domain is approximately constant over time as shown in Fig. 8(b). Therefore, it is essentially an undrained cavity expansion in elastic medium. Using cylindrical cavity expansion theory, the hydraulic fracture pressure for tensile fracture in undrained condition can be obtained as follows [26]:

Pf ¼ 2r3i  u0 þ r0t

ð53Þ

where Pf is hydraulic fracture pressure, r3i is initial applied minimum principle stress, u0 is initial pore pressure and r0t is tensile strength of soil. The fracture initiation pressure of the high permeability soil case was smaller than that of the low permeability soil case as shown in Fig. 9(a). This is due to leak-off of fluid to the soil sample when the soil is permeable. The pore fluid pressure around the borehole is increasing due to fluid injection in the borehole as

263

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

Undrained Condition - Tension Criterion

Pore Pressure (Kpa)

200 175 150 125 100 75 50 25 0 0

10

20

30

40

50

60

time (s)

Hydraulic Fracture Pressure (Kpa)

Fig. 8(b). The variation of pore pressure around the cavity due to fluid injection for undrained condition under tensile failure criterion.

Fig. 7(a). Nodal arrangement of the EFG model (half of the model).

Drained Condition - Tension Criterion

H.F.

400

300

200

100

0

0

50

100

150

200

250

300

time (s) Fig. 9(a). The variation of cavity surface traction due to fluid injection for drained condition under tensile failure criterion.

Drained Condition - Tension Criterion

Pore Pressure (Kpa)

400

Fig. 7(b). The arrangement of the background integration mesh (a quarter of the model).

300

200

100

0

0

50

100

150

200

250

300

time (s) Fig. 9(b). The variation of pore pressure around the cavity due to fluid injection for drained condition under tensile failure criterion.

Hydraulic Fracture Pressure (Kpa)

Undrained Condition - Tension Criterion

H.F.

500 400 300 200 100 0

0

10

20

30

40

50

60

time (s) Fig. 8(a). The variation of cavity surface traction due to fluid injection for undrained condition under tensile failure criterion.

shown in Fig. 9(b), which reduces the effective stress of the soil. Increasing of the traction on the borehole wall together with the

pore pressure in the soil will cause tensile failure on the borehole wall. Bjerrum et al. [27] have evaluated hydraulic fracture pressure at a cylindrical cavity in permeable soil to interpret the fracturing conditions around a driven piezometer. Vertical cracks in the radial direction from the piezometer develop when the circumferential effective stress reduces to the negative sign of the tensile strength of the soil. Hence, the hydraulic fracture pressure for tensile fracture in drained condition can be obtained as follows:

 Pinj ¼ u0 þ

1

v

  1 ½rt þ ð1  aÞr0h0 

ð54Þ

where Pinj is injection pressure, u0 is the initial pore pressure, v is Poisson’s ratio, rt is the tensile strength, r0h0 is the initial horizontal effective stresses and a is a disturbance factor for the changes in

264

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

Table 2 The comparison between hydraulic fracture initiation pressures for analytical solutions and EFG results, using tensile failure (facture) criterion and under different permeability (drainage) conditions (U: undrained, D: drained, T: tensile). Analytical

425 323

400 344

Drained Condition Undrained Condition

400 300 200 100

Undrained Condition - Tension Criterion

0 0.01

450

0.02

0.03

0.04

0.05

X (m) 440

Fig. 11. The fluid pressure variation along the crack axis and in the time step before the last time step for tensile failure criterion.

430 420 410

Drained Condition - Tension Criterion 350 51

52

53

54

55

Fig. 10(a). Fluid pressure at crack entrance during hydraulic fracture propagation for undrained condition under tensile failure criterion.

Drained Condition - Tension Criterion

Pore Pressure (Kpa)

400

time (s)

Hydraulic Fracture Pressure (Kpa)

Crack Pressure (Kpa)

EFG

U–T D–T

Hydraulic Fracture Pressure (Kpa)

Pf (kPa)

Tension Criterion 500

300 250

t=264

200

t=272

150

t=280

100 50

350

0

0.01

340

0.02

0.03

0.04

0.05

X (m) 330

Fig. 12(a). The pore pressure in domain during hydraulic fracture propagation (I axis) for drained condition under tensile failure criterion.

320 310 300 264

Drained Condition - Tension Criterion 266

268

270

272

274

276

278

350

280

Fig. 10(b). Fluid pressure at crack entrance during hydraulic fracture propagation for drained condition under tensile failure criterion.

circumferential effective stress due to piezometer installation. a varies between 0.4 and 1.1 based on soil compressibility. The comparison of hydraulic fracture initiation pressure for analytical/semi-analytical solutions and EFG results, using tensile failure (fracture) criterion and under different permeability (drainage) conditions, are presented in Table 2. The computed fracturing pressures are obtained from Fig. 8(a) and 9(a), whereas the theoretical values are calculated from Eqs. (53) and (54). The EFG results are very good agreement with closed form analytical solutions. Figs. 10(a) and 10(b) shows the fluid pressure at node A – Fig. 7 – during HF propagation for low and high permeability soil cases. Node A is separated into nodes A1 and A2 at the entrance of crack after fracturing. The figure shows variation of fluid pressure during hydraulic fracture propagation until reaching to the outer boundary of the sample. The propagation of crack is faster and takes place under higher stresses in low permeability soil than in high permeability soil. Fig. 11 shows the fluid pressure variation inside the crack, along the crack propagation axis at one step before the crack reaches to the outer boundary. For the low permeability soil case, the fluid pressure is approximately constant along the crack axis, whereas

Pore Pressure (Kpa)

time (s)

300 250

t=264

200

t=272 150

t=280

100 50 0 0.01

0.02

0.03

0.04

0.05

X (m) Fig. 12(b). The pore pressure in domain during hydraulic fracture propagation (J axis) for drained condition under tensile failure criterion.

for the high permeability soil case, there is pore pressure gradient due to infiltration of fluid from crack to the surrounding material. The progress of infiltration (leak-off) is illustrated in Figs. 12(a)–(c) in which the pore pressures along the three lines (I, J and K) in Fig. 7(a) are plotted at three different times (at the beginning of crack propagation or crack initiation, at the middle of crack propagation and at the end of crack propagation) for darined case. As expected, the pore pressure is rising with time and the magnitude of the increase is greater at locations close to the crack axis and the borehole. (ii) Heterogeneous soil

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

Drained Condition - Tension Criterion

Pore Pressure (Kpa)

350 300 250

t=264

200

t=272 150

t=280

100 50 0

0.01

0.02

0.03

0.04

0.05

X (m) Fig. 12(c). The pore pressure in domain during hydraulic fracture propagation (K axis) for drained condition under tensile failure criterion.

265

selected to be 0.001 times of its Young’s modulus. The confining pressures for this problem are as the same as the homogeneous soil in (i). The initial nodal arrangement of the EFG model is shown in Fig. 14 by black circles. A radial nodal arrangement at the crack tip (progressive with the crack position) and also an arrangement on crack surfaces after crack propagation are added after fracturing takes place. The arrangement of the background integration mesh is 10  10, and the Gauss points arrangement is 5  5 for each cell. The hydraulic fracture propagation path for the heterogeneous soil problem is presented in Fig. 14 by white triangular markers. After hydraulic fracture initiation, the crack propagates from the borehole towards the weaker layers. During crack propagation, when the crack reaches to a new layer, crack deviates from the previous path. The amount of deviation depends on the difference between characteristic parameters (strength and stiffness) of two different layers. When the crack reaches to a very stiff layer with high strength, the crack cannot propagate through this layer and it grows parallel to the layer toward the outer boundary.

5

E = 50 MPa

4

E = 30MPa

5. Conclusions

3 2 1

E = 20MPa

0 -5

-4

-3

-2

-1

0

1

2

3

4

5

-1 -2 -3

E = 10MPa

E = 5MPa

-4 -5

E = 100MPa Fig. 13. Geometry of the model for the problem of prediction of hydraulic fracture propagation path in inhomogeneous domain.

In this paper, the details of the Element Free Galerkin (EFG) mesh-less method and its numerical implementation have been presented to study the hydraulic fracture initiation and propagation using coupled hydro-mechanical numerical modelling in geotechnical engineering. EFG is an effective method to approximate and numerically determine the field variables (displacement and pore water pressure). Unlike other mesh-less methods, the EFG has a simple approach for construction of shape functions and its spatial derivatives. This is due to nature of its polynomial basis and weight function. The weak forms are developed through a global equilibrium at each time-step using the coupled expression. Essential boundary conditions can be easily implemented using penalty method. A fully implicit scheme in time domain is used to achieve stability and accuracy and avoid spurious ripple effect. The validity, feasibility and capability of the EFG method are illustrated by examples of hydraulic fracture initiation and propagation in saturated soils. While the EFG method is attractive in many problems because of the simplicity of discretization (only nodes in and on boundaries are needed) and avoidance of post-processing (results are already smooth), its greatest appeal is in problems with large change in geometry, such as progressive fracture, where finite element methods require considerable remeshing. Hence, by using the EFG method, the arbitrary discrete fracture path can be simulated elegantly. Finally, the involved fluid response in crack is hard to describe, i.e. around the crack tip. Although a constant flux is assumed in this paper around the crack for simplicity, more investigation in this area is required to study the effect of bulk fluid flow inside the crack. Research is currently under way by the authors. References

Fig. 14. The prediction of hydraulic fracture propagation path in inhomogeneous domain.

Fig. 13 shows the geometry of the heterogeneous soil. Tensile fracture criterion was adopted in high permeability soil. The domain consists of six types of materials having a range of stiffnesses from 5 MPa to 100 MPa. The tensile strength of each layer is

[1] Belytschko T, Lu YY, Gu L. Crack propagation by element free Galerkin methods. Eng Fract Mech 1995;51(2):295–315. [2] Liu GR. Mesh Free Methods-Moving beyond the finite element method 2003. CRC Press. [3] Lucy L. A numerical approach to testing the fission hypothesis. Astron J 1977;82:1013–24. [4] Gingold RA, Monaghan JJ. Smooth particle hydrodynamics: theory and applications to non-spherical stars. Mon Not R Astron Soc 1977;181:375–89. [5] Belytschko T, Lu YY, Gu L. Element free Galerkin methods. Int J Numer Meth Eng 1994;37:229–56. [6] Liu GR, Gu YT. A point interpolation method. In: Proc. 4th Asia-Pacific conference on computational mechanics, Singapore; December, 1999. p. 1009– 14. [7] Modaressi H, Aubert P. Element-free Galerkin method for deforming multiphase porous media. Int J Numer Meth Eng 1998;42:313–40.

266

M.N. Oliaei et al. / Computers and Geotechnics 55 (2014) 254–266

[8] Murakami A, Kawabata H, Aoyama S. EFGM analysis for saturated soil. In: Desai CS, et al., editors. Computer methods and advances in geomechanics, vol. 1; 2001. p. 153–6. [9] Murakami A, Arimoto S. Localized behaviour of saturated soil via element-free strategy. In: Proc. 12th Asian regional conference on soil mechanics and geotechnical engineering, 1, Balkema; 2003. p. 925–8. [10] Arimoto S, Murakami A, Setsuyasu T, Nishiyama T. Numerical analysis of saturated soil by EFG within finite strains. Trans JSME 2004;70(691):442–8. in Japanese. [11] Murakami A, Setsuyasu T, Arimoto S. Mesh-free method for soil-water coupled problem within finite strain and its numerical validity. Soils Found 2005;45(2):145–54. [12] Wang JG, Liu GR, Lin P. Numerical analysis of Biot’s consolidation process by radial point interpolation method. Int J Solids Struct 2002;39:1557–73. [13] Oliaei M, Soga K, Pak A. Some numerical issues using Element Free Galerkin mesh-less method for coupled hydro-mechanical problems. Int J Numer Anal Meth Geomech 2009;33(7):915–38. [14] Lim IL, Johnston IW, Choi SK. A finite element code for fracture propagation analysis within elasto-plastic continuum. Eng Fract Mech 1996;53(2):193–211. [15] Belytschko T, Lu YY, Gu L, Tabbara M. Element free Galerkin methods for static and dynamic fracture. Int J Solids Struct 1995;32(17–18):2547–70. [16] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. Comput Meth Appl Mech Eng 1996;139:3–47.

[17] Lancaster P, Salkauskas K. Surfaces generated by moving least squares methods. Math Comput 1981;37:141–58. [18] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992;10:307–18. [19] Liu GR, Yang KY. A penalty method for enforce essential boundary conditions in element free Galerkin method. In: Proc. 3rd HPC Asia’98, Singapore; 1998. p. 715–21. [20] Bathe KJ. Finite element procedures in engineering analysis. Englewood Cliffs, N.J.: Prentice-Hall, Inc.; 1982. p. 114–94 [chapter 4]. [21] Thomas GW. Principles of Hydrocarbon Reservoir simulation. Tabir Publishers; 1977. [22] Biot MA. General theory of three-dimensional consolidation. J Appl Phys 1941;12:155. [23] Liu GR, Yang KY. A new meshless method for stress analysis for solids and structure. In: Proc. 4th Asia-Pacific conference on computational mechanics, Singapore; 1999. [24] Gibson RE, Schiffman RL, Pu SL. Plane strain and axially symmetric consolidation of a clay layer on a smooth impervious base. Q J Mech Appl Mech 1970;23(4):505–20. [25] Ingraffea AR. Theory of crack initiation and propagation in rock. In: Atkinson BK, editor. Fracture mechanics of rock. Academic Press Geology Series; 1987. [26] Mitchell JK, Soga K. Fundamentals of soil behaviour. 3rd ed. John Wiley & Sons, Inc; 2005. [27] Bjerrum L, Nash JKTL, Kennard RM, Gibson RE. Hydraulic fracturing in field permeability testing. Geotechnique 1972;22(2):319–32.