# A two-scale approach for fluid flow in fratured porous media

Julien Réthoré, René de Borst, Marie-Angèle Abellan

► To cite this version:

Julien Réthoré, René de Borst, Marie-Angèle Abellan. A two-scale approach for fluid flow in fratured porous media. International Journal for Numerical Methods in Engineering, 2007, pp.780-800. 10.1002/nme.1962 . hal-00380792

HAL Id: hal-00380792

https://hal.science/hal-00380792v1

Submitted on 4 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A two-scale approach for fluid flow in fractured porous media

Julien Réthoré¹, René de Borst¹, ², *, † and Marie-Angèle Abellan³

¹Faculty of Aerospace Engineering, Delft University of Technology, Delft, Netherlands

²LaMCoS-UMR CNRS 5514, INSA de Lyon, Villeurbanne, France

³LTDS-ENISE-UMR CNRS 5513, Saint-Etienne, France

A two-scale numerical model is developed for fluid flow in fractured, deforming porous media. At the microscale the flow in the cavity of a fracture is modelled as a viscous fluid. From the micromechanics of the flow in the cavity, coupling equations are derived for the momentum and the mass couplings to the equations for a fluid-saturated porous medium, which are assumed to hold on the macroscopic scale. The finite element equations are derived for this two-scale approach and integrated over time. By exploiting the partition-of-unity property of the finite element shape functions, the position and direction of the fractures is independent from the underlying discretization. The resulting discrete equations are non-linear due to the non-linearity of the coupling terms. A consistent linearization is given for use within a Newton-Raphson iterative procedure. Finally, examples are given to show the versatility and the efficiency of the approach, and show that faults in a deforming porous medium can have a significant effect on the local as well as on the overall flow and deformation patterns.

KEY WORDS: fracture; porous medium; multiscale method; multiphase medium; fluid flow

## 1. INTRODUCTION

Flow of fluids in deforming porous media has been a topic of attention in engineering science ever since the seminal works of Terzaghi [1] and Biot [2], see [3] for a recent account. Indeed, the subject is crucial for understanding and predicting the physical behaviour of many systems of interest. We mention geotechnical and petroleum engineering as more classical examples, but the recognition that many human tissues are, in fact, porous media and must be analysed as such, has put flow in porous media in a totally new perspective. Because of the complicated structure and

<footnote>*Correspondence to: René de Borst, Faculty of Aerospace Engineering, Delft University of Technology, P.O. Box 5058, NL-2600 GB Delft, Netherlands.</footnote>

<footnote>†E-mail: r.deborst@tudelft.nl</footnote>

<footnote>Contract/grant sponsor: European Commission; contract/grant number: HPRN-CT-2002-00198</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. Flow characteristics in a fractured, deforming porous medium.</figcaption></figure>

functioning of human tissues, the classical two-phase theory has been extended in recent times to three- and four-phase media, taking into account ion transport and electrical charges as well [4–6]. A general framework for accommodating multi-field problems has been presented by Jouanna and Abellan [7].

In spite of the importance of the subject, flow in damaged porous media has received little attention. Yet, the presence of damage, such as cracks, faults, shear bands, can markedly change the physical behaviour. Indeed, as illustrated in Figure 1, the physics of the flow within such discontinuities will be very different from the flow of the interstitial fluid in the deforming bulk material. These differences in physics will therefore affect the flow pattern and the deformations in the vicinity of the discontinuity. As we will show at the end of the paper, the local differences in flow characteristics even influence the flow and deformations in the entire body of interest, as the porosity on this scale is altered.

In this contribution, we will develop a general numerical model for flow in fractured porous media. Building on earlier work in which we have constructed a numerical model for fracture propagation in deforming, fluid-saturated porous media [8], we now extend the theory to include flow inside such cavities directed tangentially to the discontinuity. This is achieved by a priori adopting a two-scale approach. At the fine scale the flow in the cavity created by the discontinuity is modelled as a viscous fluid using Stokes equations. Since the cross-sectional dimensions of the cavity are small compared to its length, the flow equations can be averaged over the cross-section of the cavity. The resulting equations provide the momentum and mass couplings to the standard equations for a porous medium, which are assumed to hold on the macroscopic scale.

Numerically, the two-scale model which ensues imposes some requirements on the interpolation of the displacement and pressure fields near the discontinuity. The displacement field must be discontinuous across the cavity. Furthermore, the micromechanics of the flow within the cavity require that the flow normal to the cavity is discontinuous, and in conformity with Darcy’s relation
which is assumed to hold for the porous medium, the normal derivative of the fluid pressure field must also be discontinuous from one face of the cavity to the other. For arbitrary discretizations, these requirements can still be satisfied, namely by exploiting the partition-of-unity property of finite element shape functions [9], as has been done successfully in applications to cracking in single-phase media [10–20].

To provide a proper setting, we will first briefly recapitulate the governing equations for a deforming porous medium under quasi-static loading conditions. The strong as well as the weak formulations will be considered, since the latter is crucial for incorporating the micromechanical flow model properly. This micromechanical flow model is the subject of the subsequent section, whereafter it will be demonstrated how the momentum and mass coupling of the micromechanical flow model to the macroscopic porous medium can be accomplished *via* the weak formulation. Time integration and a consistent linearization of the resulting equations, which are non-linear due to the coupling terms, complete the numerical model. Finally, example calculations are given of a body with existing faults, which can have an arbitrary direction and an arbitrary location within the body. Apart from demonstrating the effectiveness of the two-scale approach, the calculations show that the influence of the presence of discontinuities on flow and deformation patterns can be significant.

## 2. GOVERNING EQUATIONS FOR THE BULK

### 2.1. *Strong form*

The bulk is considered as a two-phase medium subject to the restrictions of small displacement gradients and small variations in the concentrations [7]. Furthermore, it is assumed that there is no mass transfer between the constituents, that convective terms and the gravity acceleration can be neglected, and that the processes which we consider, occur isothermally and quasi-statically. With these assumptions, the balances of linear momentum for the solid and the fluid phases read, e.g. [21]:

$$ \nabla \cdot \boldsymbol{\sigma}_{\pi} + \hat{\boldsymbol{p}}_{\pi} = \mathbf{0} \qquad (1) $$

with $\boldsymbol{\sigma}_{\pi}$ the stress tensor of constituent $\pi$. In the remainder of this paper we follow $\pi = s, f$, with $s$ and $f$ denoting the solid and fluid phases, respectively. Further, $\hat{\boldsymbol{p}}_{\pi}$ is the source of momentum for constituent $\pi$ from the other constituent, which takes into account the possible local drag interaction between the solid and the fluid. Evidently, the latter source terms must satisfy the momentum production constraint:

$$ \sum_{\pi=s,f} \hat{\boldsymbol{p}}_{\pi} = \mathbf{0} \qquad (2) $$

Adding both momentum balances, and taking into account Equation (2), one obtains the ‘standard’ equilibrium equation for the mixture

$$ \nabla \cdot \boldsymbol{\sigma} = \mathbf{0} \qquad (3) $$

where

$$ \boldsymbol{\sigma} = \boldsymbol{\sigma}_s + \boldsymbol{\sigma}_f \qquad (4) $$
Under the same assumptions as for the balance of momentum, one can write the balance of mass for each phase as

$$ \frac{\partial \rho_{\pi}}{\partial t} + \rho_{\pi} \nabla \cdot \mathbf{v}_{\pi} = 0 \qquad (5) $$

with $\rho_\pi$ is the apparent mass density and $\mathbf{v}_\pi$ the absolute velocity of constituent $\pi$. We multiply the mass balance for each constituent $\pi$ by its volumic ratio $n_\pi$, add them and utilize the constraint

$$ \sum_{\pi=s,f} n_{\pi} = 1 \qquad (6) $$

to give

$$ \nabla \cdot \mathbf{v}_s + n_f \nabla \cdot (\mathbf{v}_f - \mathbf{v}_s) + \frac{n_s}{\rho_s} \frac{\partial \rho_s}{\partial t} + \frac{n_f}{\rho_f} \frac{\partial \rho_f}{\partial t} = 0 \qquad (7) $$

The change in the mass density of the solid material is related to its volume change by

$$ \nabla \cdot \mathbf{v}_s = - \frac{K_s n_s}{K_t \rho_s} \frac{\partial \rho_s}{\partial t} \qquad (8) $$

with $K_s$ the bulk modulus of the solid material and $K_t$ the overall bulk modulus of the porous medium. Using the definition of the Biot coefficient, $\alpha = 1 - K_t/K_s$ [3], this equation can be rewritten as

$$ (\alpha - 1) \nabla \cdot \mathbf{v}_s = \frac{n_s}{\rho_s} \frac{\partial \rho_s}{\partial t} \qquad (9) $$

For the fluid phase, a phenomenological relation is assumed between the incremental changes of the apparent fluid mass density and of the fluid pressure $p$ [3]:

$$ \frac{1}{Q} d p = \frac{n_f}{\rho_f} d \rho_f \qquad (10) $$

with the compressibility modulus

$$ \frac{1}{Q} = \frac{\alpha - n_f}{K_s} + \frac{n_f}{K_f} \qquad (11) $$

where $K_f$ is the bulk modulus of the fluid. Inserting relations (9) and (10) into the balance of mass of the total medium, Equation (7), gives:

$$ \alpha \nabla \cdot \mathbf{v}_s + n_f \nabla \cdot (\mathbf{v}_f - \mathbf{v}_s) + \frac{1}{Q} \frac{\partial p}{\partial t} = 0 \qquad (12) $$

The governing equations, i.e. the balance of momentum of the saturated medium, Equation (3), and the balance of mass, Equation (12), are complemented by the kinematic relation,

$$ \varepsilon_s = \nabla^s \mathbf{u}_s \qquad (13) $$

with $\mathbf{u}_s$, $\varepsilon_s$ the displacement and strain fields of the solid, respectively, the superscript s denoting the symmetric part of the gradient operator, and an incrementally linear stress-strain relation for
<figure><img src="image_3.png" /><figcaption>Figure 2. Boundary conditions and internal discontinuity $\Gamma_d$ for body $\Omega$.</figcaption></figure>

the solid skeleton. The effective stress increment in the solid skeleton, $d\boldsymbol{\sigma}'_s$ is related to the strain increment $d\boldsymbol{\varepsilon}_s$ by an incrementally linear stress–strain relation for the solid skeleton,

$$d\boldsymbol{\sigma}'_s = \bar{\mathbf{D}}^{\text{tan}} : d\boldsymbol{\varepsilon}_s \qquad (14)$$

where $\bar{\mathbf{D}}^{\text{tan}}$ is the fourth-order tangent stiffness tensor of the solid material and the d symbol denotes a small increment. Since the effective stress in the solid skeleton is related to the partial stress by $\boldsymbol{\sigma}'_s = \boldsymbol{\sigma}_s/n_s$, the above relation can be replaced by

$$d\boldsymbol{\sigma}_s = \mathbf{D}^{\text{tan}} : d\boldsymbol{\varepsilon}_s \qquad (15)$$

where the notation $\mathbf{D}^{\text{tan}} = n_s \bar{\mathbf{D}}^{\text{tan}}$ has been used. In the examples, a linear-elastic behaviour of the bulk material has been assumed, and we have set $\mathbf{D}^{\text{tan}} = \mathbf{D}$, the linear-elastic stiffness tensor. For the pore fluid flow, Darcy’s relation for isotropic media is assumed to hold,

$$n_f(\mathbf{v}_f - \mathbf{v}_s) = -k_f \nabla p \qquad (16)$$

with $k_f$ the permeability coefficient of the porous medium. The boundary conditions (see Figure 2)

$$\mathbf{n}_{\Gamma} \cdot \boldsymbol{\sigma} = \mathbf{t}_{p}, \quad \mathbf{u} = \mathbf{u}_{p} \qquad (17)$$

hold on complementary parts of the boundary $\partial\Omega_t$ and $\partial\Omega_u$, with $\Gamma = \partial\Omega = \partial\Omega_t \cup \partial\Omega_u$, $\partial\Omega_t \cap \partial\Omega_u = \emptyset$, $\mathbf{t}_p$ being the prescribed external traction and $\mathbf{u}_p$ the prescribed displacement at the boundary, and

$$n_f(\mathbf{v}_f - \mathbf{v}_s) \cdot \mathbf{n}_{\Gamma} = \mathbf{q}_p \cdot \mathbf{n}_{\Gamma}, \quad p = p_p \qquad (18)$$

hold on complementary parts of the boundary $\partial\Omega_q$ and $\partial\Omega_p$, with $\Gamma = \partial\Omega = \partial\Omega_q \cup \partial\Omega_p$ and $\partial\Omega_q \cap \partial\Omega_p = \emptyset$, $\mathbf{q}_p \cdot \mathbf{n}_{\Gamma}$ and $p_p$ being the prescribed outflow of pore fluid and the prescribed pressure, respectively.

## 2.2. Weak form and coupling

To arrive at the weak form of the balance equations, we multiply the momentum balance (3) and the mass balance (12) by kinematically admissible test functions for the displacements $\eta$ of the
skeleton and for the pressure $\zeta$. Substitution into Equations (3) and (12), using Darcy's relation (16), integrating over the domain $\Omega$ and using the divergence theorem leads to the corresponding weak forms

$$ \int_{\Omega} (\nabla \cdot \boldsymbol{\eta}) \cdot \boldsymbol{\sigma} \, d\Omega + \int_{\Gamma_d} [\boldsymbol{\eta} \cdot \boldsymbol{\sigma}] \cdot \mathbf{n}_{\Gamma_d} \, d\Omega = \int_{\Gamma} \boldsymbol{\eta} \cdot \mathbf{t}_p \, d\Omega \quad (19) $$

and

$$ \begin{aligned} & - \int_{\Omega} \alpha \zeta \nabla \cdot \mathbf{v}_s \, d\Omega + \int_{\Omega} k_f \nabla \zeta \cdot \nabla p \, d\Omega - \int_{\Omega} \zeta Q^{-1} \dot{p} \, d\Omega \\ & + \int_{\Gamma_d} \mathbf{n}_{\Gamma_d} \cdot [\zeta n_f (\mathbf{v}_f - \mathbf{v}_s)] \, d\Gamma = \int_{\Gamma} \zeta \mathbf{n}_{\Gamma} \cdot \mathbf{q}_p \, d\Gamma \end{aligned} \quad (20) $$

where the shorter notation of a superimposed dot has been introduced to denote differentiation with respect to time. Because of the presence of a discontinuity inside the domain $\Omega$, the power of the external tractions on $\Gamma_d$ and the normal flux through the faces of the discontinuity are essential features of the weak formulation. Indeed, these terms enable the coupling between the bulk and the inner cavity of the discontinuity.

The mechanical coupling comes from the tractions applied by the fluid in the discontinuity onto the discontinuity faces. We assume that the tractions are the same on each side of the discontinuity and that they are imposed by the fluid pressure inside the discontinuity. Because of the continuity from the cavity to the bulk, this gives

$$ \boldsymbol{\sigma} \cdot \mathbf{n}_{\Gamma_d} = -p \mathbf{n}_{\Gamma_d} \quad \text{on } \Gamma_d \quad (21) $$

Hence, the weak form of the balance of momentum becomes

$$ \int_{\Omega} (\nabla \cdot \boldsymbol{\eta}) \cdot \boldsymbol{\sigma} \, d\Omega - \int_{\Gamma_d} [\boldsymbol{\eta}] \cdot p \mathbf{n}_{\Gamma_d} \, d\Gamma = \int_{\Gamma} \boldsymbol{\eta} \cdot \mathbf{t}_p \, d\Gamma \quad (22) $$

As a consequence, the test function for the pressure, $\zeta$, also has the same value at both faces of the cavity, and the mass transfer coupling term can be rewritten as follows:

$$ \int_{\Gamma_d} \mathbf{n}_{\Gamma_d} \cdot [\zeta n_f (\mathbf{v}_f - \mathbf{v}_s)] \, d\Gamma = \int_{\Gamma_d} \zeta n_f \mathbf{n}_{\Gamma_d} \cdot [\mathbf{v}_f - \mathbf{v}_s] \, d\Gamma = \int_{\Gamma_d} \zeta \mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d \, d\Gamma \quad (23) $$

where $\mathbf{q}_d$ denotes the fluid flux through the faces of the discontinuity.

The above identities for the coupling of the mass transfer can be interpreted as follows. Part of the fluid that enters the cavity through one of its faces flows away tangentially, that is in the cavity. Therefore, the fluid flow normal to the discontinuity is discontinuous. Because the fluid flow between the cavity and the surrounding porous medium has to be continuous at each of the faces of the discontinuity, and because the fluid velocity is related to the pressure gradient via Darcy's law, the gradient of the pressure normal to the discontinuity must be discontinuous across the discontinuity. Next, the influence at the 'macro'-scale of the tangential 'micro'-flow inside the discontinuity—the normal flux discontinuity—will be quantified.
<figure><img src="image_4.png" /><figcaption>Figure 3. Cavity geometry.</figcaption></figure>

### 3. MODEL FOR THE FLOW INSIDE THE CAVITY

Subject to the same assumptions as listed in the preceding section for the bulk, the balance of momentum for the fluid phase in the cavity reads (Figure 3)

$$ \nabla \cdot \boldsymbol{\sigma}_m = \mathbf{0} \qquad (24) $$

The subscript m denotes that a quantity is related to the ‘micro’-flow. Assuming a Newtonian fluid, the balance of momentum inside the cavity gives

$$ \mu \nabla \cdot \nabla \mathbf{v}_m = \nabla p_m \qquad (25) $$

where $\mu$ is the viscosity of the fluid and $p_m$ the pressure in the ‘micro’-flow. We assume that the dimensions of the cross section of the cavity are small in comparison to the length of the discontinuity (for the two-dimensional problem of Figure 3 this dimension is the cavity width $2h$). A quasi-one-dimensional flow is assumed inside the cavity. Hence, the vector equation (25) leads to

$$ 0 = \frac{\partial p_m}{\partial y} \qquad (26) $$

in the normal direction $\mathbf{n}_{\Gamma_d}$ and

$$ \mu \frac{\partial^2 v}{\partial y^2} = \frac{\partial p_m}{\partial x} \qquad (27) $$

in the tangential direction $\mathbf{t}_{\Gamma_d}$. $x$, $y$ give the position of a material point inside the cavity in the $\mathbf{t}_{\Gamma_d}$, $\mathbf{n}_{\Gamma_d}$ co-ordinate system. $v = \mathbf{v}_m \cdot \mathbf{t}_{\Gamma_d}$ is the velocity tangential to the discontinuity. From the first equation, $p_m$ is constant through the cavity and, for continuity reasons, it is equal to $p$

$$ p_m = p \qquad (28) $$
Integrating the second equation twice with respect to the y-co-ordinate, we obtain the velocity profile

$$v(y) = \frac{1}{2\mu} \frac{\partial p}{\partial x} (y^2 - h^2) + v_f \quad (29)$$

where $v = v_f = \mathbf{v}_f \cdot \mathbf{t}_{\Gamma_d}$ has been applied as essential boundary condition on both faces of the cavity ($y = -h$ and $y = h$).

Under the assumption of small changes in the concentrations the balance of mass in the ‘micro’-flow reads

$$\frac{\partial \rho_f}{\partial t} + \rho_f \nabla \cdot \mathbf{v}_m = 0 \quad (30)$$

We now assume that the first term can be neglected because the problem is monophasic in the cavity and the velocities are therefore much higher than in the porous medium. Accordingly, the mass balance inside the cavity simplifies to

$$\nabla \cdot \mathbf{v}_m = 0 \quad (31)$$

If $w = \mathbf{v}_m \cdot \mathbf{n}_{\Gamma_d}$ denotes the normal component of the velocity, the simplified mass balance can also be written as

$$\frac{\partial v}{\partial x} + \frac{\partial w}{\partial y} = 0 \quad (32)$$

To relate the mass transfer equation to the ‘macro’-flow, the mass balance is averaged over the width of the cavity:

$$\int_{-h}^{h} \frac{\partial v}{\partial x} dy + \int_{-h}^{h} \frac{\partial w}{\partial y} dy = 0 \quad (33)$$

Accordingly, the difference in the velocity components that are normal to both crack faces is given by

$$\llbracket w_f \rrbracket = w(h) - w(-h) = - \int_{-h}^{h} \frac{\partial v}{\partial x} dy$$

Substitution of Equation (29) into Equation (33) then leads to

$$\llbracket w_f \rrbracket = \frac{2}{3\mu} \frac{\partial}{\partial x} \left( \frac{\partial p}{\partial x} h^3 \right) - 2h \frac{\partial v_f}{\partial x} \quad (34)$$

This equation gives the amount of fluid attracted in the tangential fluid flow. It can be included in the weak form of the mass balance of the ‘macro’-flow to ensure the coupling between the ‘micro’-flow and the ‘macro’-flow. Indeed, the coupling term $\mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d$ can now be written as

$$\mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d = n_f \llbracket w_f - w_s \rrbracket = n_f \left( \frac{2}{3\mu} \frac{\partial}{\partial x} \left( \frac{\partial p}{\partial x} h^3 \right) - 2h \frac{\partial v_f}{\partial x} - 2 \frac{\partial h}{\partial t} \right) \quad (35)$$
where $w_s$ is the velocity of the solid skeleton in the normal direction, and the difference between the two crack faces is given by

$$ \llbracket w_s \rrbracket = 2 \frac{\partial h}{\partial t} $$

and $v_f$ follows from Darcy's law in the tangential direction

$$ v_f = \left( \mathbf{v}_s - \frac{k_f}{n_f} \nabla p \right) \cdot \mathbf{t}_{\Gamma_d} $$

In Equation (35), the last term ensures the continuity between the cavity and the bulk. The first term is driven by the balance of momentum of the viscous fluid.

# 4. NUMERICAL ELABORATION

## 4.1. Discretization

The model presented in the two previous sections puts some requirements on the displacement and pressure fields in the vicinity of the discontinuity. The opening of the discontinuity must be described with sufficient accuracy in order to reproduce the cavity geometry properly. The pressure field has to be continuous across the discontinuity, but its normal derivative should be discontinuous from one face of the cavity to the other. These requirements can be satisfied by exploiting the partition-of-unity property of finite element shape functions [9].

In the spirit of previous work on cracks in single-phase media [10–12] the interpolation of each component of the displacement field of the solid phase is enriched with discontinuous and singular functions

$$ u_s = \sum_{i \in \mathcal{N}} N_i \bar{u}_i + \sum_{i \in \mathcal{N}_{\text{cut}}} N_i \mathcal{H}_{\Gamma_d} \hat{u}_i + \sum_{i \in \mathcal{N}_{\text{tip}}} \sum_{j \in [1,4]} N_i \Psi_j \tilde{u}_{ij} \quad (36) $$

where $N_i$ are standard finite element shape functions supported by the set of nodes $\mathcal{N}$ included in the discretized domain $\Omega$. Nodes in $\mathcal{N}_{\text{cut}}$ have their support completely cut by the discontinuity. They hold additional degrees of freedom $\hat{u}_i$ corresponding to the discontinuous function $\mathcal{H}_{\Gamma_d}$ defined by

$$ \mathcal{H}_{\Gamma_d}(\mathbf{x}) = \frac{\mathbf{x} \cdot \mathbf{n}_{\Gamma_d}}{\|\mathbf{x} \cdot \mathbf{n}_{\Gamma_d}\|} \quad (37) $$

Nodes whose support contain the tips of the discontinuity are included in $\mathcal{N}_{\text{tip}}$. Four singular functions $\Psi_j$ and their corresponding degrees of freedom $\tilde{u}_{ij}$ are supported by this set of nodes (Figure 4). The singular functions are

$$ \{\Psi_j(r, \theta)\}_{j \in [1,4]} = \left\{ \sqrt{r} \cos \frac{\theta}{2}, \sqrt{r} \cos \frac{\theta}{2} \sin \theta, \sqrt{r} \sin \frac{\theta}{2}, \sqrt{r} \sin \frac{\theta}{2} \sin \theta \right\} \quad (38) $$

where $r, \theta$ are the polar co-ordinates of the tip of the discontinuity. Symbolically, Equation (36) can be written as

$$ \mathbf{u}_s = \mathbf{N}\mathbf{U} \quad (39) $$
<figure><img src="image_5.png" /><figcaption>Figure 4. Enrichment strategy: circles denote $N_{cut}$, square $N_{tip}$, and grey field circle $N_{pres}$.</figcaption></figure>

where the matrix **N** contains the standard interpolation polynomials $N_i$ as well as the discontinuous function $\mathcal{H}_{\Gamma_d}$ and the singular functions $\Psi_j$, and the array **U** contains the displacement degrees of freedom $\bar{u}_i$, $\hat{u}_i$, and $\tilde{u}_{ij}$. Furthermore, we define the matrix $\llbracket \mathbf{N} \rrbracket$ which contains the jumps in the interpolation polynomials, and the matrix **B**, which contains the bounded part of spatial derivatives of **N**.

For the pressure approximation, the standard finite element interpolation is enriched with the function $\mathcal{D}_{\Gamma_d}$, which is the distance to the discontinuity. $\mathcal{D}_{\Gamma_d}$ is continuous through the discontinuity, but its normal derivative is discontinuous and equal to $\mathcal{H}_{\Gamma_d}$

$$ \mathbf{n}_{\Gamma_d} \cdot \nabla \mathcal{D}_{\Gamma_d} = \mathcal{H}_{\Gamma_d} \qquad (40) $$

All nodes whose support is cut, or partially cut by the discontinuity, $N_{pres}$, hold additional pressure degrees of freedom, cf. [12]:

$$ p = \sum_{i \in \mathcal{N}} H_i \bar{p}_i + \sum_{i \in \mathcal{N}_{\text{pres}}} H_i \mathcal{D}_{\Gamma_d} \hat{p}_i \qquad (41) $$

$H_i$ are the finite element shape functions used as partition of unity for the pressure interpolation. In a fashion similar to the displacement interpolation, we write Equation (41) symbolically as

$$ p = \mathbf{H}\mathbf{P} \qquad (42) $$

where **H** contains the shape functions $H_i$ as well as the distance function $\mathcal{D}_{\Gamma_d}$, and the array **P** contains the degrees of freedom $\bar{p}_i$ and $\hat{p}_i$.
The choices for $N_i$ and $H_i$ are driven by modelling requirements. Indeed, the modelling of the fluid flow inside the cavity needs the second derivatives of the pressure, see Equation (35). Hence, the order of the finite element shape functions $H_i$ has to be high enough, otherwise the coupling between the fluid flow in the cavity and the bulk will not be achieved. Further, the order of the finite element shape functions $N_i$ must be greater than or equal to the order of $H_i$ for consistency in the discrete balance of momentum equation. The examples in the next section illustrate these arguments.

## 4.2. Discrete equations and resolution

In a Bubnov-Galerkin sense we choose the test functions for the displacements and the pressures, $\eta$ and $\zeta$, respectively, in the same space as the displacement and pressure interpolations defined by Equations (36) and (41), and require that Equation (19) and Equation (20) hold for all admissible $\eta$ and $\zeta$. This gives:

$$ \int_{\Omega} \mathbf{B}^{\mathrm{T}} \boldsymbol{\sigma} \, \mathrm{d}\Omega + \mathbf{F}_{\mathrm{inter}} = \mathbf{F}_{\mathrm{ext}} \quad (43) $$

$$ - \int_{\Omega} \alpha \mathbf{H}^{\mathrm{T}} \mathbf{m}^{\mathrm{T}} \dot{\nabla} \mathbf{u}_{\mathrm{s}} \, \mathrm{d}\Omega + \int_{\Omega} k_{\mathrm{f}} \nabla \mathbf{H}^{\mathrm{T}} \nabla p \, \mathrm{d}\Omega - \int_{\Omega} Q^{-1} \mathbf{H}^{\mathrm{T}} \dot{p} \, \mathrm{d}\Omega + \mathbf{Q}_{\mathrm{inter}} = \mathbf{Q}_{\mathrm{ext}} \quad (44) $$

where, for two dimensions, $\mathbf{m} = [1, 1, 0]$, and the external force and flux vectors, $\mathbf{F}_{\mathrm{ext}}$ and $\mathbf{Q}_{\mathrm{ext}}$, are given by

$$ \mathbf{F}_{\mathrm{ext}} = \int_{\Gamma} \mathbf{N}^{\mathrm{T}} \mathbf{t}_{\mathrm{p}} \, \mathrm{d}\Gamma \quad (45) $$

$$ \mathbf{Q}_{\mathrm{ext}} = \int_{\Gamma} \mathbf{H}^{\mathrm{T}} \mathbf{n}^{\mathrm{T}} \mathbf{q}_{\mathrm{p}} \, \mathrm{d}\Gamma \quad (46) $$

The time integration is carried out using a backward finite difference scheme

$$ \left( \frac{\mathrm{d}(.)}{\mathrm{d}t} \right)^{t+\Delta t} = \frac{(.)^{t+\Delta t} - (.)^t}{\Delta t} = \frac{\Delta(.)}{\Delta t} \quad (47) $$

where $\Delta t$ is the time increment, while $(.)^t$ and $(.)^{t+\Delta t}$ denote the unknowns at $t$ and $t + \Delta t$, respectively. Substitution in the semi-discrete equation (44) leads to:

$$ - \int_{\Omega} \alpha \mathbf{H}^{\mathrm{T}} \mathbf{m}^{\mathrm{T}} \nabla \Delta \mathbf{u}_{\mathrm{s}} \, \mathrm{d}\Omega + \Delta t \int_{\Omega} k_{\mathrm{f}} \nabla \mathbf{H}^{\mathrm{T}} \nabla p \, \mathrm{d}\Omega - \int_{\Omega} Q^{-1} \mathbf{H}^{\mathrm{T}} \Delta p \, \mathrm{d}\Omega + \Delta t \mathbf{Q}_{\mathrm{inter}} = \Delta t \mathbf{Q}_{\mathrm{ext}} \quad (48) $$

The interfacial force vector $\mathbf{F}_{\mathrm{inter}}$ is derived from Equation (22) by integrating $[\eta] \cdot p \mathbf{n}_{\Gamma_d}$ along $\Gamma_d$:

$$ \mathbf{F}_{\mathrm{inter}} = - \left( \int_{\Gamma_d} [\mathbf{N}]^{\mathrm{T}} \mathbf{n}_{\Gamma_d} \mathbf{H} \, \mathrm{d}\Gamma \right) \mathbf{P} \quad (49) $$
In a similar fashion, Equation (23) is integrated along $\Gamma_d$ to give

$$ \mathbf{Q}_{\text{inter}} = \int_{\Gamma_d} \mathbf{H}^{\text{T}} \mathbf{n}_{\Gamma_d}^{\text{T}} \mathbf{q}_d \, d\Gamma \qquad (50) $$

The interfacial flux vector $\mathbf{Q}_{\text{inter}}$ is then elaborated by incorporating the discrete form of Equation (35), yielding:

$$ \begin{aligned} \mathbf{n}_{\Gamma_d}^{\text{T}} \mathbf{q}_d = & \frac{n_f}{12\mu} (\mathbf{t}_{\Gamma_d}^{\text{T}} (\nabla \nabla \mathbf{H} \mathbf{t}_{\Gamma_d}) \mathbf{P}) (\mathbf{n}_{\Gamma_d}^{\text{T}} \llbracket \mathbf{N} \rrbracket \mathbf{U})^3 \\ & + \frac{n_f}{4\mu} (\mathbf{t}_{\Gamma_d}^{\text{T}} \nabla \mathbf{H} \mathbf{P}) (\mathbf{n}_{\Gamma_d}^{\text{T}} \llbracket \nabla \mathbf{N} \rrbracket \mathbf{t}_{\Gamma_d} \mathbf{U}) (\mathbf{n}_{\Gamma_d}^{\text{T}} \llbracket \mathbf{N} \rrbracket \mathbf{U})^2 \\ & - (\mathbf{n}_{\Gamma_d}^{\text{T}} \llbracket \mathbf{N} \rrbracket \mathbf{U}) (n_f \langle \mathbf{N} \rangle \dot{\mathbf{U}} - k_f \nabla \mathbf{H} \mathbf{P})^{\text{T}} \mathbf{t}_{\Gamma_d} \\ & - n_f \mathbf{n}_{\Gamma_d}^{\text{T}} \llbracket \mathbf{N} \rrbracket \dot{\mathbf{U}} \end{aligned} \qquad (51) $$

where $\langle \cdot \rangle$ is the average of . taken on the — and the + sides of the discontinuity.

As $\mathbf{Q}_{\text{inter}}$ is non-linear, an iterative procedure must be used to compute the solution at each time step. For this purpose, a residual vector $\mathbf{R}^i$ is defined at iteration $i$:

$$ \mathbf{R}^i = \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{K}_{up}^{\text{T}} & \mathbf{K}_{pp}^{(1)} \end{bmatrix} \left( \frac{\Delta \mathbf{U}}{\Delta \mathbf{P}} \right)^i + \begin{bmatrix} \mathbf{K}_{uu} & \mathbf{K}_{up} \\ \mathbf{0} & \Delta t \mathbf{K}_{pp}^{(2)} \end{bmatrix} \left( \begin{matrix} \mathbf{U} \\ \mathbf{P} \end{matrix} \right)^i + \left( \frac{\mathbf{F}_{\text{inter}}}{\Delta t \mathbf{Q}_{\text{inter}}} \right)^i - \left( \frac{\mathbf{F}_{\text{ext}}}{\Delta t \mathbf{Q}_{\text{ext}}} \right) \qquad (52) $$

where the stiffness matrices are defined as

$$ \mathbf{K}_{uu} = \int_{\Omega} \mathbf{B}^{\text{T}} \mathbf{D}^{\text{tan}} \mathbf{B} \, d\Omega \qquad (53) $$

$$ \mathbf{K}_{up} = - \int_{\Omega} \alpha \mathbf{B}^{\text{T}} \mathbf{m} \mathbf{H} \, d\Omega \qquad (54) $$

$$ \mathbf{K}_{pp}^{(1)} = - \int_{\Omega} Q^{-1} \mathbf{H}^{\text{T}} \mathbf{H} \, d\Omega \qquad (55) $$

$$ \mathbf{K}_{pp}^{(2)} = - \int_{\Omega} k_f \nabla \mathbf{H}^{\text{T}} \nabla \mathbf{H} \, d\Omega \qquad (56) $$

In a full Newton-Raphson algorithm, the iterative matrix $\mathbf{K}^i$ is the Jacobian matrix of the residual $\mathbf{R}$

$$ \mathbf{K}^i = \begin{bmatrix} \mathbf{K}_{uu} & \mathbf{K}_{up} + \frac{\partial \mathbf{F}_{\text{inter}}}{\partial \mathbf{P}} \\ \mathbf{K}_{up}^{\text{T}} + \Delta t \frac{\partial \mathbf{Q}_{\text{inter}}}{\partial \mathbf{U}} & \mathbf{K}_{pp}^{(1)} + \Delta t \mathbf{K}_{pp}^{(2)} + \Delta t \frac{\partial \mathbf{Q}_{\text{inter}}}{\partial \mathbf{P}} \end{bmatrix} \qquad (57) $$

where all quantities are evaluated at iteration $i$.
The coupling terms, $F_{inter}$ and $Q_{inter}$, cause the Jacobian matrix to become unsymmetric. To restore symmetry, the contributions due to the coupling terms are omitted in the Jacobian matrix, and the iterations in the example calculations in the next section have been carried out with the following stiffness matrix:

$$ \mathbf{K}^i = \begin{bmatrix} \mathbf{K}_{uu} & \mathbf{K}_{up} \\ \mathbf{K}_{up}^T & \mathbf{K}_{pp}^{(1)} + \Delta t \mathbf{K}_{pp}^{(2)} \end{bmatrix} \qquad (58) $$

## 5. EXAMPLE CALCULATIONS

A two-dimensional specimen under plane-strain conditions is considered. The specimen is a square-shaped fractured block, see Figure 5. A normal fluid flux $q_0 = 10^{-4} \text{ms}^{-1}$ starting at $t=0 \text{ s}$ is imposed at the bottom face while the top face is assigned a drained condition with zero pressure. Both left and right faces have undrained boundary conditions. No mechanical load is applied, but essential boundary conditions have been applied in order to remove rigid body motions. The block is 10m × 10m and consists of a porous material with a fluid volumic fraction $n_f = 0.3$. The absolute mass densities are $\rho'_s = \rho_s/n_s = 2000 \text{ kg/m}^3$ for the solid phase and $\rho'_f = \rho_f/n_f = 1000 \text{ kg/m}^3$ for the fluid phase. The solid constituent is assumed to behave in a linear elastic manner with a Young's modulus $E = 9 \text{ GPa}$ and a Poisson ratio $\nu = 0.4$. The Biot coefficient $\alpha$ has been set equal to 1, and the Biot modulus has been assigned a value $Q = 10^{18} \text{ GPa}$ so as to simulate a quasi-incompressible fluid. This is not a limitation of the model, but more clearly brings out the influence of a fault. The bulk material has a permeability $k_f = 10^{-9} \text{ m}^3/\text{N s}$ while the viscosity of the fluid is $\mu = 10^{-3} \text{ N/m}^2 \text{ s}$.

<figure><img src="image_6.png" /><figcaption>Figure 5. Fractured block with boundary conditions.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. Fine and coarse mesh composed of quadrilateral elements.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Deformed mesh for a crack angle of 30° (magnified × 1000).</figcaption></figure>

### 5.1. One inclined fault

Initially, only one fault is considered. This fault, at the centre of the specimen, is 2 m long, and has an angle $α$ with the horizontal axis. A reference simulation has been run for a period of 10 s using 75 time steps. For this reference simulation, the coarse mesh shown in Figure 6 has been used with quadrilateral elements equipped with quadratic shape functions.

The results obtained for $α = 30°$ are shown in Figures 7–9. Because of the imposed fluid flux at the bottom, the pressure increases in the specimen and inside the fault, which subsequently opens,
<figure><img src="image_9.png" /><figcaption>Figure 8. Normal derivative of the pressure $n_{\Gamma_d} \cdot \nabla p$ for a crack angle of 30°. The top picture, where the data have been plotted as a piecewise-constant field, shows the results for the coarse mesh. Below, the results for the fine mesh have been smoothed and are depicted with filled iso-values.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9. Fluid velocity field for a crack angle of 30°.</figcaption></figure>

and fluid can flow inside, see Figure 7. The discontinuity in the normal derivative of the pressure is illustrated in Figure 8. Its jump has an opposite sign at both tips of the fault. The pore fluid flows into the cavity at the left tip where the gradient of the fault opening is positive. At the right tip, the gradient of the fault opening is negative and the fluid in the cavity flows back into the bulk.

Figure 9 gives the absolute velocity field of the fluid. The ‘macro’-flow is oriented from the bottom to the top of the specimen. The ‘micro’-flow in the cavity is oriented by the direction of the fault. In the cavity, the velocity of the fluid is very high compared to the velocity in the bulk because there is no resisting solid skeleton. As a consequence, the global response of the specimen, which is defined in terms of the ratio outflow over inflow, $q_{out}/q_{in}$, is affected by the presence of the fault. Indeed, part of the fluid in the vicinity of the fault is, at least temporarily, stored inside the cavity, which causes the global flow to slow down in a macroscopic sense. This is illustrated in Figure 10, which shows the ratio $q_{out}/q_{in}$ for different fault angles. These results are compared with those obtained without a fault. The time scale has been normalized by the 5% response time $\tau_{5\%}$ of the non-fractured specimen.

The results of Figure 10 show that the variation of the response time can be up to about 4% for a horizontal fault. Obviously, the variation of the time response decreases when the fault angle increases. Indeed, for the degenerated case where the fault is vertical, the pressure gradient is tangential to the fault. Hence, even if the fault is opening, the tangential flow which is supplied by the discontinuity in the normal derivative of the pressure, is not activated.
<figure><img src="image_11.png" /><figcaption>Figure 10. The ratio of the outward flow over influx for several fault angles. Results are for the coarse mesh with quadrilateral elements and quadratic shape functions.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Comparison of the outward fluid flow for different quadrilateral elements.</figcaption></figure>

Figures 11 and 12 show the influence of the order of interpolation. Figure 11 justifies the use of quadrilateral elements with linear shape functions: the refinement of the mesh allows the use of a lower order of interpolation. In Figure 11, a comparison is made between the results obtained with the reference simulation and a simulation with the fine mesh of Figure 6 and bilinear shape functions. The results are almost identically similar and the latter type of interpolation will therefore be used in the simulation of the next paragraph. Such 'convergence' is not obtained with triangular elements, see Figure 12. Indeed, to be consistent with the model, the interpolation must have non-zero second derivatives. Otherwise, the coupling term between the bulk and the cavity is not
<figure><img src="image_13.png" /><figcaption>Figure 12. Comparison of the outward fluid flow for different triangular elements.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. Three levels of mesh refinement with triangular elements.</figcaption></figure>

computed correctly. With quadrilateral elements, a constant second derivative is obtained by virtue of the bilinear shape functions, which is evidently not the case for triangular elements. Figure 12 proves that mesh refinement, see Figure 13, cannot remedy this deficiency.

## 5.2. *Ten randomly generated faults*

A set of 10 faults is randomly generated, see Figure 14. The length of the faults is between 1 and 3 m, while the fault angle varies from $\alpha = -10^\circ$ to $30^\circ$. Figure 15 shows the influence of the faults on the norm of the pressure gradient. The global fluid flow is strongly affected by the ‘micro’-flows inside the faults. From Figure 15 it is observed that the main effect is due to the two longest faults. Evidently, the effect of a fault on the ‘macro’-flow increases with its length because more fluid can flow inside the cavity. The outward flow obtained for this 10 randomly generated faults is plotted in Figure 16. The variation of the response time is about 18%. For this simulation, the fine quadrilateral mesh has been used with bilinear shape functions. Indeed, when a large number of faults are considered, the iterative matrix becomes ill-conditioned when using quadratic shape functions. This is avoided by using bilinear shape functions.
<figure><img src="image_15.png" /><figcaption>Figure 14. Ten randomly generated faults.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. $L_2$-Norm of the pressure gradient.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16. Outward fluid flow for 10 randomly generated faults.</figcaption></figure>

## 6. CONCLUDING REMARKS

A methodology has been proposed to insert discontinuities such as cracks, faults, or shear bands, in a porous medium. The discontinuities can be located arbitrarily, not related to the underlying discretization. For the fluid flow in the fractured porous medium a two-scale approach has been chosen, where the flow of the fluid inside the discontinuity (the ‘micro’-scale) is modelled independently from the flow of the pore fluid in the surrounding porous medium (the ‘macro’-scale). The mechanical and the mass transfer couplings between the two scales are obtained by inserting the homogenized ‘constitutive’ relations of the ‘micro’-flow into the weak form of the balance equations of the bulk. The assumptions made for the fluid flow in and near the discontinuity require the addition of special enrichment functions for the displacement and the pressure fields. These conditions are satisfied by exploiting the partition of unity property of the finite element polynomial shape functions. The examples confirm the efficiency of the method and also guide the choice of the order of interpolation. Finally, the influence has been studied which faults have on the fluid flow in a porous medium. It has been shown that the effect can be extensive.

## REFERENCES

1. Terzaghi K. *Theoretical Soil Mechanics*. Wiley: New York, 1943.

2. Biot MA. *Mechanics of Incremental Deformations*. Wiley: Chichester, 1965.

3. Lewis RW, Schrefler BA. *The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media* (2nd edn). Wiley: Chichester, 1998.

4. Snijders H, Huyghe JM, Janssen JD. Triphasic finite element model for swelling porous media. *International Journal for Numerical Methods in Fluids* 1995; **20**:1039–1046.

5. Huyghe JM, Janssen JD. Quadriphasic mechanics of swelling incompressible media. *International Journal of Engineering Science* 1997; **35**:793–802.

6. Van Loon R, Huyghe JM, Wijlaars MW, Baaijens FPT. 3D FE implementation of an incompressible quadriphasic mixture model. *International Journal for Numerical Methods in Engineering* 2003; **57**:1243–1258.

7. Jouanna P, Abellan MA. Generalized approach to heterogeneous media. In *Modern Issues in Non-Saturated Soils*, Gens A, Jouanna P, Schrefler B (eds). Springer: New York, 1995; 1–128.
8. de Borst R, Réthoré J, Abellan MA. A numerical approach for arbitrary cracks in a fluid-saturated medium. *Archive of Applied Mechanics* 2006; **75**:595-606.

9. Babuska I, Melenk JM. The partition of unity method. *International Journal for Numerical Methods in Engineering* 1997; **40**:727-758.

10. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. *International Journal for Numerical Methods in Engineering* 1999; **45**:601-620.

11. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 1999; **46**:131-150.

12. Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. *International Journal for Numerical Methods in Engineering* 2001; **50**:993-1013.

13. Wells GN, Sluys LJ. Discontinuous analysis of softening solids under impact loading. *International Journal for Numerical and Analytical Methods in Geomechanics* 2001; **25**:691-709.

14. Wells GN, Sluys LJ, de Borst R. Simulating the propagation of displacement discontinuities in a regularized strain-softening medium. *International Journal for Numerical Methods in Engineering* 2002; **53**:1235-1256.

15. Wells GN, de Borst R, Sluys LJ. A consistent geometrically non-linear approach for delamination. *International Journal for Numerical Methods in Engineering* 2002; **54**:1333-1355.

16. Remmers JJC, de Borst R, Needleman A. A cohesive segments method for the simulation of crack growth. *Computational Mechanics* 2003; **31**:69-77.

17. Samaniego E, Belytschko T. Continuum-discontinuum modelling of shear bands. *International Journal for Numerical Methods in Engineering* 2005; **62**:1857-1872.

18. Réthoré J, Gravouil A, Combescure A. An energy-conserving scheme for dynamic crack growth using the extended finite element method. *International Journal for Numerical Methods in Engineering* 2005; **63**:631-659.

19. Réthoré J, Gravouil A, Combescure A. A combined space-time extended finite element method. *International Journal for Numerical Methods in Engineering* 2005; **64**:260-284.

20. Areias PMA, Belytschko T. Two-scale shear band evolution by local partition of unity. *International Journal for Numerical Methods in Engineering* 2006; **66**:878-910.

21. Abellan MA, de Borst R. Wave propagation and localisation in a softening two-phase medium. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**:5011-5019.