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.
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.
[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:
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:
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 inkeilegavlen2019porepy , 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 .
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.
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.
- (1) Aavatsmark, I. (2002) An Introduction to Multipoint Flux Approximations for Quadrilateral Grids. Computational Geosciences
- (2) Berge, R. L., Berre, I., Keilegavlen, E., Nordbotten, J. M. and Wohlmuth, B. (2019) Finite volume discretization for poroelastic media with fractures modeled by contact mechanics. International Journal for Numerical Methods in Engineering
- (3) Salimzadeh, S., Paluszny, A., Nick, Hamidreza M. and Zimmerman, R. W. (2018) A three-dimensional coupled thermo-hydro-mechanical model for deformable fractured geothermal systems. Geothermics
- (4) Keilegavlen, E., Berge, R., Fumagalli, A., Starnoni, M., Stefansson, I., Varela, J. and Berre, I (2019) PorePy: An Open-Source Software for Simulation of Multiphysics Processes in Fractured Porous Media. arXiv preprint arXiv:1908.09869
- (5) Geuzaine, C., Remacle, J.-F. (2009) Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering
- (6) Nordbotten, J. M. (2016) Stable Cell-Centered Finite Volume Discretization for Biot Equations. SIAM Journal on Numerical Analysis
- (7) Hüeber S, Stadler G, Wohlmuth B.I. (2018) A primal-dual active set algorithm for three-dimensional contact problems with coulomb friction. SIAM J Sci Comput
- (8) Run scripts at https://github.com/IvarStefansson/Finite-volume-discretisation-of-fracture-deformation-in-thermo-poroelastic-media.git