- Email: [email protected]

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol

A sequential implicit discrete fracture model for three-dimensional coupled ﬂow-geomechanics problems in naturally fractured porous media ⁎

Mostafa Moradia, Amir Shamlooa, , Alireza Daneh Dezfulib a b

Department of Mechanical Engineering, Sharif University of Technology, Azadi Ave., Tehran, Iran Department of Mechanical Engineering, Faculty of Engineering, Shahid Chamran University, Ahvaz, Iran

A R T I C L E I N F O

A BS T RAC T

Keywords: Fractured porous media Discrete-fracture model Coupled ﬂow-geomechanics Finite Element method

A sequential implicit numerical method based on discrete-fracture model and the Galerkin Finite Element method, for time-dependent coupled ﬂuid ﬂow and geomechanics problems in fractured subsurface formations is presented. Discrete-fracture model has been used to explicitly represent the fracture network inside porous media. The Galerkin Finite Element method with adaptive unstructured gridding is implemented to numerically solve the spatially three-dimensional and time-dependent problem. The presented method is validated with previously obtained solutions. Two problems are numerically solved by applying the presented methodology in a three-dimensional fractured petroleum reservoir under diﬀerent production and injection conditions. The modeling results demonstrate that changes in fracture aperture and permeability due to geomechanical stress, have a signiﬁcant impact on well performance and production rate in fractured reservoirs.

1. Introduction The stress dependent ﬂuid ﬂow is one of the challenging problems in many ﬁelds of geosciences. Interactions between solid media and the ﬂuid which is stored in it, lead to an inherently time-dependent coupled ﬂuid ﬂow and geomechanics problem. These interactions are due to changes in ﬂuid pore pressure (via injection or removal of the ﬂuid), which in turn cause the solid media to deform. As a result of this deformation, the characteristics of the solid media such as permeability (i.e. the ability of the media to conduct ﬂuid) change. The changes in solid media's characteristics alters the way that the ﬂuid ﬂows within it. This chain of events inside the porous media demonstrates the nature of the coupled ﬂuid ﬂow and geomechanics problem (Rutqvist and Stephansson, 2003). Disregarding the eﬀects of geomechanical parameters on ﬂuid ﬂow in petroleum reservoirs can lead to signiﬁcant economic and environmental damages. For instance, excessive water ﬂooding for oil recovery enhancement, without considering the geomechanical factors and how they activate or deactivate ﬂuid pathways has caused massive oil leaks in the past (Gu et al., 2014). Therefore it is important to pay careful attention to geomechanical eﬀects in fractured formations. Many previous studies exist in the literature, concerning the stresssensitivity eﬀects in conventional porous media (Pedrosa, 1986; Kilmer et al., 1987; Zhang and Ambastha, 1994; Chin et al., 2000; Wang et al., 2010; Qanbari and Clarkson, 2012, 2014; Shaoul et al., 2015).

⁎

However, coupled ﬂuid ﬂow and geomechanics study in fractured formations is a relatively new ﬁeld of study in reservoir simulation. Diﬀerent approaches have been taken by several researchers to present a method for solving the stress dependent ﬂuid ﬂow problem in fractured porous media. In one-way coupled method, the changes in ﬂuid properties aﬀect the geomechanics, but geomechanical information is not shared with the ﬂuid ﬂow. On the other hand, two-way coupled methods allow ﬂuid ﬂow and geomechanics to interact with each other. Numerical simulation of two-way coupled problems can be carried out using two distinct approaches: one is the fully implicit scheme, in which the ﬂow and geomechanics equations are solved simultaneously, and the second approach is called the sequential method. Fully implicit methods oﬀer a more accurate solution, but have large computational costs. However, sequential methods are more ﬂexible from a computational perspective, but might suﬀer from convergence and stability limitations (Castelletto et al., 2015; Gutierrez et al., 2001; Kim et al., 2012; Jalali and Dusseault, 2012; Settari and Mourits, 1998). Another approach which is called loose coupling, is a middle ground between one-way and two-way coupling method. In loose coupling, two sets of equations, one for ﬂuid ﬂow and the other for geomechanics are solved independently, but information is shared between the two at speciﬁed time intervals (Minkoﬀ et al., 2003). There are diﬀerent numerical methods to solve the ﬂuid ﬂow and the geomechanics problem in porous media. For the ﬂuid ﬂow problem,

Corresponding author. E-mail address: [email protected] (A. Shamloo).

http://dx.doi.org/10.1016/j.petrol.2016.12.027 Received 13 September 2016; Received in revised form 20 November 2016; Accepted 20 December 2016 0920-4105/ © 2016 Elsevier B.V. All rights reserved.

Please cite this article as: Moradi, M., Journal of Petroleum Science and Engineering (2016), http://dx.doi.org/10.1016/j.petrol.2016.12.027

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

ﬁnite element based methods have been widely used to simulate the coupled problem (Karimi-Fard and Firoozabadi, 2003; Monteagudo and Firoozabadi, 2004; Martin et al., 2005). Many other works based on the ﬁnite volume method have also been published (Karimi-Fard, 2004). For the solid mechanics part, ﬁnite element based methods have also been used extensively (Zienkiewicz and Taylor, 2005). Finite element based methods and hybrid ﬁnite element-ﬁnite volume methods have been employed in solving the coupled ﬂuid ﬂow and geomechanics problems (Gu et al., 2014; Garipov et al., 2016). Integral ﬁnite diﬀerence method is another approach that has been used in such problems (Hu et al., 2013). Another relatively new scheme, is the extended ﬁnite element method (XFEM) which is suitable for problems with local features that cannot be ﬁxed by mesh reﬁnement (Khoei et al., 2006; Lamb et al., 2013; Shamloo et al., 2005). Many researchers have investigated the coupled problem in fractured media. Some of these works have utilized the dual-porosity model (Chen and Tuefel, 2000; Segura et al., 2016). Discrete fracture model has been more widely used recently in coupled ﬂuid ﬂow and geomechanics problems (Garipov et al., 2016; McClure et al., 2016; Norbeck et al., 2014; Philip et al., 2005). Instead of solving the coupled system of equations, Moinfar et al. (2013) used in-situ stress conditions to estimate the stress ﬁeld within the reservoir, and utilized an experimental relation to compute the fracture aperture deformation due to the stress ﬁeld (Moinfar et al., 2013). Rutqvist et al. (2014) implemented a linked multicontinuum and crack tensor approach to simulate coupled ﬂuid ﬂow, geomechanics and solute transport in fractured rock (Rutqvist et al., 2014). Garipov et al. (2016) presented a fully implicit coupled ﬂuid ﬂow and geomechanics formulation for fractured formations. In their work, contact problem has been included to study the eﬀects of geomechanical stress on fracture aperture and hydrocarbon production rate (Garipov et al., 2016). In the present work, we aim to implement the discrete fracture model along with the Galerkin ﬁnite element model to present a sequential implicit time-dependent numerical method to solve the coupled ﬂuid ﬂow and geomechanics problem in subsurface naturally fractured formations. The presented method is sequential implicit, which means that ﬁrst, the ﬂuid ﬂow equation is solved implicitly. After obtaining the unknown ﬂuid ﬂow variable, i.e. pore pressure, geomechanical equations are explicitly solved based on the values acquired for pore pressure. We use the discrete fracture model (Karimi-Fard and Firoozabadi, 2003; Monteagudo and Firoozabadi, 2004) to represent the fracture system inside the porous media. In this work, an experimental relation is used to correlate the fracture aperture deformation with the calculated stress ﬁeld. Firstly, the theory and governing equations for ﬂuid ﬂow and geomechanical stress, fracture modeling and fracture aperture relation with stress ﬁeld are presented. In the third section, the numerical method and computational grid structure are explained. The method is then validated with the results obtained from similar research in section four. After model validation, some examples of the application of the presented method are oﬀered.

Fig. 1. Shear dilation.

• •

In-situ stress ﬁeld is constant at all times. Gravity eﬀects are neglected.

2.2. Governing equations The mechanical behavior of Porous media is described by the poroelasticity theory. According to this theory, the total stress is the sum of stress caused by ﬂuid pore pressure and the stress inside the rock skeleton (Biot, 1941; Chen et al., 1995).1

σij =σij′+αPδij

⎧ 1,& i = j δij=⎨ ⎩ 0, & i ≠ j

Where σ is the total stress tensor, σ′ is the biot eﬀective stress, P is ﬂuid pore pressure and α is the biot coeﬃcient. Fluid ﬂow equation is governed by the law of conservation of mass:

α

∂ϵ ∂P +S =−∇.v ∂t ∂t

(2)

Where ϵ is strain, v is ﬂuid velocity and S is storativity which is equal to:

S =ϕCf +(α −ϕ) Cr

(3)

Where, ϕ is the porosity of porous media, Cf is ﬂuid compressibility and Cr is rock skeleton compressibility. In Eq. (2), there are two unknown variables, P and v . However according to Darcy's law, ﬂuid velocity is related to pore pressure:

k v=− ∇P μ

(4)

Where k is the permeability of porous media and μ is the viscosity of ﬂuid. Combining, Eqs. (3) and (4), results in Eq. (5).

α

∂ϵ ∂P +S ∂t ∂t

k =∇.( ∇P ) μ

(5)

The stress which acts on a representative volume of the porous media is the total stress. The conservation of momentum in rock skeleton can be written as:

∇.σ −f =0 2. Mathematical formulation

(6)

Where f is the vector of external forces. Assuming that deformations are small and porous media is linearly elastic, eﬀective stress can be written in terms of strain:

In this section, the underlying theory and equations governing the coupled ﬂuid ﬂow and geomechanics problem are presented. The assumptions made in the method are also presented.

⎛ 2 ⎞ σij′=−⎜K − G⎟ ϵ v −2G ϵij i=jσij′=−2G ϵij i=j ⎝ 3 ⎠

2.1. Assumptions

• • •

(1)

(7)

In Eq. (7), K is the bulk modulus and G is the shear modulus of rock skeleton, ϵ v is the volumetric strain and is equal to:

The ﬂuid is single-phase and slightly compressible and has constant properties. The solid media is linearly elastic and its properties also remain constant. The fractures are two dimensional rectangular planes that have variable apertures that are functions of the applied stress ﬁeld.

ϵ v=ϵxx +ϵ yy +ϵzz

1

2

Tensile stress is positive.

(8)

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

In Eq. (13), Kn is the normal fracture stiﬀness. 2- Shear stress, τf , acting along the fracture makes the fracture surfaces to slide relative to each other and the normal component of this displacement is the change in fracture aperture. This mechanism is called “shear dilation” and it is shown in Fig. 1. Shear dilation does not always occur when a shear stress is applied to the fracture. Shear dilation starts when the shear stress is equal to the shear strength, σc . In other words, when τf < σc , shear stress causes tangential displacement to occur between the fracture surfaces, however normal displacement does not happen until τf > σc (Barton et al., 1983). The relation between shear stress and tangential displacement is:

Strain tensor is related to rock skeleton displacement vector by Eq. (9).

ϵ=1/2(∇u+∇T u)

(9)

Where u is the rock skeleton displacement vector. Rewriting Eq. (6) in terms of the displacement vector and then combining it with Eq. (1) lead to the equation shown below:

⎛ 1 ⎞ ⎜K + G⎟ ∇ϵkk +G∇2 u−α∇P=f ⎝ 3 ⎠

(10)

So far, there are four independent variables (i.e. three components of displacement vector and pore pressure) and four equations (Eqs. (5) and (10)). Note that ϵkk is a function of u, and therefore is not an independent unknown variable.

dt =

τf (14)

Ks

2.3. Discrete fracture model

Nevertheless, the tangential displacement that is associated with shear dilation, is related to shear stress as follows:

Flow modeling in fractured formations is substantially more complicated than conventional porous media. The main diﬃculty associated with modeling fractured porous media is that fractures have a strong eﬀect on the direction of ﬂow and the way ﬂuid is transported within the porous media. The diﬀerence between the size scale of fracture width and porous media renders simple numerical methods costly and diﬃcult (Jaﬀre et al., 2011). Therefore it is necessary to implement a special method to overcome this diﬃculty. There are diﬀerent approaches for modeling the system of fractures inside the porous media. Continuum models rely on averaging and homogenization methods, and therefore cannot model the geometry of the fracture system accurately. But on the hand, the discrete-fracture model represents each fracture individually. This allows us to include and model the system of fractures inside the porous media more accurately. In discrete-fracture model, the domain is subdivided into two domains: fracture system and porous media (matrix) and it is assumed that fractures are two-dimensional entities within a three-dimensional space. The system of Eqs. (10) and (5) are solved in these two subdomains (Monteagudo and Firoozabadi, 2004).

⎧ τ f − σc ,& τ > σ ⎪ f c dt , sd =⎨ Ks ⎪ ⎩ 0, & τf < σc

Ω = Ωm +df Ωf

dn, sd =dt , sd tan(Ψ)

(16)

In Eq. (16), Ψ is the shear dilation angle. The above mentioned relations are valid if the fracture is linearly elastic with respect to normal and tangential stresses applied to it. However, if stress levels exceed a certain value, plastic deformation occurs and diﬀerent relation should be used to correlate displacement with the stress ﬁeld. Using the Mohr-Coulomb criterion, fracture failure point can be determined. According to the Mohr-Coulomb criterion for rock media, maximum shear stress that can be applied to fracture, without failure is as follows (Olsson and Barton, 2001).

τf , max =σf tan(ϕ)

(17)

Where τf , max is the maximum shear stress before failure, and ϕ is the fracture internal friction angle. Assuming perfect plastic fracture behavior after failure, the following relations can be written for fracture aperture.

(11)

In Eq. (11) Ω represents the whole domain, m and f denote the matrix and fracture system, respectively, and d f is the fracture aperture. Since the fracture system,Ω f is two-dimensional and Ω m is a three-dimensional domain, in order to resolve the discrepancy between the dimensions, fracture system domain is multiplied by fracture aperture, d f . 2.4. Fracture aperture and mechanical stress As mentioned earlier, fracture system properties are functions of the reservoir mechanical stress ﬁeld. Two distinct mechanisms cause the fracture aperture to change: 1- Variation of normal stress applied to fracture, directly alters the fracture aperture. The change in fracture aperture, Δdn is related to normal stress applied to the fracture, σf according to the well-known Barton and Bandis equation (Barton et al., 1985):

Δdn=

Dn 1+

Kn, i Dn σf

(12)

Where Δdn is the normal fracture closure, Dn is the maximum fracture aperture closure, and Kn, i is the initial normal fracture stiﬀness. Normal fracture stiﬀness is a measure of fracture's sensitivity to normal stress and can be calculated by Eq. (13).

⎡ ⎤−2 σf ⎥ Kn=Kn, i ⎢1− ⎣ Kn, i Dn + σf ⎦

(15)

Where, dt , sd is the tangential displacement associated with shear dilation, and Ks is the shear stiﬀness of fracture. The normal component of displacement caused by shear dilation is (Barton et al., 1985):

(13)

Fig. 2. Solution procedure in each time step.

3

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

d f =d f ,0−Δdn+min(dn, sd , Ucs ) dn, sd ⎧ 0 τf ≤σf ⎪ dt , sd tan(Ψ) ⎪ σ ≤ τ =⎨ ⎛ c f < τf , max ⎪ ⎜ τ f , max − σc ⎞⎟ tan(Ψ)+ σf − σf , max τf =τf , max ⎪ ⎝ Ks ⎠ Kn tan(Ψ) ⎩

(18)

Where σf , max is the normal stress corresponding to maximum shear stress, and Ucs is the maximum shear displacement. After fracture aperture is calculated, fracture permeability can also be calculated by Eq. (19).

kf =

d f2 (19)

12

Fig. 3. 2D validation problem.

3. Numerical method and solution strategy

Table 1 input data.

200m × 400m

In this section, the numerical method used to solve the system of Eqs. (5) and (10) is presented. Furthermore, solution strategy and procedure are also explained.

5 × 10−14m2 0.15

3.1. Galerkin Finite Element method

Rock properties Domain dimensions Matrix permeability, km Matrix Porosity, ϕ Compressibility, cr

6.29 × 10−10Pa−1 0.83

Biot’s coeﬃcient, α Young’s modulus, E Poisson’s ratio, ν Internal friction angle, ϕ Shear dilation angle, Ψ Shear strength, σc

5.8 × 109Pa 0.3 24.9° 5

Initial normal stiﬀness, Kn, i

1.2 × 1011Pa /m

Shear stiﬀness, Ks

1.2 × 1011Pa /m

No-load fracture aperture, df , 0

2.5 × 10−4m

Maximum fracture aperture closure, Dn

2.4 × 10−4m

Fluid properties Viscosity, μ

1.8 × 10−4Pa. s

Compressibility, cf

5.39 × 10−9Pa−1

Initial and Boundary conditions Initial pore pressure, Pi

10 × 106Pa

Injecting well pressure, Pinj

20 × 106Pa

Producing well pressure, Pprod

10 × 106Pa

In-situ stress, σis

20 × 106Pa

In this work, we use the Galerkin variational method with Finite Element discretization to solve the system of Eqs. (5) and (10). Numerical methods based on the Finite Element scheme are numerically more demanding compared to other numerical methods, such as the Finite Diﬀerence method, but since they can be implemented in conjunction with unstructured gridding, Finite Element methods are suitable for three-dimensional problems with complex geometry. Another advantage of the Finite Element based methods is that mesh reﬁnement can be restricted to speciﬁc areas that require more precise grid structure, such as the fractures network and wellbore zone (Karimi-Fard and Firoozabadi, 2003). According to the Galerkin variational method, unknown variables, in this case pore pressure, P and displacement vector, u are approximated based on the nodal values in each element.

5 × 106Pa

n

P= ∑ Pj Nj u = j =1

n

∑ ui,j Nj i = x, y, z

Where Nj is the shape function of j

Fig. 4. Numerical solution for model validation.

4

(20)

j =1 th

node. ux , u y and uz are

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 5. Fracture distribution for Section 5.1.

Table 2 input data for Section 5.1. Rock properties Domain dimensions No-load fracture aperture, d f ,0

1 × 10−3m

Maximum fracture aperture closure, Dn

0.98 × 10−3m

{ε} = [B]{u}e

(29)

{σ} = [E ]{ϵ}

(30)

After calculating the above-mentioned matrices for each element, they are combined to form the global system of equations.

200m × 200m × 100m

Ne, f

Ne, m

[M ] =

∑

[T ]em +d f

e =1

∑ [T ]ef e =1

Ne, f

Ne, m

Table 3 Fracture properties for different rock types.

[M ] =

∑

[M ]em +d f

e =1

Fracture properties

Case 1

Case 2

Case 3

Initial normal stiﬀness, Kn, i [GPa /m] Shear stiﬀness, Ks [GPa /m]

1200 600

120 60

12 6

(31)

∑ [M ]ef e =1

(32)

Ne, m

[K ] =

∑

[K ]e (33)

e =1 Ne, m

[Q ] =

displacements in x , y and z direction, respectively. Multiplying Eqs. (5) and (10) by a weight function and integrating over the domain of each element, results in Eqs. (21) and (22). In the Galerkin method this weight function is the shape function itself.

⎡⎛

⎡

⎤

∫Ω ⎢⎣α∇. ∂∂ut +S ∂∂Pt +qsc′ ⎥⎦ Nd Ω

=

[K ]{u}−[Q]{P}={ f } (21)

[M ]

[M ]

(23)

d {P} d{u} + [Q]T + [T ]{P} = {qsc} dt dt

e

(24)

Where [K ], [Q], [M ] and [T ] are coeﬃcient matrices, that are functions of element geometry, ﬂuid characteristics, and mechanical properties of both rock and fractures:

[K ]e =

∫Ω [B]T [E ][B] d Ω

(25)

[Q]e =

∫Ω α [B]T m [N ] d Ω

(26)

[M ]e =

∫Ω S [N ]T [N ] d Ω

(27)

[T ]e =

∫Ω kμ [∇N ]T [∇N ] d Ω

Where [m] = {1 1 1 0 0

0}T

d {P} d{u} +[Q]T +[T ]{P}={qsc} dt dt

(35) (36)

Where {P} and {u} are the global nodal values of pore pressure and displacement vector, respectively. Several approaches can be taken to solve Eqs. (35) and (36). One approach is to solve both system of equation simultaneously using a fully implicit scheme (Garipov et al., 2016). Since the equations are coupled together, this method will yield the most accurate result. However, because fracture properties are time-dependent, the coeﬃcient matrices, must be updated in each time step. If we denote the number of total nodes by n , this means that in every time step a system of 4n must be solved simultaneously, in order to ﬁnd the nodal values of pore pressure and displacement vector. Another approach is to ﬁrst solve one system of equations and then solve the other set of equations based on the values acquired from the ﬁrst system of equations. Since the time-dependency of coeﬃcients has its roots in the ﬂuid ﬂow equation, this method will allow us to only change the coeﬃcient matrices of the ﬂuid ﬂow equation. This in turn, will reduce the computational cost signiﬁcantly. The method introduced in this work is similar to the latter approach. As shown by Fig. 2, in each time step, ﬁrst the ﬂuid ﬂow equation is solved implicitly to ﬁnd pore pressure, P . Then, based on the calculated values for pore pressure, P , the second set of equations is solved to ﬁnd the nodal values of the displacement vector, u. After calculating the displacement vector, the stress ﬁeld and then new fracture apertures are calculated using Eq. (18). Iteration is done within each time step until the unknown variables, P and u, converge.

(22)

If displacement vector and pore pressure in above equations, are replaced by the approximations of Eq. (20), they can be written in matrix form:

[K ]{u} − [Q]{P} = { f } e

(34)

Subscripts m and f denote porous media and fractures, respectively. The global system of equations is:

⎞

⎛

∫Ω ∇.⎜⎝ kμ ∇P⎟⎠ Nd Ω

[Q]e

e =1

⎤

⎞

∫Ω ⎢⎣ ⎜⎝K + 13 G⎟⎠ ∇ϵ+G∇2u−α∇P⎥⎦ Nd Ω= ∫Ω fNd Ω

∑

(28) . Also, by deﬁnition: 5

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 6. Steady-state pressure distribution for diﬀerent types of rock.

4. Model validation

The same process is repeated for next time steps.

In this section, in order to validate the accuracy of the presented method, the numerical solution of a two-dimensional problem is compared to the solution obtained by Gu et al. (2014). First, the problem is described and then the numerical solutions are presented. The problem consists of a two-dimensional domain with two horizontal wells at each side, as shown by Fig. 3. Four fractures are located inside the domain according to Fig. 3. Injecting and producing wells, each operate at constant pressure. Pore pressure is uniform at t = 0 and in-situ stress applied to the outer boundaries of the domain is constant at all times. Also, the fractures have similar properties at t = 0 . The input data used in this section is listed in Table 1. Fig. 4 shows the numerical solution for the fracture apertures along the width of the domain in steady-state conditions. Comparison between the results obtained in the present work and the work done by Gu et al., 2014. demonstrates that the numerical method is valid and accurate. It can be seen in Fig. 4 that the variation in fracture aperture

3.2. Computational grid Acquiring accurate solutions in fractured porous media requires an appropriate grid structure. Geological features such as fractures, faults and well locations must be considered carefully, in order to create such an appropriate grid. This involves adjusting the grid resolution near such features. Optimum grid resolution can improve solution stability and accuracy (Bahrainian and Daneh Dezfuli, 2014). In the present work, Gmsh software has been used to generate three-dimensional unstructured grid (Geuzaine and Remacle, 2009). Gmsh software implements a Delaunay-based triangulation scheme for grid generation. Here, an adaptive unstructured grid, using tetrahedral elements. Since fractured porous media have complex geometries, unstructured grid is more suitable for such formations. 6

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 7. Fracture aperture distribution for diﬀerent types of rock.

introduced and the validity of its solution has been established. In this section, two problems are investigated by this method, in order to demonstrate the capabilities of the presented model.

along the width of the domain is almost a straight line with constant slope. This is due to the fact that in steady-state conditions, pressure gradient between the injecting and the producing well is approximately constant, which means that pore pressure variation is also a line with constant slope. Therefore, according to the governing equations, the change in fracture aperture would also be a line with constant slope.

5.1. Fractured reservoir with a single production well Fig. 5 shows the fracture distribution inside the reservoir. The reservoir is 200 m×200 m and has a formation thickness of 100 m. Fractures no. 1 to no. 7 are connected to each other and fracture no. 5 is a hydraulic fracture that intersects the production well. The produc-

5. Application in reservoir simulation So far, a coupled ﬂuid ﬂow and geomechanics model has been

Fig. 8. Left: Volumetric and Right: Cumulative production rate vs. time for diﬀerent rock types.

7

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

than shear dilation eﬀects. Fig. 8 shows the production rate and cumulative production variations with time, respectively. It is evident that approximately at day 4 ﬂuid ﬂow reaches steady-state and after that production rate remains constant and cumulative production is a line with constant slope. In case 1, in which fractures are most rigid, there is very little variation in production rate compared to the non-sensitive case. But in cases 2 and 3, reduction in production rate. The reason for this rate reduction is that as pore pressure of the reservoir decreases, the aperture of the fractures also decrease, making it harder for the ﬂuid to travel through the fracture network. Since, pressure gradient between well and reservoir is similar to non-sensitive conditions, more resistance against ﬂuid ﬂow result in a decrease in the production rate.

Fig. 9. Schematic of reservoir in Section 5.2.

5.2. Two-layer fractured reservoir with production and injection wells

Table 4 Input data for Section 5.2.

In this case, we consider a fracture reservoir that consists of two geological layers. The lower layer is 100 times more permeable than the upper shale layer. As shown in Fig. 9, there are two horizontal wells located at opposite sides of the reservoir. The injection well is placed inside the lower layer and production well is located in the upper shale layer. Fracture network geometry and distribution is similar to case 5.1, and ﬂuid and rock properties are given in Table 4. Outer boundaries of the reservoir are impermeable and no-ﬂow boundary condition is imposed. Initial reservoir pore pressure is 2 × 107Pa , production well pressure is 107Pa , and injection well pressure is equal to 3 × 107Pa . Similar to Section 5.1, the outer boundaries do not move and the displacement vector is equal to zero at the outer boundaries. Fluid properties are similar to Section 4 and rock properties are given in Table 2. Other rock properties are similar to Section 4. In-situ stress is constant and equal to 5 × 107Pa . Fig. 10 displays the pressure distribution at 4 diﬀerent times in a section of reservoir that intersects both wells. At t = 1hr there is still no signiﬁcant amount of pressure diﬀusion and ﬂuid ﬂow has not yet reached to the fracture system. Once we get to t = 1day , it is evident from the ﬁgure that ﬂuid is being transported via the fracture system from the injection well to the production well. At this stage, the matrix does not have a signiﬁcant contribution in producing the ﬂuid, but since the lower layer is more permeable, it plays a role in the injection process. At t = 10days the fracture system is almost completely pressurized by the injection well and its pressure is equal to the injection well. Here, it can be seen that matrix is starting to contribute in the production of ﬂuid. When the ﬂow reaches steady-state conditions, almost all of the reservoir has the same pressure as the injection well. The reason is that the lower more permeable layer is much more capable of transporting the ﬂuid, but since the upper shale layer is relatively impermeable, ﬂuid starts to build up inside the reservoir. Therefore the reservoir pore pressure also builds up until it reaches the injection well pressure. Changes in injection and production rates with time are displayed in Fig. 11. Four distinct time spans can be distinguished on the injection and production rate curves. It can be seen that since the injection well is situated in the more permeable layer, the injection rate is bigger than the production rate at all times. In the ﬁrst time span (I), both rates are rapidly decreasing, due to the high pressure gradients around the wellbores at early times. In the second stage (II), production rate stays unchanged. The reason for this is that the decrease in reservoir pore pressure around the production rate is compensated by the ﬂow coming from the injection rate and traveling through the fracture network towards the production well. Therefore, pressure gradient near the production well remains constant and as a result production rate also does not change. As ﬂow reaches the third stage (III), the increase in pore pressure due to injection surpasses the decrease in pore pressure due to production in near production well zone. This leads to an increase in ﬂow rate in the production well. In

Rock properties Domain dimensions Lower layer permeability, km, l

200m × 200m × 100m

5 × 10−14m2

Upper layer permeability, km, u

5 × 10−16m2

No-load fracture aperture, df , 0

1 × 10−3m

Maximum fracture aperture closure, Dn

0.98 × 10−3m

Fluid properties Viscosity, μ

1.8 × 10−3Pa. s

Initial and Boundary conditions Initial pore pressure, Pi

20 × 106Pa

Injecting well pressure, Pinj

30 × 106Pa

Producing well pressure, Pprod

10 × 106Pa

In-situ stress, σis

20 × 106Pa

tion well is positioned at (x, y)=(100m, 100m ) and is perforated from z = −25m to z = 25m . Initial reservoir pore pressure is 5 × 107Pa , the production well pressure is 107Pa , and pore pressure at the outer boundaries of the reservoir is 5 × 107Pa . It is assumed that outer boundaries of the reservoir do not move, thus the displacement vector is equal to zero at the outer boundaries. Fluid properties are similar to Section 4 and rock properties are given in Table 2. Other rock properties are similar to Section 4. In-situ stress is constant everywhere and equal to 5 × 107Pa . We consider three types of rock with three diﬀerent values for normal and shear fracture stiﬀness, which are given in Table 3. Fig. 6 shows the comparison of steady-state pressure distribution for nonsensitive and sensitive reservoir. It can be seen that as we decrease the initial normal fracture stiﬀness, Kn, i and shear stiﬀness, Ks , pressure diﬀusion also decreases. This is due to the fact that with lower stiﬀness, fracture are more easily aﬀected by the stress ﬁeld, and therefore for lower values of Kn, i , fracture aperture reduction is more dramatic. This reduction in fracture aperture is caused by the production process. As the production process goes on, reservoir pore pressure begins to deplete, changing the stress distribution within the reservoir. This change in stress ﬁeld leads to the reduction in fracture aperture. Conversely, for injection processes, fractures might open up, allowing the ﬂuid ﬂow to travel through the fracture system with less diﬃculty. Reservoir pressure depletion causes normal closure of fracture aperture and for less values of normal fracture stiﬀness, this closure becomes greater. However, as discussed in previous sections, shear dilation acts against this mechanism. In other words, lower shear stiﬀness makes it easier for the fracture planes to slide on one another. This means that for lower values of Ks , shear dilation eﬀect is more signiﬁcant. But as we see in Figs. 6 and 7, for less rigid rocks, fracture aperture closes more. The reason behind this is that, in this case, direct normal displacement of fractures due to normal stress is much greater 8

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 10. Pressure distribution for Section 5.2.

the fourth stage (IV), ﬂuid ﬂow becomes steady and injection rate equals the production rate. Fig. 12 demonstrates the variations in average fracture aperture with time for every fracture. Average fracture aperture for all fractures rapidly increases until t = 100days , and after that it remains roughly the same. The increase in fracture aperture is due the pressure build up during the injection process. Since injection rate is always greater than the production rate, pore pressure of the reservoir increases, which as a result, leads to an increase in fracture aperture.

along with an adaptive unstructured grid was used to numerically solve the coupled problem. An experimental model was incorporated to relate the fracture aperture with the reservoir stress ﬁeld. The accuracy of the presented model has been validated by previously published results. The method has been used to solve two problems regarding the ﬂuid ﬂow in fractured porous media. The following points can be inferred from the numerical results: 1. Fracture aperture is dependent on the applied stresses to the fracture surfaces. But, the degree of this dependency is highly related to mechanical properties of rock and the production/injection conditions. 2. Normal displacement of the fracture aperture due to normal stress is much more signiﬁcant than normal displacement due to shear dilation eﬀects. 3. Production rate in fracture reservoirs is very sensitive to changes in fracture aperture and permeability. Therefore, the geomechanical

6. Conclusion In this work, a coupled ﬂuid ﬂow and geomechanics numerical method based on the discrete-fracture model and Finite Element method is introduced. This method, utilizes a sequential implicit scheme for solving the poroelasticity and ﬂuid ﬂow (continuity) equations. The Galerkin method with Finite Element discretization 9

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al.

Fig. 11. Left: injection and production rates vs. time, right: cumulative injection and production vs. time. Gutierrez, M., Lewis, R.W., Masters, I., 2001. Petroleum reservoir simulation coupling ﬂuid ﬂow and geomechanics. Soc. Pet. Eng. 4, 164–172. Gu, S., Liu, Y., Chen, Z., 2014. Numerical study of dynamic fracture aperture during production of pressure-sensitive reservoirs. Int. J. Rock. Mech. Min. Sci. 70, 229–239. Hu, L., Winterﬁeld, P.H., Fakcharoenphol, P., Wu, Y., 2013. A novel fully-coupled ﬂow and geomechanics model in enhanced geothermal reservoirs. J. Pet. Sci. Eng. 107, 1–11. Jaﬀre, J., Mnejja, M., Roberts, J.E., 2011. A discrete fracture model for two-phase ﬂow with matrix-fracture interaction. Procedia Comput. Sci. 4, 967–973. Jalali, M.R., Dusseault, M.B., 2012. Coupling geomechanics and transport in naturally fractured reservoirs. Int. J. Min. Geo-Eng. (IJMGE) 46 (2), 105–131. Karimi-Fard, M., Durlofsky, L., Aziz, K., 2004. An eﬃcient discrete fracture model applicable for general-purpose reservoir simulators. SPE J. 9 (2), 227–236. Karimi-Fard, M., Firoozabadi, A., 2003. Numerical simulation of water injection in fractured media using the discrete-fracture model and the Galerkin method. SPE Reserv. Eval. Eng. 6 (2), 117–126. Khoei, A.R., Shamloo, A., Anahid, M., Shahim, K., 2006. The extended ﬁnite element method (X-FEM) for powder forming problems. J. Mater. Process. Technol. 177 (1– 3), 53–57. Kilmer, N.H., Morrow, N.R., Pitman, J.K., 1987. Pressure sensitivity of low permeability sandstones. J. Pet. Sci. Eng. 1 (1), 65–81. Kim, J., Moridis, G., Yang, D., Rutqvist, J., 2012. Numerical studies on two-way coupled ﬂuid ﬂow and geomechanics in hydrate deposits. Soc. Pet. Eng. 17 (2), 485–501. Lamb, A.R., Gorman, G.J., Elsworth, D., 2013. A fracture mapping and extended ﬁnite element scheme for coupled deformation and ﬂuid ﬂow in fractured porous media. Int. J. Numer. Anal. Methods Geomech. 37 (17), 2916–2936. Martin, V., Jaﬀre, J., Roberts, J., 2005. Modeling fractures and barriers as interfaces for ﬂow in porous media. SIAM J. Sci. Comput. 26 (5), 1667–1691. McClure, M.W., Babazadeh, M., Shiozawa, S., Huang, J., 2016. Fully coupled hydromechanical simulation of hydraulic fracturing in 3d discrete-fracture networks. Soc. Pet. Eng. 21, 1302–1320. Minkoﬀ, S.E., Stone, C.M., Bryant, S., Peszynska, M., Wheeler, M.F., 2003. Coupled ﬂuid ﬂow and geomechanical deformation modeling. J. Pet. Sci. Eng. 38, 37–56. Moinfar, A., Sepehrnoori, K., Johns, R., Varavei, A., 2013. Coupled geomechanics and ﬂow simulation for an embedded discrete fracture model. In proceedings of: Society of Petroleum Engineers – SPE Reservoir Simulation Symposium 2, 1238–1250. Monteagudo, J.E.P., Firoozabadi, A., 2004. Control-volume method for numerical simulation of two-phase immiscible ﬂow in two- and three-dimensional discretefractured media. Water Resour. Res. 40, 1–20. Norbeck, J., Huang, H., Podgorney, R., Horne, R., 2014. An integrated discrete fracture model for description of dynamic behavior in fractured reservoirs. In proceedings of: 39th workshop on geothermal reservoir engineering, Stanford university, California. Olsson, R., Barton, N., 2001. An improved model for hydromechanical coupling during shearing of rock joints. Int. J. Rock. Mech. Min. Sci. 38 (3), 317–329. Pedrosa, O.A., 1986. Pressure transient response in stress-sensitive formations. In: Proceedings of the SPE California Regional Meeting, Oakland, California, USA, 203214. Philip, Z.G., Jennings, J.W., Olson, J.E., Laubach, S.E., Holder, J., 2005. Modeling coupled. Soc. Pet. Eng. 8 (4), 300–309. QanbariF., ClarksonC.R., 2012. Rate transient analysis of stress-sensitive formations during transient linear ﬂow period. In: Proceedings of the SPE Canadian Unconventional Resources Conference, Calgary, Alberta, Canada. Qanbari, F., Clarkson, C.R., 2014. Analysis of transient linear ﬂow in stress-sensitive formations. SPE Reserv. Eval. Eng. 17 (1), 98–104. Rutqvist, J., Stephansson, O., 2003. The role of hydromechanical coupling in fractured rock engineering. Hydrogeol. J. 11 (1), 7–40. Segura, J.M., Paz, C.M., de Bayser, M., Zhang, J., Bryant, P.W., González, N.A., Rodrigues, E., Vargas, P.E., Carol, I., Lakshmikantha, M.R., Das, K.C., Sandha, S.S.,

Fig. 12. Fracture aperture vs. time.

eﬀects on fracture network should be taken into account for accurate prediction of well performance. References Bahrainian, S.S., Daneh Dezfuli, A., 2014. A geometry-based adaptive unstructured grid generation algorithm for complex geological media. Comput. Geosci. 68, 31–37. Barton, N., Bandis, S., Bakhtar, K., 1985. Strength, deformation and conductivity coupling of rock joints. Int. J. Rock. Mech. Min. Sci. Geomech. Abstr. 22 (3), 121–140. Biot, M.A., 1941. General theory of three‐dimensional consolidation. J. Appl. Phys. 12, 155–164. Castelletto, N., White, J.A., Tchelepi, H.A., 2015. Accuracy and convergence properties of the ﬁxed-stress iterative solution of two-way coupled poromechanics. Int. J. Numer. Anal. Methods Geomech. 39, 1593–1618. Chen, H.-Y., Teufel, L., 2000. Coupling ﬂuid-ﬂow and geomechanics in dual-porosity modeling of naturally fractured reservoirs – model description and comparison. In: Proceedings of SPE International Petroleum Conference and Exhibition in Mexico, Villahermosa, Mexico Chen, H.-Y., Teufel, L.W., Lee, R.L., 1995. Coupled ﬂuid ﬂow and geomechanics in reservoir study - i. theory and governing equations. In: proceedings of: SPE Annual Technical Conference and Exhibition, Dallas, Texas, USA Chin, L.Y., Raghavan, R., Thomas, L.K., 2000. Fully coupled analysis of well responses in stress-sensitive reservoirs. SPE Reserv. Eval. Eng. 3, 435–443. Garipov, T.T., Karimi-Fard, M., Tchelepi, H.A., 2016. Discrete fracture model for coupled ﬂow and geomechanics. Comput. Geosci. 20 (1), 149–160. Geuzaine, C., Remacle, J.-F., 2009. Gmsh: a three-dimensional ﬁnite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Methods Eng. 79 (11), 1309–1331.

10

Journal of Petroleum Science and Engineering (xxxx) xxxx–xxxx

M. Moradi et al. Cerqueira, R., Mello, U., 2016. Coupling a ﬂuid ﬂow simulation with a geomechanical model of a fractured reservoir. American Rock Mechanics Association. In: Proceedings of the 50th US Rock Mechanics/Geomechanics Symposium, Houston, Texas, USA, 26-29. Settari, A., Mourits, F.M., 1998. A coupled reservoir and geomechanical simulation system. Soc. Pet. Eng. 3 (3), 219–226. Shamloo, A., Azami, A.R., Khoei, A.R., 2005. Modeling of pressure-sensitive materials using a cap plasticity theory in extended ﬁnite element method. J. Mater. Process. Technol. 164–165, 1248–1257. Shaoul, J.R., Ayush, A., Park, J., Pater, C.J.D., 2015. The eﬀect of stress sensitive permeability reduction on the evaluation of post-fracture welltests in tight gas and

unconventional reservoirs. In: proceedings of the SPE European damage conference and exhibition, Budapest, Hungary. Wang, H.-X., Wang, G., Chen, Z.-X.J., Wong, R.C.K., 2010. Deformational characteristics of rock in low permeable reservoir and their eﬀect on permeability. J. Pet. Sci. Eng. 75 (1–2), 240–243. Zhang, M.Y., Ambastha, A.K., 1994. New insights in pressure-transient analysis for stress-sensitive reservoirs. In: Proceedings of the 69th SPE Annual Technical Conference and Exhibition, New Orleans, L.A., USA, pp. 617–628. Zienkiewicz, O., Taylor, R., 2005. The Finite Element Method for Solid and Structural Mechanics. Elsevier.

11

Copyright © 2022 COEK.INFO. All rights reserved.