Finite volume discretisation of fracture deformation in thermo-poroelastic media

by   Ivar Stefansson, et al.
University of Bergen

This paper presents a model where thermo-hydro-mechanical processes are coupled to a deformation model for preexisting fractures. The model is formulated within a discrete-fracture-matrix framework where the rock matrix and the fractures are considered as individual subdomains, and interaction between them takes place on the matrix-fracture interfaces. A finite volume discretisation implemented in the simulation toolbox PorePy is presented and applied in a simulation showcasing the effects of the different mechanisms on fracture deformation governed by contact mechanics, as well as their different timescales.



There are no comments yet.


page 7


A fully coupled numerical model of thermo-hydro-mechanical processes and fracture contact mechanics in porous media

A range of phenomena in the subsurface is characterised by the interplay...

Learning crystal plasticity using digital image correlation: Examples from discrete dislocation dynamics

Digital image correlation (DIC) is a well-established, non-invasive tech...

Numerical modelling of convection-driven cooling, deformation and fracturing of thermo-poroelastic media

Convection-driven cooling in porous media influences thermo-poro-mechani...

Polyconvex anisotropic hyperelasticity with neural networks

In the present work, two machine learning based constitutive models for ...

On the stability of the generalized, finite deformation correspondence model of peridynamics

A class of peridynamic material models known as constitutive corresponde...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

We consider the simulation of fully coupled thermo-hydro-mechanical (THM) dynamics in fractured porous media, where the fractures can undergo sliding if the shear forces on the fracture planes are sufficient to overcome frictional resistance. These processes are highly relevant for several subsurface applications, including geothermal energy extraction, storage of CO

and energy and groundwater management. Our simulation approach is based on three main ingredients: First, conservation of mass, energy and momentum is preserved under discretisation by the employment of a fully coupled finite volume approach for the governing equations. Second, the network of fractures, which act as main conduits for fluid flow and energy transport, is explicitly represented in the simulation model. Specifically, the fractures are represented as lower-dimensional manifolds that are embedded in the host porous medium, thus the simulation model is defined on a mixed-dimensional geometry. Third, the sliding of fractures is modelled as a frictional contact problem, which is solved by an active set approach. The discretisation of the contact problem benefits from the finite volume approach, which directly provides discrete representations of displacements as well as of mechanical, fluid and thermal forces on the fracture surfaces. Furthermore, the explicit degrees of freedom for fluid pressure inside the fractures allow us to capture the critical interplay between elevated fluid pressures and fracture deformation.

2 Model

We consider a mixed-dimensional geometry which is decomposed in subdomains of different dimensions representing the host porous medium and the lower-dimensional planar fractures, and separated by interfaces, see keilegavlen2019porepy for details. Variables and governing equations are defined on subdomains and interfaces, with full flexibility to vary the type and number of variables and equations between geometric objects. The framework accommodates heterogeneous and multiphysics models, with a natural treatment of modelling and discretisation of mixed-dimensional problems.

We denote a subdomain by [i] and its boundary by [i], and identify the variables defined within it by subscript . Where convenient, we will denote the higher-dimensional matrix domain by [h] and lower-dimensional fracture domains by [l], as indicated in Fig. 1.

Figure 1: Schematic representation of a fracture [l] and a matrix subdomain [h] to the left. The two subdomains are separated by the interface, whose two sides are denoted by and . Also shown are the projection operators used for transfer of variables between the subdomains and the interface. To the right, we show the block structure of the matrix resulting from a fully implicit discretisation of Eqs. 1 through 5.

[i] may be divided into the external boundary and the internal fracture boundary , which coincides geometrically with both the immersed fracture domain [l] and the interface between [h] and [l]. This interface is denoted by [j], and the associated variables identified by subscript . The two sides of the fracture are denoted by and , as shown in Fig. 1. Projection of variables from the interface to the subdomains is performed by hj and , respectively, whereas hj and lj project from the subdomains to the interface.

For [h], the primary variables are displacement , fluid pressure and temperature . The fluid flux and the advective and conductive heat fluxes are denoted by , and , respectively. Assuming all external source and sink terms to be zero, conservation of momentum, mass and energy in [h] can be defined as salimzadeh2018thm


Similarly, conservation of mass and energy in [l] and fracture intersection domains [x] is given by


Here, the specific volume accounts for the extension in the collapsed dimension(s), while subscript 0 denotes the initial value and and indicate fluid and solid parameters, respectively.

Denoting the trace operator by ⋅, the conditions on [j] are the three flux relationships and the traction balance:


Eqs. 1 through 3 are complemented by the internal boundary conditions , , and , and standard Dirichlet and Neumann conditions on the external boundaries.

The fracture deformation is described by relations between the contact traction on the fracture surface, , and the jump in displacement over the fracture. Denoting the displacements on the two sides of the interface by and , the displacement jump is , and denotes its increment. is also related to the aperture and specific volume: for [l], we set . For [x] we use the mean and product of for the adjacent cells of the intersecting fractures when computing and , respectively.

Since the fracture deformation depends on the traction caused by the contact between the two surfaces, we subtract the contribution from the pressure [l] on the fracture surfaces. Thus, the interface tractions on the two fracture surfaces and the traction balance are


Denoting tangential and normal components of vectors on the fracture by subscripts and

, respectively, the fracture deformation for [l] is governed by three non-penetration relations and three friction law constraints:


Further detail on the fracture deformation is found in berge2019hmdfm ; hueeber2008 , whereas the definition of all parameters may be found in the repository at stefansson2019runscripts .

3 Discretisation

Applying Implicit Euler for the temporal discretisation, the scalar conservation Eqs. 1 and 2 are discretised using a Multi-Point Flux Approximation aavatsmark2002mpfa for the conductive terms and a first order upwind scheme for the advective term. The momentuum conservation equation and the terms in the scalar conservation laws are discretised using the FV scheme introduced in nordbotten2016biot . The scheme, termed Multi-Point Stress Approximation (MPSA), is based on local momentum conservation and is formulated in terms of discrete cell centred pressures and displacement unknowns. Originally developed for the pure hydro-mechanical problem, the coupled discretisation approach can readily be extended to the THM case.

Thanks to the structure provided by the mixed-dimensional framework, the discretisation of the coupling fluxes of Eq. 3 consists of two simple tasks. Discrete projection operators transfer variables from higher-dimensional faces and lower-dimensional cells to the interface cells. The interface fluxes are then discretised directly using the projected variables, see keilegavlen2019porepy .

The traction balance and fracture deformation relations of Eqs. 4 and 5 are formulated in terms of displacement and traction on the fracture surfaces. The former is included as a primary interface variable, and thus directly available. While the latter is not a primary variable, the FV framework is formulated in terms of face tractions, implying that in the discrete setting, the surface traction can be reconstructed from the primary variables. Specifically, we apply available discretisation operators to get contributions to the stress from displacements, pressures and temperatures in [h], the interface variable [j], and conditions on external boundaries.

For the discretisation of the fracture deformation relations we first reformulate Eq. 5 as two nonlinear complementary functions and compute their derivatives. A semismooth Newton scheme is applied on the basis of the three sets


which correspond to fracture cells which are open, sticking and sliding, respectively. denotes a numerical parameter and is the friction bound. For each fracture cell , this results in the following cell-wise constraints when computing iterate from the current iterate :


The coefficients , and are functions of , and , and can thus be computed from the previous iterate and time step. For further details of the discretisation and implementation of the fracture deformation equations, we refer to berge2019hmdfm .

In terms of implementation, we mention that the mixed-dimensional framework allows us to discretise each term for each subdomain or interface independently. We may thereby break the highly complex task of discretising the contact conditions with a coupled THM stress down in manageable tasks, and also reuse discretisations for different subdomains. For the global discretisation matrix, this manifests as a two-level block structure. The first level has the subdomains on the diagonal and interfaces on the off-diagonals. The second level corresponds to the primary variables, and has coupling terms between different variables on the off-diagonals.

The model is implemented for two- and three-dimensional problems in the open source simulation toolbox PorePy presented in

keilegavlen2019porepy , and run scripts for the example simulation presented in the following section may be found in the repository stefansson2019runscripts . The simplicial spatial grid is constructed such that the lower-dimensional cells coincide with higher-dimensional faces through a back-end to Gmsh gmsh2009 .

4 Results

To demonstrate the applicability of the model and discretisation, we present simulations of THM and fracture deformation effects for a 2d domain containing seven fractures, see Fig. 2 for the geometry and numbering of the fractures. The setup is designed to expose the method to a wide range of physical driving forces, and thus probe the stability and performance of the simulation model.

Starting out from a homogeneous initial state for all primary variables, the simulation consists of four phases, where we study the effect of sequentially adding different driving forces. In phase I, the deformation is caused by a boundary displacement of applied at the top. To allow the system to reach equilibrium, this phase lasts from to . In phase II, a pressure gradient is applied from left to right. At the end of the phase, at =, the pressure has virtually reached a steady state, see Fig. 2. In phase III, we reduce the temperature at the left boundary of from to and in phase IV we increase it to , thus exploring both thermal expansion and contraction. At the end of each of the two last phases, at and , the domain has has reached a close to uniform temperature.

For the end of each of the simulation phases, Fig. 3 shows the deformation state, i.e. whether and are nonzero for each fracture cell. These show that the state changes for all fractures, and that for each phase, at least three fractures change their state. Fig. 4 shows the norm of and on each fracture throughout the simulation.

Because of the time scale difference, the pressure phase is shown in the left plot and the temperature phase in the right one. The former shows gradual and moderate deformation for the fractures which have nonzero jumps at the end of phase I, and onset of sliding for fracture 3. The latter shows more complex deformation, displaying non-monotone jump evolution for several fractures, e.g. fracture 3 first opening and subsequently closing, and fracture 7 undergoing the reverse process.

Fig. 4 also displays the number of Newton iterations required for convergence for each of the time steps. While fairly stable results are demonstrated, some increase may be observed whenever the fracture state changes more markedly, e.g. when the cooling is introduced at the onset of phase IV.

Figure 2: The equilibrated pressure state at the end of phase II.
Figure 3: The deformation state of the fractures at the end of the four phases.
Figure 4: Norm of displacement jumps at each fracture and number of Newton iterations for each time step. Phase II is shown to the left, with the state at the end of phase I shown at , and phases III and IV to the right. Normal jumps are shown in dashed lines and tangential jumps in solid lines.

5 Conclusion

A FV framework for simulation of thermo-poroelasticity fully coupled to fracture deformation is presented. The contact mechanics problem is naturally discretised in terms of displacement and traction at the fracture faces, exploiting the availability of these in the FV formulation of the THM problem. A numerical example exhibiting complex THM interactions and fracture dynamics demonstrates both the range of the processes captured by the model, and its applicability for challenging problems.