Simulation of hydro-mechanically coupled processes in rough rock fractures using an immersed boundary method and variational transfer operators

TL;DR

This paper introduces an immersed boundary method with variational transfer operators for high-resolution simulation of hydro-mechanical processes in rough fractures, enabling accurate modeling of complex boundary interactions without remeshing.

Contribution

It adapts an immersed boundary approach with variational transfer operators to efficiently simulate fluid-structure interactions in rough fractures, simplifying complex boundary handling.

Findings

Successfully resolves rough fracture boundaries in fluid flow simulations.

Captures flow changes during fracture closure under normal load.

Simulates fracture opening due to increased fluid pressure.

Abstract

Hydro-mechanical processes in rough fractures are highly non-linear and govern productivity and associated risks in a wide range of reservoir engineering problems. To enable high-resolution simulations of hydro-mechanical processes in fractures, we present an adaptation of an immersed boundary method to compute fluid flow between rough fracture surfaces. The solid domain is immersed into the fluid domain and both domains are coupled by means of variational volumetric transfer operators. The transfer operators implicitly resolve the boundary between the solid and the fluid, which simplifies the setup of fracture simulations with complex surfaces. It is possible to choose different formulations and discretization schemes for each subproblem and it is not necessary to remesh the fluid grid. We use benchmark problems and real fracture geometries to demonstrate the following capabilities of the presented approach: (1) Resolving the boundary of the rough fracture surface in the fluid; (2) Capturing fluid flow field changes in a fracture which closes under increasing normal load; and (3) Simulate the opening of a fracture due to increased fluid pressure.

Figures

Equations

$${\rho }_f{\dot{v}}_f+{\rho }_f(v_f\cdot \nabla )v_f-{\mu }_f\nabla \cdot {\sigma }_f(p_f,v_f)$$

$$\nabla \cdot v_f$$

$$v_f$$

$${\sigma }_f\cdot n_f$$

$${\rho }_s{\ddot{u}}_s-\text{div}{\sigma }_s(u_s)$$

$$u_s$$

$${\sigma }_s\cdot n_s$$

$${\dot{u}}_s$$

$${\sigma }_f=-p_f\text{Id}+\frac{{\nabla v_f}^t+\nabla v_f}{2}.$$

$${\mathcal{W}}^v$$

$${\mathcal{W}}^p$$

$$\displaystyle\int_{\Omega_{f}}\rho_{f}\,\big{(}\dot{v}_{f}+(v_{f}\cdot\nabla)\,v_{f}\big{)}\,\cdot\delta v\,dV+\int_{\Omega_{f}}\sigma_{f}:\nabla(\delta v)\,dV$$

$${\int }_{{\Omega }_f}{\rho }_f{f}_{\text{fsi}}\cdot \delta vdV-{\int }_{{\Gamma }_f^N}h\cdot \delta vdA+{\int }_{{\Omega }_f}\nabla v\cdot \delta pdV$$

$$\mathcal{V}^{u}=\big{\{}\delta u\in[H^{1}(\Omega_{s}^{t})]^{d}\,|\,\delta u=0\text{ on }(\Gamma_{s}^{t})^{D}\big{\}},\;d=2,3$$

$${\rho }_s{\int }_{{\Omega }_s^t}{\ddot{u}}_s\cdot \delta udV+{\int }_{{\Omega }_s^t}{\sigma }_s:\nabla (\delta u)dV-{\rho }_s{\int }_{{\Omega }_s^t}{f}_s\cdot \delta udV-{\int }_{({\Gamma }_s^t{)}^N}t\cdot \delta udA=0.$$

$$\displaystyle\int\limits_{\Omega_{f}}\big{(}\Pi(v)-v\big{)}\mu\,d\omega=0\quad\forall\mu\in M.$$

$$\underset{{\Omega }_f}{\int }{w}_h{\mu }_{h}d\omega =\underset{{\Omega }_f}{\int }{v}_h{\mu }_{h}d\omega ,\forall {\mu }_h\in M_h,$$

$$\underset{{\Omega }_f}{\int }\sum _{i=1}^{n_{W_h}}{w}_i{\varphi }_i^{W_h}{\varphi }_k^{M_h}d\omega =\underset{{\Omega }_f}{\int }\sum _{j=1}^{n_{V_h}}{v}_j{\varphi }_j^{V_h}{\varphi }_k^{M_h}d\omega ,k=1,....,n_{M_h}$$

$$\iff \sum _{i=1}^{n_{W_h}}{w}_i\underset{{\Omega }_f}{\int }{\varphi }_i^{W_h}{\varphi }_k^{M_h}d\omega =\sum _{j=1}^{n_{V_h}}{v}_j\underset{{\Omega }_f}{\int }{\varphi }_j^{V_h}{\varphi }_k^{M_h}d\omega ,k=1,....,n_{M_h}.$$

$$\mathbf{w}$$

$$D$$

$$B$$

$$D\mathbf{w}=B\mathbf{v},\text{or}\mathbf{w}=D^{-1}B\mathbf{v}=T\mathbf{v}.$$

$${\mathbf{T}}_{v_s}$$

$${\mathbf{T}}_I$$

$${\mathbf{T}}_{p_f}^{\ast }$$

$${\mathbf{T}}_{v_f}^{\ast }$$

$${\int }_{{\Omega }_f\cap {\Omega }_s}{ϵ}^{-1}\cdot \left({\Pi }_{{\dot{u}}_s}({\dot{u}}_s)-v_f\right)\delta vdw.$$

$$(S(u),\delta u):={\rho }_s{\int }_{{\Omega }_s^t}{\ddot{u}}_s\cdot \delta udV+{\int }_{{\Omega }_s^t}{\sigma }_s:\nabla (\delta u)dV,$$

$$R_s(\delta u)={\rho }_s{\int }_{{\Omega }_s^t}{f}_s\cdot \delta udV-{\int }_{({\Gamma }_s^t{)}^N}t\cdot \delta udA+{\int }_{{\Omega }_f\cap {\Omega }_s}{\sigma }_f:\nabla (\delta u)dV.$$

$$\displaystyle(F(v,p),(\delta v,\delta p))=\rho_{f}\int_{\Omega_{f}}\,\big{(}\dot{v}_{f}+(v_{f}\cdot\nabla)\,v_{f}\big{)}\,\delta v\,dW+\int_{\Omega_{f}}\sigma_{f}:\nabla(\delta v)\,dW$$

$${\int }_{{\Omega }_f}\delta v\nabla vdW$$

$$R_f(\delta v)={\int }_{{\Gamma }_f^N}h\cdot \delta vdA+{\rho }_f{\int }_{{\Omega }_f}{f}_{\text{fsi}}\cdot \delta vdW.$$

$$\left[\begin{array}[]{cc}\mathbf{S}&\mathbf{0}\\ \mathbf{P}_{s}&\mathbf{F}-\mathbf{P}_{f}\end{array}\right]\left[\begin{array}[]{c}\mathbf{y}_{s}\\ \mathbf{y}_{f}\end{array}\right]=\left[\begin{array}[]{c}\mathbf{R}_{s}\\ \mathbf{R}_{f}\end{array}\right],$$

$$\mathbf{P}_{s}:=\epsilon^{-1}\left[\begin{array}[]{c}\mathbf{I}_{\Omega_{f}\cap\Omega_{s}}\cdot\mathbf{M}_{f}\cdot\mathbf{T}_{v_{s}}\\ \mathbf{0}\end{array}\right],$$

$$\mathbf{P}_{f}:=\epsilon^{-1}\left[\begin{array}[]{cc}\mathbf{I}_{\Omega_{f}\cap\Omega_{s}}\cdot\mathbf{M}_{f}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}\end{array}\right].$$

$$(\mathbf{I}_{\Omega_{f}\cap\Omega_{s}})_{k_{l}}:=\left\{\begin{array}[]{ll}1,\;\text{if}\,k=l,\text{and}\,T_{I}(i)>t,\\ 0,\,\text{otherwise}\end{array}\right.$$

Full text

Keywords:** Fluid flow Fracture mechanics Non-matching meshes pseudo-$L^{2}$-projection Immersed boundary Hydro-mechanical coupling Geothermal energy

1 Introduction

Hydro-mechanical (HM) processes in discontinuities such as rock fractures play an important role in a number of reservoir engineering applications, for example enhanced geothermal systems, oil- and gas recovery or CO2 sequestration. They govern productivity and thereby economical feasibility of reservoir operations [41, 26, 1]. Global reservoir permeability is governed by the hydraulic properties of single fractures as well as their degree of connectedness and geometric configuration in a fracture network [7, 56, 38, 10, 17]. The hydraulic properties of a single fracture are subject to constant change by various physical processes, which are often tightly coupled and highly non-linear. Given the impact of fracture hydraulics on resulting flow rates, HM processes in fractures are an actively researched topic. In the following text we exclusively focus on HM processes in single fractures.

For the mechanical behavior of a fracture, the highly variable fracture surface topographies yield a complex distribution of contact area and a non-linear mechanical response to normal or shear loading [3, 36, 19, 25, 49, 40]. Mechanical loading normal to the fracture surface results in convergent fracture closure behavior with increasing load, as more and more area of the fracture comes into contact and local compressive stresses at contact areas increase [3, 37, 25, 54, 44, 47, 46, 22].

For the hydraulic behavior of fractures, rough fracture surface topographies result in highly heterogeneous aperture fields, which characterize the mechanical gap between two fracture surfaces [56, 55]. These heterogeneous aperture fields are naturally subject to change with mechanical loading and yield complex fluid flow patterns through fractures [42, 6, 43, 49, 28, 46]. While an increase in mechanical loading leads to fracture closure, this can also result in elevated fluid pressures, as larger fluid pressure gradients are required to sustain fluid flow rates through closing fractures [44]. If the fluid pressure in the fracture increases sufficiently, it can counteract compressive mechanical load and even cause fractures to open.

Given the small scale surface topography and aperture variations, as well as the tightly coupled physical processes, simulation of HM phenomena can become computationally very expensive. Therefore, most numerical and analytical studies rely on simplifying assumptions when studying HM processes [4, 28, 27, 11, 14, 12]. Commonly used assumptions are Darcy or lubricated flow, global or local parallel plate models for fluid flow in rock fractures, where averaged hydraulic and mechanical parameters are used globally or locally (i.e., on smaller surface patches) to capture the non-linear behavior of HM processes in rock fractures. The number of studies that model HM processes on real rough rock geometries, however, is small. The aim of this study is to present an approach to model HM processes in complex three dimensional geometries without the need to simplify them or the physical relations.

To this end, we use here ideas of the immersed boundary (IB) method which can be considered as a representative of the wider class of fictitious domain (FD) methods. The main feature of FD methods is to substitute the solution of a PDE on a complex domain by embedding or immersing it in a simpler domain on which it is easier to solve the problem. It is outside the scope of this article to give an in-depth review of the intertwined history of FD and IB methods. To this end we refer the reader to [18, 39, 15, 2, 33, 5, 24, 16] and the references therein.

An FD method with fluid structure interaction (FSI) to simulate particulate flow was introduced by [15]. There, the interaction is enforced with Lagrange multipliers distributed over each particle to match the velocities of the fluid and the moving rigid bodies. This method was later extended by [53] to deal with flexible bodies. [2] proposed a fictitious domain/mortar element method to impose the continuity of solid and fluid velocities along the solid-fluid boundary, which was later generalized by [16] and [5] by enforcing the velocity constraint not only over the boundary, but over the whole overlapping region. In the IB method, the immersed structure is represented by a force density term in the Navier-Stokes equations. The solid and the fluid problem are solved on different, possibly non-conforming grids, whereby a Lagrangian formulation for the solid problem and an Eulerian formulation for the fluid problem is used. With respect to the methods based on moving meshes such as the Arbitrary Lagrange Eulerian method [31], no remeshing of the fluid grid is necessary and in addition, different discretization schemes may be chosen for the fluid and the solid problem (see e.g. [8, 30]).

Here, we adapt an immersed boundary approach used in cardiac simulations by [29]. The fluid and solid problem are solved on separate domains, whereby the solid problem is formulated in Lagrangian and the fluid problem in Eulerian coordinates. The equality of fluid and solid velocities is enforced by a penalty method. The complete problem is solved using a staggered approach with a fixed point iteration, solving the fluid and solid problem in sequence (Fig. 1). For the coupling of the fluid and solid problem we use variational transfer operators, more precisely pseudo-$L^{2}$-projections. They are defined as projections which satisfy an orthogonality condition in the sense of the $L^{2}$-scalar product and ensure a stable mapping between the non-matching meshes [9, 23].

One motivation behind the development of this type of methods was to avoid the remeshing of a moving fluid mesh in case of large solid displacements. This is necessary in other FSI approaches to ensure numerical stability of the fluid part of the simulation [16]. In addition, no explicit representation of the fluid-solid boundary is required, which is a particular advantage for high resolution HM simulations. There, small scale changes in the aperture field result in large changes of the fluid flow patterns. The ability to implicitly reflect these changes in the fluid-solid boundary after every deformation, however small, immensely simplifies the setup of these simulations. Lastly, as remarked earlier, the immersed boundary approach gives some freedom which formulation, mesh type or method, can be used for solving each subproblem. These characteristics make this FSI formulation a highly suitable method to compute hydro-mechanically coupled processes in rough fractures in a coupled manner while simultaneously resolving small scale roughness on the fracture surfaces.

This paper is organized as follows. In section 2 we give an overview of the method and the governing equations. We introduce the solid and fluid problem in their strong and weak formulation, followed by the definition and construction of the transfer operators, the discretization in space and time, and finally, an explanation of the coupling and an overview of the algorithm. In section 3 we show numerical examples of the approach. First, we validate the method with 2D and 3D benchmarking problems and then we present numerical examples with real rough surface geometries.

2 Method & governing equations

We use an adapted immersed boundary approach following [16, 30, 29, 32]. In our approach, the FSI problem is solved separately for the fluid and the solid in a fixed point iteration using finite elements. We start by solving the fluid problem formulated by the Navier-Stokes equations for incompressible flow in Eulerian coordinates. The fluid problem is coupled with the solid problem by enforcing the equality of fluid and solid velocities on the overlap between the fluid and the solid using a penalty method. We then map the fluid velocities and pressure back to the solid, convert them into a fluid stress tensor and subtract the resulting momentum from the solid balance. After solving the solid problem, formulated with linear elasticity and Lagrangian coordinates, we check for convergence of the complete iteration and either stop or continue it by mapping the solid velocities back to the fluid (see Fig. 1 or Alg. 1). The mappings between the domains are realized as volumetric variational transfer operators explained in subsection 2.4. Note, that while we employ the same principles as [16], [30] or [29], our adaption differs from these methods, as we do not use Lagrange multipliers, but a penalty method to enforce the equality of solid and fluid velocities.

We denote the fluid domain with $\Omega_{f}\subset\mathbb{R}^{3}$, the moving solid domain depending on time $t$ with $\Omega^{t}_{s}\subset\Omega_{f}$ , the boundary of the fluid domain with $\Gamma_{f}$ and the boundary of the solid domain with $\Gamma^{t}_{s}$. Both the fluid and the solid boundary are decomposed into a Neumann boundary with superscript N and a Dirichlet boundary with superscript D, such that we have $\Gamma_{f}=\Gamma_{f}^{D}\dot{\cup}\Gamma_{f}^{N}$ and $\Gamma^{t}_{s}=(\Gamma^{t}_{s})^{D}\dot{\cup}(\Gamma^{t}_{s})^{N}$. The outer normal vectors on the boundaries are denoted $n_{f}$ and $n_{s}$. For function spaces we use $L^{2}(\Omega)$ for the space of Lebesgue measurable and square integrable functions and $H^{1}(\Omega)$ for the standard space of weakly differentiable functions in $L^{2}(\Omega)$, where $\Omega\in\{\Omega_{f},\Omega_{s}^{t}\}$.

2.1 Strong formulation

In this section we introduce the strong formulation of the fluid problem and the dynamic solid problem. The fluid problem is formulated using the Navier Stokes equations for incompressible flow and the solid problem with linear elasticity.

We denote with $\sigma_{f}$ the fluid stress tensor, with $\rho_{f}$ the fluid density, with $v_{f}$ the fluid velocity and with $\mu_{f}$ the dynamic fluid viscosity. For the solid problem we introduce the stress tensor $\sigma_{s}$ of the solid, the solid displacement $u_{s}$ and the solid density $\rho_{s}$. For time derivatives, we use the dot notation, i.e. $\dot{u}_{s}:=\frac{\partial}{\partial t}u_{s}$. The strong formulation of the FSI problem then reads as:

[TABLE]

Equations 1 and 2 are the Navier Stokes equations with Dirichlet and Neumann (Eq. 4) boundary conditions (Eq.3). Equation 5 is the equation of motion for the solid, which consists of a linear elastic material obeying Hooke’s law. Again we have Dirichlet (Eq. 6) and Neumann (Eq. 7) boundary conditions. Equation 8 couples the fluid and solid velocities on the overlap $\Omega^{t}_{s}\cap\Omega_{f}$. The relationship between the Cauchy stress $\sigma_{f}$ for a Newtonian fluid and fluid velocity and pressure is given by:

[TABLE]

In the following, we will interchangeably use $v_{s}:=\dot{u}_{s}=\frac{\partial}{\partial t}u_{s}$ for the solid velocity. For the sake of preserving the generality and extendability of our approach, we use time dependent formulations here. However, it is worth mentioning, that the problem could also be formulated in steady state.

2.2 Weak Formulation

We describe the weak formulation for two reasons in particular: First, in finite elements it is ultimately an approximation of the weak formulation that is implemented and secondly, it is by the weak formulation that the second part of the coupling described later in subsection 2.5 can be understood.

We begin by introducing the function spaces for the velocities $v$ and pressure $p$ for the fluid problem:

[TABLE]

For test functions $\delta v\in\mathcal{W}^{v}$ and $\delta p\in\mathcal{W}^{p}$ we can then write the fluid problem in the following weak form:

[TABLE]

For the solid problem we use the following space for the displacements:

[TABLE]

and write it in the weak form:

[TABLE]

2.3 Discretization in space and time

For both the fluid and the solid problem we use a Galerkin finite element method. We mesh $\Omega_{f}$ with a structured rectangular mesh of hexahedral elements and denote the corresponding space of bilinear Lagrange elements for the fluid velocities by $W^{v}$ and the space of bilinear Lagrange elements for the fluid pressure with $W^{p}$. Furthermore, we employed Streamline Upwind Petrov Galerkin (SUPG) and the Pressure Stabilized Petrov Galerkin (PSPG) formulations for stabilization. The solid domain $\Omega_{s}^{t}$ has an unstructured mesh with tetrahedral elements and the corresponding spaces of linear Lagrange elements for the solid pressure, displacements and velocities are denoted $V^{p}$, $V^{u}$ and $V^{v}$. The degrees of freedom for each space are denoted with $n^{\bullet},\,\bullet\in\{W^{v},W^{p},V^{u},V^{p},V^{v}\}$. Both the solid and the fluid problem are time dependent. For the fluid problem we used backward differentiation formula, BDF2 and for the solid problem we used a Newmark scheme with parameters $\beta=\frac{1}{2}$ and $\gamma=\frac{1}{4}$.

2.4 Variational Transfer Operators

One of the main characteristic of our method is that we solve the fluid and the solid problem on different domains with non-matching meshes. The two problems need to be coupled however, i.e. we have to transfer velocities and pressures between the two domains. This is achieved by using variational transfer operators, which provide a stable and flexible way to map variables between arbitrary meshes.

To define the operators, we assume that the solid domain $\Omega_{s}^{t}$ is immersed into the fluid domain $\Omega_{f}$. We further assume that we have a function space $V$ on $\Omega_{s}^{t}$ with elements $v$ and a function space $W$ on $\Omega_{f}$ with elements w, such that $V$ is a subset of $W$, i.e. $V\subset W$ and that the $L^{2}$-scalar product $(v,w)_{L^{2}}=\int_{\Omega_{f}}v\,w\,d\omega$ is well defined. The variational transfer operator is then defined as a mapping $\Pi:V\rightarrow W$, such that for every $v\in V$ the following $L^{2}$-orthogonality holds for all elements $\mu$ of a yet to be defined multiplier space $M$:

[TABLE]

Naturally, in applications we need a discrete matrix representation $T$ of $\Pi$. $T$ is constructed by reformulating Eq.15 and introducing the discrete approximations of the spaces $V,W$ and $M$. These discrete spaces are $V_{h},W_{h}$ and $M_{h}$ with bases $(\phi_{i}^{V_{h}})_{i=1,...,n^{V_{h}}}$, $(\phi_{i}^{W_{h}})_{i=1,...,n^{W_{h}}}$ and $(\phi_{i}^{M_{h}})_{i=1,...,n^{M_{h}}}$, where $n^{V_{h}}$, and $n^{W_{h}}=n^{M_{h}}$ are the dimensions of each space. Now, let $v_{h}$ be an element of $V_{h}$, then with $w_{h}:=\Pi(v_{h})$, we can reformulate Eq. 15 as:

[TABLE]

and using the basis representations of $v_{h},w_{h}$ and $\mu_{h}$ we get:

[TABLE]

With

[TABLE]

we can now define the discrete operator as $T:=D^{-1}B$ and obtain the relationships:

[TABLE]

Note that $D$ and $B$ are essentially volumetric mass matrices, whereby $D$ is a quadratic and $B$ a rectangular matrix . We have some liberty, which multiplier space $M_{h}$ we select. With the most straight forward choice $M_{h}=W_{h}$ we would end up computing the inverse of $D$, which is an unwanted computational expense. Instead we choose a space $M_{h}$ spawned by biorthogonal basis functions [51], which leads to a diagonal matrix $D$ that is easy to invert. More specifically, in this use, operators of this type are called pseudo-$L^{2}$-projections. They have a solid mathematical base and $H^{1}$-stability and $L^{2}$-approximation property are shown in [9]. Following this principle, we assemble the projections:

[TABLE]

In practice, the assembly of the discrete variational operators is highly complex. To compute the individual entries $d_{ik}$ and $b_{ik}$ in Eq. 20 & 21, we need to know all intersections between the elements of the fluid and the solid mesh, in order to determine the common support of the bases $\phi^{W_{h}},\phi^{V_{h}}$ and $\phi^{M_{h}}$. In parallel computing environments these elements might lie on different compute nodes, so additional communication between the nodes is necessary. It is only with the recent development of the MOONoLith library [23], that the assembly routines have become efficient enough to make this immersed boundary approach practicable.

2.5 Coupling

The fluid and the solid problem are fully coupled. First, we enforce the equality of fluid and solid velocities on the overlap of the immersed solid and the fluid using a penalty method. Secondly, we transfer the fluid velocity and pressure to the solid problem and use them to compute a fluid force term. Note that the coupling relies on the transfer operators which are assembled by computing volumetric integrals. Hence no explicit representation of the boundary between the fluid and the solid is ever required.

To couple the fluid with the solid we apply a penalty method, which enforces the equality of fluid and solid velocities. Let $\epsilon$ be then penalty parameter and $\Pi_{\dot{u}_{s}}$the non-discrete form of $\mathbf{T}_{\dot{u}_{s}}$. Then we add the following term to the weak formulation of the fluid problem (Eq. 12):

[TABLE]

In its discrete form, the overlap between the solid and the fluid is determined by constructing the transfer operator $\mathbf{T}_{I}:V^{i}\rightarrow W^{i}$, which maps an indicator variable $i$ from the solid to the fluid mesh, whereby $V^{i}$ and $W^{i}$ are the spaces of (bi-)linear Lagrange elements on each domain. Let $\mathcal{T}^{W}$ be the mesh of the space $W^{i}$. The overlap $\Omega_{f}\cap\Omega_{s}$ is defined as all nodes $p$ in $\mathcal{T}^{W}$ with $T_{I}(i)\geq\alpha$ ($\alpha$ being a mesh dependent threshold value).

To couple the solid with the fluid, we transfer the fluid velocity $v_{f}$ and pressure $p_{f}$ to the solid simulation after every iteration. We then recompute the Cauchy stress according to Eq. 9 and subtract the term for the contributions of the fluid forces, i.e. the second term $\int_{\Omega_{f}}\sigma_{f}:\nabla(\delta v)\,dV$ in Eq. 12, from the weak formulation of the solid (Eq. 14). In addition, we should also subtract the dynamic contributions of the fluid phase from the solid phase. But since in geophysical applications the Young’s Modulus of the solid is normally on the order of GPa and solid velocities are consequentially small, we omit the dynamic contributions of the fluid on the solid.

Note that using a penalty method can be seen as a simplification of the coupling described in [16] and [29]. There, Lagrange multipliers are used, which can be interpreted as forces acting on the fluid. From a theoretical point of view, using Lagrange multipliers is superior to penalty methods as the forces can be mapped back from the fluid to the solid using the transposed of the mapping operators, resulting in a variationally consistent scheme. There are, up to the knowledge of the authors, no equivalent results for using a penalty method. On the other hand, the penalty method has the advantage that no saddle point problem needs to be solved, instead we can use a standard iterative solver for the regularized problem. Furthermore, we will show in the section 3 with numerical experiments, that using the penalty method does not impact the capability of resolving the solid-fluid boundary.

2.6 Algorithm

The coupled problem is solved in a fixed point iteration in up to $L_{\mathrm{FP}}$ steps. We formulate one step $l$ by introducing the coupling terms to the weak formulations of Eq. 12 - 14. For the solid problem we define the operator $S(u):=\mathcal{V}^{u}\rightarrow(\mathcal{V}^{u})^{*}$ as:

[TABLE]

and the right hand side $R_{s}(\delta u)$ as:

[TABLE]

Note, that we have added the contributions from the fluid $\int_{\Omega_{f}\cap\Omega_{s}}\sigma_{f}:\nabla(\delta v)\,dV$ to $R_{s}$. This will become relevant later in the context of a fixed point iteration, as we will use the results of a previous iteration $l-1$ to form $\sigma_{f}$. For the fluid problem we define $F(v,p):=\mathcal{W}^{v}\times\mathcal{W}^{p}\rightarrow(\mathcal{W}^{v}\times\mathcal{W}^{p})^{*}$ as:

[TABLE]

and the right hand side $R_{f}(\delta v)$ as:

[TABLE]

Using the discretizations from Section 2.3 we obtain the discrete representatives $\mathbf{S},\mathbf{F},\mathbf{R_{s}}$ and $\mathbf{R}_{f}$ from $S$(Eq.28 ), $F$(Eq.30), $R_{s}$(31), and $R_{f}$(Eq.31), respectively. We then discretize the penalty term from Eq. 27 and split it into two parts $\mathbf{P}_{s}$ and $\mathbf{P}_{f}$. Then we denote with $\mathbf{y}_{s}$ and $\mathbf{y}_{f}$ the solutions of the solid and fluid problem, whereby $\mathbf{y_{f}}$ contains fluid velocities and pressure, and write the coupled problem as:

[TABLE]

where $\mathbf{P}_{s}$is defined as:

[TABLE]

and $\mathbf{P}_{f}$ as:

[TABLE]

Here, $\mathbf{I}_{\Omega_{f}\cap\Omega_{s}}$ is a diagonal matrix containing zeros and ones according to the values of $T_{I}(i)$, i.e.:

[TABLE]

where $l,k=1,...,n^{W^{v}}$ and $\mathbf{M}_{f}$ is the mass matrix obtained from the discretization on $W^{v}$. The matrices $\mathbf{0}$ contain all the entries associated with the degrees of freedom from the fluid pressure. By adding the superscripts $l,l-1$ for the iteration steps of the fixed-point iteration, we can formulate the staggered approach to solve the coupled problem by algorithm 1.

We note that the fixed-point iteration is carried out in each time-step. For ease of presentation, here we have dropped any indices connected to the time discretization. We furthermore note that the transfer of the solid velocities is realized by means of the multiplication of $\mathbf{y}_{s}$ with $\mathbf{P}_{s}$, as the latter contains the discrete projection $\mathbf{T}_{v_{s}}$. The fixed-point iteration shown in Alg. 1 can be interpreted as a block-Gauss-Seidel iteration. We will see later in sec. 3.1 that in practice convergence of the overall FSI problem is only mildly dependent on $\epsilon$.

2.7 Implementation

The variational transfer operators are implemented as userobjects in the MOOSE framework [13], employing MOONoLith [23], Utopia [57], and libMesh [21]. The MOOSE framework is also used to form the weak formulations. All components are designed for parallel computing.

3 Numerical examples

The following section validates the presented numerical modeling approach with 2D parallel plate (Section 3.1.1) and 3D tunnel flow (Section 3.1.2) benchmark problems, before applying the approach to real fracture geometries (Section 3.2).

Throughout the simulations, we use the same material parameters, which are representative of a hard granodiorite rock, saturated with water. For the rock we use a Poisson ratio of 0.33, a Young’s modulus of 1$\times 10^{10}$ Pa, and a density of 2.75 $g/cm^{3}$. For water at about 20∘C we use a kinematic viscosity of 1 ${mm}^{2}/s$ and a density of 1 $g/cm^{3}$. The penalty parameter $\epsilon$ was set to $10^{-6}$, the tolerance $tol_{\mathrm{FP}}$ to $10^{-8}$ and $L_{\mathrm{FP}}$ to 10. Simulations were conducted on the cluster of the Institute of Computational Science in Lugano, Switzerland, using up to 10 compute nodes (2 x Intel Xeon E5-2650 v3 @ 2.30GHz) with 16 CPU’s each.

3.1 Benchmark problems

We designed a two benchmark simulations to assess the method’s capability of detecting boundaries and simulating realistic fluid flow. The first benchmark is a 2D parallel plates example and a second one a simulation of tunnel or Poiseuille flow in 3D. The results of the simulations are compared to analytical solutions and/or the solution of an equivalent setup for fluid flow with rigid boundaries.

Apart from the dimension and the used geometries, the equivalent Navier-Stokes simulation experiments have the same setup. We apply zero-Dirichlet boundary conditions on the entire solid boundary, thereby creating a rigid object. On the fluid domain we apply a pressure gradient of 0.1 MPa from the left inlet to the right outlet, forcing the fluid to flow (Fig. 5).

The benchmark problems were also used to investigate the influence of the penalty parameter. We refer to Fig. 2, where we ran the 3D simulation for tunnel flow with the different penalty parameters $\epsilon=10^{-1},10^{-3},10^{-6}$, and compared the resulting discrete solutions to the analytical solution. Eventually, the penalty parameter $\epsilon$ was set to $10^{-6}$. In addition to solving the FSI problem with different penalty parameters, we also ran the 3D tunnel flow for different mesh widths $h$ (with $h=h_{f}\approx h_{s}$) to assess how the convergence of the overall and the fluid and the solid problem depended on $h$ and $\epsilon$. The results collected in Fig. 3 show only a mild dependency on $\epsilon$ in terms of the number of necessary fixed-point iterations. For the inner linear problems we used Krylov-space methods, more precisely a biconjugate gradient method (BiCG) for the solid problem and GMRES for the fluid problem. We note that by using standard preconditioners for the solid problem, such as, e.g., algebraic multigrid methods, the increase in the iteration numbers could be avoided.

3.1.1 2D Parallel plate

In this experiment, the fluid domain has a length of 3 mm and a width of 1 mm. The immersed solid domain consists of two parallel plates with a thickness of 0.1 mm each, located at the boundary of the fluid domain, forming a channel (Fig. 4a). The solid mesh consists of 1’270 nodes with mesh width $h_{s}$ of 0.025 mm. The fluid mesh has 1’281 nodes and a mesh width $h_{f}$ of 0.08 mm. As a comparison and proxy for an analytical result, we set up an equivalent Navier-Stokes finite element simulation. Its length is 3 mm and in order to mimic the same solid boundary as in our benchmark experiment, the width is 0.8 mm. We apply the same pressure gradient and set the fluid velocity to zero across the lateral boundaries (Fig. 4b). The resulting fluid velocities in the center of the channel (i.e. at x=1.5 mm) are compared between the immersed boundary approach and the classic setup and show good agreement (Fig. 7). See also Fig. 6 for global comparison.

3.1.2 3D Tunnel-flow

To test the immersed boundary approach in three dimensions, we designed a benchmark problem of 3D tunnel flow. As this is a classic set-up of Poiseuille flow, we can compare the results of the immersed boundary approach with the analytical solution.

The fluid domain has a length of 3 mm and a width of 1 mm (Fig. 8b). The structured mesh has 203’401 nodes with a mesh width $h_{f}$ of 0.05 mm. The solid domain has the same length and width as the fluid domain (Fig. 8a), but along the x-axis, we have cut out a cylinder with a diameter of 0.4 mm. The unstructured mesh has 63’866 nodes and a mesh width $h_{s}$ of around 0.02 mm.

To compare the result with the analytical solution, the velocity profile in the x-plane in 2D is extracted and plotted together with an analytical solution derived from Poiseuille’s law in 3D. Fig. 9a shows the alignment of the velocity profile of the immersed boundary approach with the analytical solution. This observation is confirmed, when we compare the 1D velocity profile at y = 0.5 mm with the same analytical solution (Fig. 9b).

In conclusion of this benchmark section, we observe that the benchmark problems demonstrate the ability of the immersed boundary approach to resolve the boundary between the solid and fluid and to simulate the fluid flow accordingly in two and three dimensions.

3.2 Rough Fracture Surfaces

This section investigates fluid flow through rough rock fracture geometries and the respective responses to increasing normal load or fluid injection pressures. For this purpose, two numerical experiments are conducted: (1) Fracture closure with increasing confining stresses normal to the fracture and (2) fracture opening with increasing fluid pressures in the fracture.

We obtain the geometries for the solid domain $\Omega_{s}^{t}$ from the contact simulations in [34]. There, numerical experiments are conducted on modeled rock fracture geometries under increasing normal loads from 0 to 22 MPa, using linear elasticity and linearized contact conditions (see e.g. [20]). The contact conditions between the nonmatching surface meshes were thereby resolved using a mortar approach with dual Lagrange multipliers [52, 35, 34]. The geometries originate from a granodiorite rock specimen from the Grimsel Test Site located in the Swiss Alps, which were previously used in laboratory experiments [45, 48, 46]. The specimen is of cylindrical shape, with a fracture normal to the cylinder axis in the center of the specimen. Fracture topographies were obtained through photogrammetry scans, which allowed the creation of numerical meshes for the upper and lower half of the specimen in [34]. Because we want to focus on the processes near the fracture, we do not simulate the entire cylinder. Instead, we extract a rectangular cubical region around the fracture (red part in Fig. 10), in order to immerse the resulting solid in a rectangular structured mesh for the fluid. This helps setting up the fluid problem, as the role of the boundaries and the boundary conditions to be used become clearer. In addition, the resolution in the fracture region is increased, as nodes to model peripheral regions are not needed, as would be the case if the entire cylinder would be used.

The original cylinder from the contact simulations by [34], had a diameter of 122 mm and a height of 240 mm. From this we extract solid geometries with lengths and widths of 81 mm. Depending on the experiment type, the height is either 9 mm or 23 mm.

3.2.1 Fluid flow in closing fracture

These simulations model the change of flow patterns in a fracture which is subjected to increasing normal loads. Here, we use thinner geometries with a height of 9 mm, which are obtained from different stages of the closure process (i.e., confining pressures ranging from 0.25 to 20 MPa), that is, the geometry itself reflects the fracture under different confining pressures, but does not have a normal loading in the HM-simulation, in order to reduce the complexity of the setup.

Solid meshes resolution ranges between 100’637 and 123’231 nodes. The meshes are refined along the fracture surface to achieve a higher resolution at the fluid-solid boundary (Fig. 11a). The structured fluid mesh has the same length, width and height as the solid mesh. It has 1’212’741 nodes and is more refined in the z-direction to gain a better resolution parallel to the fracture plane (Fig. 11b).

We apply Dirichlet boundary conditions to the top and bottom boundaries of the solid (Fig. 12a). On the fluid, we apply non-slip conditions at the top, bottom, front and back boundaries. We set the fluid pressure on the left (1 MPa) and right side (0 MPa) of the fluid domain (Fig. 12b), thereby establishing a fluid pressure gradient across the fracture.

We visualize the setup in Fig. 13, where we increase the aperture of the fracture in order to clarify the concept. The picture shows the velocity streamlines of the fluid, originating from a sphere in the center of the open fracture geometry. Fluid flow is tortuous against the rough fracture surface with no intersections between the fluid and the solid.

Fig. 14 shows the fluid flow rate averaged in the z-direction (to ease visual comparison) and the aperture distribution at a confining pressure of 0.25 MPa. The aperture field is strongly heterogeneous (Fig. 14a), resulting in the characteristic preferential fluid flow paths [42, 6, 43, 28] that can be seen in Fig. 14b. The two plots are combined in Fig. 14c to illustrate the dependence of the fluid flow field on the aperture field more clearly. As expected, flow channels form in regions with large apertures (white) or around large regions of fracture wall contact (black).

Fig. 15 depicts the flow path evolution under increasing normal loads. We again average flow rates in the z-direction and view aperture values from above to ease visual comparison. Here, we show the spatial distribution of fluid flow rates, as in Fig. 14b, under normal loads of 0.25, 8 and 10 MPa. For small loads, the majority of the fracture is essentially open (15b), resulting in fluid flow across large areas of the fracture (15a). Nonetheless, areas with large aperture widths experience particular high fluid flow rates, while regions, where the fracture walls are already in contact, are void of fluid flow. With increasing normal load, fracture areas, which were open for small normal loads, are now closed and do not experience significant fluid flow rates. The importance of heterogeneous aperture width fields and their fine resolution is apparent across all loading stages, as large- and small-aperture regions always govern fluid flow fields.

To illustrate the rapid changes in fracture permeability that a fracture under load experiences, Fig. 16 shows the non-linear behavior of the fluid inflow rate into the fracture versus normal load across the fracture [50, 37, 36, 49, 44, 46]. The inflow rate is calculated by integrating the fluid velocity over the x-plane cross-section at the left side (inflow direction) of the y-plane. The total fluid flow rate through the fracture changes most drastically at initial normal load increases during small normal loads, where fractures are known to experience the largest aperture closure relative to increases in normal load [3, 25, 54, 46, 22]. While fracture closure is more pronounced during early normal loading stages, fluid flow paths through regions of large aperture widths persist even at higher loading stages, as their large-aperture regions may not close entirely (Fig. 15).

3.2.2 Fracture opening

These simulations aim to demonstrate fracture opening with the immersed boundary numerical modeling approach. In particular, we test fracture opening under increasing fluid pressure in the fracture. For this purpose, we compare results for constant fluid pressures in the fracture of 10 and 20 MPa respectively.

A fracture under a normal load of 20 MPa serves as our starting geometry, i.e. in a quasi closed state (see last row in Fig. 15). Compared to the previous section, the solid mesh is thicker, as the focus in this section is on the mechanical response of the fracture due to changes in fluid pressure. The overall height of the solid is $\sim$23 mm, resulting in more finite element nodes along the z-axis, which helps in resolving the mechanical displacements.

We apply Dirichlet boundary conditions at the top and bottom boundaries of the solid domain and we set a uniform fluid pressure, that is the same fracture inlet pressure and outlet pressure, at the left and right boundaries of the fluid domain (Fig. 17). In addition, no-slip boundary conditions are applied at the top, bottom, front and back boundaries of the fluid domain.

Figure 18 shows the response of the fracture to a fluid pressure of 10 MPa (left column) and 20 MPa (right column). The top row of the figure is a cross-section through the fracture. Blue colors denote displacements downward and red colors upward. Due to the Dirichlet boundary conditions, applied to the top and bottom boundaries of the solid domain (Fig. 17), displacements converge to zero as these solid boundaries are reached. The middle and bottom row of Figure 18 show displacements in the z-direction for the upper and for the lower fracture surfaces or planes, respectively. This demonstrates, that the immersed boundary approach can simulate the fluid pressure acting on the solid, thereby causing a mechanical response in the solid.

4 Conclusions

We presented fully coupled hydro-mechanical (HM) numerical simulations of fluid-saturated, rough-walled fractures. To this end, we adapted an immersed boundary method. The approach allows us to solve the coupled problem without an explicit representation of the complex boundary between the solid and the fluid. Instead, the coupling between the two domains relies entirely on the computation of volume integrals. This immensely simplifies the setup of HM-simulations. Additionally, our implementation solely relies on software components that are designed for parallel computing.

We validated our approach with two and three dimensional benchmarks and also with geometries from a real fracture under load, which was obtained from a laboratory rock core specimen. The numerical experiments show, that the method replicates Poiseuille flow in three dimensions, that it replicates the change of fluid flow patterns in a fracture under increasing normal loads (confining pressures) across the fracture and that it can also simulate the deformation of the solid when the fluid pressure inside the fracture changes. The latter was demonstrated by increasing the fluid pressure inside the fracture, causing fracture opening.

The methodology enables the computationally efficient investigation of complex fluid flow patterns in rough-walled fractures that exhibit complex fracture surface geometries. The method can also numerically simulate the interaction between the fluid and the solid. In particular, the ability to simulate the interaction between the solid and the fluid, due to increases in the fluid pressure inside the fracture, can make this method a valuable tool for investigating HM processes, for example, during hydraulic stimulation.

We employ variational transfer operators to map fluid pressure and velocity between the fluid and the solid domains. The same concept could also be used to extend HM simulations, using staggered approaches and coupling, to include thermal processes.

Acknowledgements: We gratefully acknowledge funding by the Swiss Competence Center for Energy Research - Supply of Electricity (SCCER-SoE), by Innosuisse - Swiss Innovation Agency under Grant Number 28305.1 and the Swiss Federal Office of Energy (SFOE) under Grant Number SI/500676-02. The Werner Siemens-Stiftung (Foundation) is thanked for its support of the Geothermal Energy and Geofluids group.