# A precis of two-scale approaches for fracture in porous media

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

► To cite this version:

René de Borst, Julien Réthoré, Marie-Angèle Abellan. A precis of two-scale approaches for fracture in porous media. Lecture Notes on Composite Materials, Springer, pp.149-171, 2008, 10.1007/978-1-4020-8772-1_5. http://dx.doi.org/10.1007/978-1-4020-8772-1

HAL Id: hal-00383662

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

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 PRECIS OF TWO-SCALE APPROACHES
FOR FRACTURE IN POROUS MEDIA

R. de Borst

Department of Mechanical Engineering

Eindhoven University of Technology

NL-5600 MB Eindhoven

The Netherlands &

LaMCoS, UMR CNRS 5514

INSA de Lyon

Villeurbanne, France

r.d.borst@tue.nl

J. Rethore

LMT, E.N.S. de Cachan, Cachan, France

M.-A. Abellan

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

**Abstract**

A derivation is given of two-scale models that are able to describe deformation and flow in a fluid-saturated and progressively fracturing porous medium. From the micromechanics of the flow in the cavity, identities are derived that couple the local momentum and the mass balances to the governing equations for a fluid-saturated porous medium, which are assumed to hold on the macroscopic scale. By exploiting the partition-of-unity property of the finite element shape functions, the position and direction of the fracture is independent from the underlying discretisation. The finite element equations are derived for this two-scale approach and integrated over time. The resulting discrete equations are nonlinear due to the cohesive crack model and the nonlinearity of the coupling terms. A consistent linearisation is given for use within a Newton-Raphson iterative procedure. Finally, examples are given to show the versatility and the efficiency of the approach.

**Keywords:**

multiscale analysis, partition-of-unity approach, fracture, porous media, multi-phase media, cohesive cracks
# 1. Introduction

Since the pioneering work of Terzaghi [1] and Biot [2] the flow of fluids in deforming porous media has received considerable attention. Recently, Lewis and Schrefler [3] have given an account of this topic which is crucial for understanding and predicting the physical behaviour of many systems of interest, for example, in geotechnical and petroleum engineering, but also for soft tissues. Because of the complicated structure and functioning of human tissues, the classical two-phase theory has been extended to three and four-phase media, taking into account ion transport and electrical charges [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, and shear bands, can markedly change the physical behaviour. Furthermore, the fluid can transport contaminants which can dramatically reduce the strength of the solid skeleton. To account for such phenomena, the fluid flow must be studied also in the presence of discontinuities in the solid phase. The physics of the flow within such discontinuities can be very different from that of the interstitial fluid in the deforming bulk material. These differences affect the flow pattern and therefore also the deformations in the vicinity of the discontinuity. As we will show at the end of the paper, the local differences in flow characteristics can even influence the flow and deformations in the entire body of interest.

In this contribution, we will describe a general numerical methodology to capture deformation and flow in progressively fracturing porous media, summarising and unifying the recent work reported in [8–11]. The model allows for flow inside an evolving crack to be in the tangential direction. This is achieved by a priori adopting a two-scale approach. At the fine scale the flow in the cavity created by the (possibly cohesive) crack can be modelled in various ways, e.g. as a Stokes flow in an open cavity, or using a Darcy relation for a damaged porous material. Since the cross-sectional dimension of the cavity is small compared to its length, the flow equations can be averaged over the width 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, at the macroscopic scale, is assumed to hold for the surrounding
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 discretisations, these requirements can be satisfied by exploiting the partition-of-unity property of finite element shape functions [12], as has been done successfully in applications to cracking in single-phase media [13–23].

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 formulation is crucial for incorporating the micromechanical flow model properly. This micromechanical flow model is discussed in the next section, where it will be demonstrated how the momentum and mass couplings of the micromechanical flow model to the surrounding porous medium can be accomplished in the weak formulation. Time integration and a consistent linearisation of the resulting equations, which are nonlinear due to the coupling terms and the cohesive crack model complete the numerical model. Finally, example calculations are given of a body with stationary and propagating cracks.

## 2. Balance equations

We consider a two-phase medium subject to the restriction of small displacement gradients and small variations in the concentrations [7]. Furthermore, the assumptions are made that there is no mass transfer between the constituents and that the processes which we consider, occur isothermally. With these assumptions, the balances of linear momentum for the solid and the fluid phases read

$$ \nabla \cdot \boldsymbol{\sigma}_{\pi} + \hat{\mathbf{p}}_{\pi} + \rho_{\pi} \mathbf{g} = \frac{\partial(\rho_{\pi} \mathbf{v}_{\pi})}{\partial t} + \nabla \cdot (\rho_{\pi} \mathbf{v}_{\pi} \otimes \mathbf{v}_{\pi}) \quad (1) $$

with $\boldsymbol{\sigma}_{\pi}$ the stress tensor, $\rho_{\pi}$ the apparent mass density, and $\mathbf{v}_{\pi}$ the absolute velocity of constituent $\pi$. As in the remainder of this paper, $\pi = s, f$, with $s$ and $f$ denoting the solid and fluid phases, respectively. Further, $\mathbf{g}$ is the gravity acceleration and $\hat{\mathbf{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{\mathbf{p}}_{\pi} = 0 \quad (2) $$

We now neglect convective, gravity and acceleration terms, so that the momentum balances reduce to

$$ \nabla \cdot \boldsymbol{\sigma}_{\pi} + \hat{\mathbf{p}}_{\pi} = \mathbf{0} \quad (3) $$
Adding both momentum balances, and taking into account Eq. (2), one obtains the momentum balance for the mixture:

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

where the stress is, as usual, composed of a solid and a fluid part, $\boldsymbol{\sigma} = \boldsymbol{\sigma}_s + \boldsymbol{\sigma}_f$.

In a similar fashion as for the balances of momentum, one can write the balance of mass for each phase as

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

Again neglecting convective terms, the mass balances can be simplified to give

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

We multiply the mass balance for each constituent $\pi$ by its volumetric ratio $n_{\pi}$, add them and utilise the constraint

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

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 (8) $$

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}{K_t} \frac{n_s}{\rho_s} \frac{\partial \rho_s}{\partial t} \qquad (9) $$

with $K_s$ the bulk modulus of the solid material and $K_t$ the overall bulk modulus of the porous medium. Using 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 (10) $$

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

$$ \frac{1}{Q} dp = \frac{n_f}{\rho_f} d\rho_f \qquad (11) $$

with the overall compressibility, or Biot modulus

$$ \frac{1}{Q} = \frac{\alpha - n_f}{K_s} + \frac{n_f}{K_f} \qquad (12) $$
where $K_f$ is the bulk modulus of the fluid. Inserting relations (10) and (11) into the balance of mass of the total medium, Eq. (8), 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 \quad (13) $$

The field equations, i.e. the balance of momentum of the saturated medium, Eq. (4), and the balance of mass, Eq. (13), are complemented by the boundary conditions

$$ \mathbf{n}_{\Gamma} \cdot \boldsymbol{\sigma} = \mathbf{t}_{p}, \quad \mathbf{v}_{s} = \mathbf{v}_{p} \quad (14) $$

which hold on complementary parts of the boundary $\partial\Omega_t$ and $\partial\Omega_v$, with $\Gamma = \partial\Omega = \partial\Omega_t \cup \partial\Omega_v$, $\partial\Omega_t \cap \partial\Omega_v = \emptyset$, $\mathbf{t}_p$ being the prescribed external traction and $\mathbf{v}_p$ the prescribed velocity, and

$$ n_f(\mathbf{v}_f - \mathbf{v}_s) \cdot \mathbf{n}_{\Gamma} = q_p, \quad p = p_p \quad (15) $$

which 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$, $q_p$ and $p_p$ being the prescribed outflow of pore fluid and the prescribed pressure, respectively. The initial conditions specify the displacements $\mathbf{u}_{\pi}$, the velocities $\mathbf{v}_{\pi}$, and the pressure field at $t = 0$,

$$ \mathbf{u}_{\pi}(\mathbf{x}, 0) = \mathbf{u}_{\pi}^{0}, \quad \mathbf{v}_{\pi}(\mathbf{x}, 0) = \mathbf{v}_{\pi}^{0}, \quad p(\mathbf{x}, 0) = p^{0} \quad (16) $$

### 3. Constitutive equations

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

$$ d\boldsymbol{\sigma}'_s = \bar{\mathbf{D}}^{\mathrm{tan}} : d\boldsymbol{\epsilon}_s \quad (17) $$

where $\bar{\mathbf{D}}^{\mathrm{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}^{\mathrm{tan}} : d\boldsymbol{\epsilon}_s \quad (18) $$

where the notation $\mathbf{D}^{\mathrm{tan}} = n_s \bar{\mathbf{D}}^{\mathrm{tan}}$ has been used. In the examples, a linear-elastic behaviour of the bulk material has been assumed, and we have set $\mathbf{D}^{\mathrm{tan}} = \mathbf{D}$, the linear-elastic stiffness tensor.

At the discontinuity $\Gamma_d$ a discrete relation holds between the interface tractions $\mathbf{t}_d$ and the relative displacements $\boldsymbol{\delta}$:

$$ \mathbf{t}_d = \mathbf{t}_d(\boldsymbol{\delta}, \kappa) \quad (19) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. Decohesion curves for a ductile (left) and a quasi-brittle solid (right)</figcaption></figure>

with $\kappa$ a history parameter. After linearisation, necessary to use a tangential stiffness matrix in an incremental-iterative solution procedure, one obtains

$$ \dot{\mathbf{t}}_d = \mathbf{T}\dot{\mathbf{\delta}} \qquad (20) $$

with $\mathbf{T}$ the material tangent stiffness matrix of the discrete traction-separation law:

$$ \mathbf{T} = \frac{\partial \mathbf{t}_d}{\partial \mathbf{\delta}} + \frac{\partial \mathbf{t}_d}{\partial \kappa} \frac{\partial \kappa}{\partial \mathbf{\delta}} \qquad (21) $$

A key element is the presence of a work of separation or fracture energy, $\mathcal{G}_c$, which governs crack growth and enters the interface constitutive relation (19) in addition to the tensile strength $f_t$. It is defined as the work needed to create a unit area of fully developed crack:

$$ \mathcal{G}_c = \int_{\delta_n=0}^{\infty} \sigma d\delta_n \qquad (22) $$

with $\sigma$ the stress across the fracture process zone. It thus equals the area under the decohesion curves as shown in Fig. 1. Evidently, cohesive-zone models as defined above are equipped with an internal length scale, since the quotient $\mathcal{G}_c/E$, with $E$ a stiffness modulus for the surrounding continuum, has the dimension of length.

In a standard manner we adopt Darcy's relation for the flow in the bulk of the porous medium,

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

with $k_f$ the permeability coefficient. The equations for the flow in the cavity, which close the initial value problem, will be detailed in a subsequent section.

## 4. Weak form of the balance equations

To arrive at the weak form of the balance equations, we multiply the momentum balance (4) and the mass balance (13) by kinetically admissible test functions for the displacements of the skeleton, $\eta$, and for the pressure, $\zeta$. After
substitution of Darcy's relation for the flow in the porous medium, Eq. (23), integrating over the domain $\Omega$ and using the divergence theorem then 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 (24) $$

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 (25) $$

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 momentum and mass couplings between the discontinuity – the microscopic scale – and the surrounding porous medium – the macroscopic scale.

The momentum coupling stems from the tractions across the faces of the discontinuity and the pressure applied by the fluid in the discontinuity onto the faces of the discontinuity. We assume stress continuity from the cavity to the bulk, so that we have

$$ \boldsymbol{\sigma} \cdot \mathbf{n}_{\Gamma_d} = \mathbf{t}_d - p \mathbf{n}_{\Gamma_d} \quad (26) $$

with $\mathbf{t}_d$ the cohesive tractions, which are given by Eq. (19). Therefore, 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 (\mathbf{t}_d - p \mathbf{n}_{\Gamma_d}) \, d\Gamma = \int_{\Gamma} \boldsymbol{\eta} \cdot \mathbf{t}_p \, d\Gamma \quad (27) $$

Since the tractions have a unique value across the discontinuity, the pressure $p$ must have the same value at both faces of the discontinuity, and, consequently, this must also hold for the test function for the pressure, $\zeta$. Accordingly, the mass transfer coupling term for the water can be rewritten as follows:

$$ \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} \zeta \mathbf{n}_{\Gamma_d} \cdot 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 (28) $$

where $\mathbf{q}_d = n_f [\mathbf{v}_f - \mathbf{v}_s]$ represents the fluid flux through the faces of the discontinuity.

The above identity 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 cavity 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 crack.

## 5. Micro-macro coupling

To quantify the influence of the “micro”-flow inside the discontinuity on the “macro”-scale, the balances of mass and momentum at the microscale must be considered. Different assumptions can now be made for the cavity. First, we will assume an open cavity which is totally filled with a Newtonian fluid. Then, the balance of mass for the “micro”-flow in the cavity reads

$$ \dot{\rho}_f + \rho_f \nabla \cdot \mathbf{v} = 0 $$

subject to the assumptions of small changes in the concentrations and that convective terms can be neglected, cf. Eq. (6). We 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. With this assumption, and focusing on a two-dimensional configuration, the mass balance inside the cavity simplifies to

$$ \frac{\partial v}{\partial x} + \frac{\partial w}{\partial y} = 0 \qquad (29) $$

with $v = \mathbf{v} \cdot \mathbf{t}_{\Gamma_d}$ and $w = \mathbf{v} \cdot \mathbf{n}_{\Gamma_d}$ the tangential and normal components of the fluid velocity in the discontinuity, respectively, see Fig. 2. Accordingly, the difference in the fluid velocity components that are normal to both crack faces is given by

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

<figure><img src="image_3.png" /><figcaption>Figure 2. Geometry and local coordinate system in cavity</figcaption></figure>
To proceed, the velocity profile of the fluid flow inside the discontinuity must be known. From the balance of momentum for the fluid in the cavity and the assumption of a Newtonian fluid, the following velocity profile results

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

where an integration has been carried out from $y = -h$ to $y = h$ and $\mu$ the viscosity of the fluid. The essential boundary $v = v_f$ has been applied at both faces of the cavity, and stems from the relative fluid velocity in the porous medium at $y = \pm h$:

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

Substitution of Eq. (31) into Eq. (30) and again integrating with respect to $y$ 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 (32)$$

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. Since the difference in the normal velocity of both crack faces is given by $\llbracket w_s \rrbracket = 2 \frac{\partial h}{\partial t}$, the mass coupling term becomes

$$\begin{aligned} \mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d &= n_f \llbracket w_f - w_s \rrbracket \\ &= n_f \left( \frac{2h^3}{3\mu} \frac{\partial^2 p}{\partial x^2} + \frac{2h^2}{\mu} \frac{\partial p}{\partial x} \frac{\partial h}{\partial x} - 2h \frac{\partial v_f}{\partial x} - 2 \frac{\partial h}{\partial t} \right) \end{aligned} \quad (33)$$

Another possibility for closure of the initial value problem is to assume that the cavity is partially filled with solid material. The mass balance for the fluid inside the cavity then reads

$$\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$$

Because the width of the cavity $2h$ is negligible compared to its length, the mass balance is again enforced in an average sense over the cross section. For the first term, we obtain

$$\int_{-h}^{h} \alpha \nabla \cdot \mathbf{v}_s dy = \int_{-h}^{h} \alpha \left( \frac{\partial v_s}{\partial x} + \frac{\partial w_s}{\partial y} \right) dy = \int_{-h}^{h} \alpha \frac{\partial v_s}{\partial x} dy - \alpha \llbracket w_s \rrbracket \quad (34)$$

where $v_s$ and $w_s$ are the component of the solid velocity tangential and normal to the crack, respectively. Because $\alpha$ depends on the capillarity pressure only, it can be assumed as constant over the cross section. Assuming furthermore that
$v_s$ varies linearly with $y$, and defining $\langle v_s \rangle = \frac{1}{2}(v_s(h) + v_s(-h))$, the integral can be solved analytically

$$ \int_{-h}^{h} \alpha \frac{\partial v_s}{\partial x} dy = 2\alpha h \left\langle \frac{\partial v_s}{\partial x} \right\rangle \quad (35) $$

Applying the same operations to the second term of Eq. (30), the following expression ensues

$$ \int_{-h}^{h} n_f \nabla \cdot (\mathbf{v}_f - \mathbf{v}_s) dy = n_f [w_f - w_s] + \frac{\partial}{\partial x} \left( \int_{-h}^{h} n_f (v_f - v_s) dy \right) - \left( n_f (v_f(h) - v_s(h)) \frac{\partial h}{\partial x} - n_f (v_f(-h) - v_s(-h)) \frac{\partial (-h)}{\partial x} \right) \quad (36) $$

We now introduce Darcy's relation projected onto the axis tangential to the crack,

$$ n_f(v_f - v_s) = -k_d \frac{\partial p}{\partial x} \quad (37) $$

with $k_d$ the permeability of the damaged, porous material inside the cavity. The dependence of the permeability inside the cavity, $k_d$, on $y$ will be negligible compared to that on $x$ and accordingly, $k_d$ can be assumed not to depend on $y$. However, the decohesion inside the cavity will affect the permeability, and therefore $k_d = k_d(h)$. Upon substitution of Eq. (37) into Eq. (36), the following relation is obtained

$$ \int_{-h}^{h} n_f \nabla \cdot (\mathbf{v}_f - \mathbf{v}_s) dy = n_f [w_f - w_s] - \frac{\partial}{\partial x} \left( \int_{-h}^{h} k_d \frac{\partial p}{\partial x} dy \right) + \left( k_d \frac{\partial p(h)}{\partial x} \frac{\partial h}{\partial x} - k_d \frac{\partial p(-h)}{\partial x} \frac{\partial (-h)}{\partial x} \right) \quad (38) $$

At the faces of the cavity the permeability equals that of the bulk, and for continuity reasons $\frac{\partial p(h)}{\partial x} = \frac{\partial p(-h)}{\partial x}$, so that the second term becomes

$$ \int_{-h}^{h} n_f \nabla \cdot (\mathbf{v}_f - \mathbf{v}_s) dy = n_f [w_f - w_s] - 2h \frac{\partial k_d(h)}{\partial x} \frac{\partial p}{\partial x} - 2k_d h \frac{\partial^2 p}{\partial x^2} \quad (39) $$

For the third term, neglecting variations of the pressure over the height of the cavity, one obtains

$$ \int_{-h}^{h} \frac{1}{Q} \frac{\partial p}{\partial t} dy = \frac{2h}{Q} \frac{\partial p}{\partial t} \quad (40) $$
Recalling that the term $n_f[\overline{w_f - w_s}]$ can be identified as the coupling term $\mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d$ of the weak form of the mass balance, we obtain

$$ \mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d = -\frac{2h}{Q} \frac{\partial p}{\partial t} + 2\alpha \frac{\partial h}{\partial t} - 2\alpha h \left\langle \frac{\partial v_s}{\partial x} \right\rangle + 2h \frac{\partial k_d(h)}{\partial x} \frac{\partial p}{\partial x} + 2k_d h \frac{\partial^2 p}{\partial x^2} \quad (41) $$

## 6. Discontinuities in a two-phase medium

A finite element method that can accommodate the propagation of discontinuities through elements was proposed by Belytschko and co-workers [13,14], exploiting the partition-of-unity property of finite element shape functions [12]. Since finite element shape functions $\varphi_j$ form partitions of unity, $\sum_{j=1}^n \varphi_j = 1$ with $n$ the number of nodal points, the components $v_i$ of a velocity field $\mathbf{v}$ can be interpolated as

$$ v_i = \sum_{j=1}^{n} \varphi_j \left( \dot{a}_j + \sum_{k=1}^{m} \psi_k \dot{a}_{jk} \right) \quad (42) $$

with $\bar{a}_j$ the “regular” nodal degrees-of-freedom for the displacements, $\psi_k$ the enhanced basis terms, and $\hat{a}_{jk}$ the additional displacement degrees-of-freedom at node $j$ which represent the amplitude of the $k$th enhanced basis term $\psi_k$. Next, we consider a domain $\Omega$ that is crossed by a single discontinuity at $\Gamma_d$ (see Fig. 3). The velocity field $\mathbf{v}_s$ can then be written as the sum of two continuous velocity fields $\bar{\mathbf{v}}_s$ and $\hat{\mathbf{v}}_s$:

$$ \mathbf{v}_s = \bar{\mathbf{v}}_s + \mathcal{H}_{\Gamma_d} \hat{\mathbf{v}}_s \quad (43) $$

where $\mathcal{H}_{\Gamma_d}$ is the Heaviside step function centred at the discontinuity. The decomposition in Eq. (43) has a structure similar to the interpolation in Eq. (42). Accordingly, the partition-of-unity property of finite element shape functions enables the direct incorporation of discontinuities, including cracks and shear bands, in finite element models such that the discontinuous character of cracks

<figure><img src="image_4.png" /><figcaption>Figure 3. Body composed of continuous displacement fields at each side of the discontinuity $\Gamma_d$</figcaption></figure>
and shear bands is preserved. With the standard small-strain assumption that the strain-rate field of the solid, $\dot{\epsilon}_s$, is derived from the symmetric part of the gradient of the velocity field, we obtain

$$ \dot{\epsilon}_s = \nabla^s \bar{\mathbf{v}}_s + \mathcal{H}_{\Gamma_d} \nabla^s \hat{\mathbf{v}}_s \quad (44) $$

away from the discontinuity $\Gamma_d$, with the superscript $s$ denoting the symmetric part of the gradient operator.

With respect to the interpolation of the pressure $p$ we note that the fluid flow normal to the discontinuity can be discontinuous. Since the fluid velocity is related to the pressure gradient via Darcy's relation, the gradient of the pressure normal to the discontinuity can therefore also be discontinuous across the discontinuity. Accordingly, the enrichment of the interpolation of the pressure must be such that the pressure itself is continuous, but has a discontinuous normal first spatial derivative. The distance function $\mathcal{D}_{\Gamma_d}$ defined as

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

satisfies this requirement, and accordingly, the interpolation of $p$ will be such that

$$ p = \bar{p} + \mathcal{D}_{\Gamma_d} \hat{p} \quad (46) $$

We now discretise the fields $\mathbf{v}_s$ and $p$ and the test functions $\eta$ and $\zeta$ for the velocity and the pressure, respectively, in a Bubnov-Galerkin sense:

$$ \begin{aligned} \mathbf{v}_s &= \mathbf{N}(\dot{\mathbf{a}} + \mathcal{H}_{\Gamma_d} \dot{\mathbf{a}}), & \mathbf{\eta} &= \mathbf{N}(\bar{\mathbf{w}} + \mathcal{H}_{\Gamma_d} \hat{\mathbf{w}}) \\ p &= \mathbf{H}(\bar{\mathbf{p}} + \mathcal{D}_{\Gamma_d} \hat{\mathbf{p}}), & \mathbf{\zeta} &= \mathbf{H}(\bar{\mathbf{z}} + \mathcal{D}_{\Gamma_d} \hat{\mathbf{z}}) \end{aligned} \quad (47) $$

where the matrix $\mathbf{N}$ contains the shape functions $N_i$ used as partition of unity for the interpolation of the velocity field $\mathbf{v}_s$, and $\mathbf{H}$ contains the shape functions $H_i$ used as partition of unity for the interpolation of the pressure field $p$. $\dot{\mathbf{a}}$ and $\dot{\mathbf{a}}$ are the nodal arrays assembling the amplitudes that correspond to the standard and enhanced interpolations of the velocity field, and $\bar{\mathbf{p}}$ and $\hat{\mathbf{p}}$ assemble the amplitudes that correspond to the standard and enhanced interpolations of the pressure field. The choice for the interpolants in $\mathbf{N}$ and $\mathbf{H}$ is driven by modelling requirements. Indeed, the modelling of the fluid flow inside the cavity requires the computation of second derivatives of the pressure, see Eqs. (33) or (41). Hence, the order of the finite element shape functions $H_i$ has to be sufficiently high, 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 that of $H_i$ for consistency in the discrete momentum balance.
Inserting Eq. (47) into Eqs. (27) and (28) and requiring that the result holds for all admissible $\bar{w}$, $\bar{z}$, $\hat{w}$ and $\hat{z}$ gives

$$ 
\begin{aligned}
\mathbf{f}_{\bar{a}}^{\text{int}} &= \mathbf{f}_{\bar{a}}^{\text{ext}} \\
\mathbf{f}_{\bar{a}}^{\text{int}} &= \mathbf{f}_{\bar{a}}^{\text{ext}} \\
\mathbf{f}_{\bar{p}}^{\text{int}} &= \mathbf{f}_{\bar{p}}^{\text{ext}} \\
\mathbf{f}_{\hat{p}}^{\text{int}} &= \mathbf{f}_{\hat{p}}^{\text{ext}}
\end{aligned}
\qquad (48) $$

with the external force vectors:

$$ 
\begin{aligned}
\mathbf{f}_{\bar{a}}^{\text{ext}} &= \int_{\Gamma} \mathbf{N}^{\text{T}} \mathbf{t}_{p} \mathrm{d}\Gamma, & \mathbf{f}_{\bar{a}}^{\text{ext}} &= \int_{\Gamma} \mathcal{H}_{\Gamma_d} \mathbf{N}^{\text{T}} \mathbf{t}_{p} \mathrm{d}\Gamma \\
\mathbf{f}_{\bar{p}}^{\text{ext}} &= \int_{\Gamma} \mathbf{H}^{\text{T}} \mathbf{n}_{\Gamma}^{\text{T}} \mathbf{q}_{p} \mathrm{d}\Gamma, & \mathbf{f}_{\hat{p}}^{\text{ext}} &= \int_{\Gamma} \mathcal{D}_{\Gamma_d} \mathbf{H}^{\text{T}} \mathbf{n}_{\Gamma}^{\text{T}} \mathbf{q}_{p} \mathrm{d}\Gamma
\end{aligned}
\qquad (49) $$

and the internal force vectors:

$$ 
\begin{aligned}
\mathbf{f}_{\bar{a}}^{\text{int}} &= \int_{\Omega} \nabla \mathbf{N}^{\text{T}} \boldsymbol{\sigma} \mathrm{d}\Omega \\
\mathbf{f}_{\bar{a}}^{\text{int}} &= \int_{\Omega} \mathcal{H}_{\Gamma_d} \nabla \mathbf{N}^{\text{T}} \boldsymbol{\sigma} \mathrm{d}\Omega + \int_{\Gamma_d} \mathbf{N}^{\text{T}} (\mathbf{t}_d - p\mathbf{n}_{\Gamma_d}) \mathrm{d}\Gamma \\
\mathbf{f}_{\bar{p}}^{\text{int}} &= - \int_{\Omega} \alpha \mathbf{H}^{\text{T}} \mathbf{m}^{\text{T}} \mathbf{v}_s \mathrm{d}\Omega + \int_{\Omega} k_f \nabla \mathbf{H}^{\text{T}} \nabla p \mathrm{d}\Omega - \int_{\Omega} Q^{-1} \mathbf{H}^{\text{T}} \dot{p} \mathrm{d}\Omega \\
\mathbf{f}_{\hat{p}}^{\text{int}} &= - \int_{\Omega} \mathcal{D}_{\Gamma_d} \alpha \mathbf{H}^{\text{T}} \mathbf{m}^{\text{T}} \mathbf{v}_s \mathrm{d}\Omega + \int_{\Omega} k_f \nabla (\mathcal{D}_{\Gamma_d} \mathbf{H})^{\text{T}} \nabla p \mathrm{d}\Omega \\
&\quad - \int_{\Omega} \mathcal{D}_{\Gamma_d} Q^{-1} \mathbf{H}^{\text{T}} \dot{p} \mathrm{d}\Omega + \int_{\Gamma_d} \mathbf{H}^{\text{T}} \mathbf{n}_{\Gamma}^{\text{T}} \mathbf{q}_d \mathrm{d}\Gamma
\end{aligned}
\qquad (50) $$

where, for two dimensions, $\mathbf{m}^{\text{T}} = [1, 1, 0]$.

To carry out the time integration in (54) a backward finite difference scheme is adopted

$$ \left(\frac{d(\cdot)}{dt}\right)^{t+\Delta t} = \frac{(\cdot)^{t+\Delta t} - (\cdot)^t}{\Delta t} = \frac{\Delta(\cdot)}{\Delta t} \qquad (51) $$

where $\Delta t$ is the time increment, while $(\cdot)^t$ and $(\cdot)^{t+\Delta t}$ denote the unknowns at $t$ and $t + \Delta t$, respectively. Substitution of this identity into the set of semi-discrete Eq. (48) results in:

$$ 
\begin{aligned}
\mathbf{f}_{\bar{a}}^{\text{int},t+\Delta t} &= \mathbf{f}_{\bar{a}}^{\text{ext},t+\Delta t} \\
\mathbf{f}_{\bar{a}}^{\text{int},t+\Delta t} &= \mathbf{f}_{\bar{a}}^{\text{ext},t+\Delta t} \\
\mathbf{f}_{\bar{p}}^{\text{int},t+\Delta t} &= \Delta t \mathbf{f}_{\bar{p}}^{\text{ext},t+\Delta t} \\
\mathbf{f}_{\hat{p}}^{\text{int},t+\Delta t} &= \Delta t \mathbf{f}_{\hat{p}}^{\text{ext},t+\Delta t}
\end{aligned}
\qquad (52) $$
where

$$ \mathbf{f}_{\bar{p}}^{\text{int},t+\Delta t} = - \int_{\Omega} \alpha \mathbf{H}^{\text{T}} \mathbf{m}^{\text{T}} \mathbf{u}_s d\Omega + \int_{\Omega} k_f \nabla \mathbf{H}^{\text{T}} \nabla p d\Omega \Delta t \\ - \int_{\Omega} Q^{-1} \mathbf{H}^{\text{T}} \Delta p d\Omega $$

$$ \mathbf{f}_{\hat{p}}^{\text{int},t+\Delta t} = - \int_{\Omega} \mathcal{D}_{\Gamma_d} \alpha \mathbf{H}^{\text{T}} \mathbf{m}^{\text{T}} \mathbf{u}_s d\Omega + \int_{\Omega} k_f \nabla (\mathcal{D}_{\Gamma_d} \mathbf{H})^{\text{T}} \nabla p d\Omega \Delta t \\ - \int_{\Omega} \mathcal{D}_{\Gamma_d} Q^{-1} \mathbf{H}^{\text{T}} \Delta p d\Omega + \int_{\Gamma_d} \mathbf{H}^{\text{T}} \mathbf{n}_{\Gamma}^{\text{T}} \mathbf{q}_d d\Gamma \Delta t $$

while the other internal force vectors remain unchanged except that they are now evaluated at $t + \Delta t$.

For use in a Newton-Raphson solution method, the set (52) must be linearised. To this end, the stress and the pressure are decomposed as

$$ \sigma_j = \sigma_{j-1} + d\sigma, \quad p_j = p_{j-1} + dp \quad (53) $$

with the subscripts $j-1$ and $j$ signifying the iteration numbers. The linearisation of the set (52) then yields

$$ \begin{bmatrix} \mathbf{K}_{\bar{a}\bar{a}} & \mathbf{K}_{\bar{a}\hat{a}} & \mathbf{C}_{\bar{a}\bar{p}} & \mathbf{C}_{\bar{a}\hat{p}} \\ \mathbf{K}_{\hat{a}\hat{a}} & \mathbf{K}_{\hat{a}\hat{a}} & \mathbf{C}_{\hat{a}\bar{p}} & \mathbf{C}_{\hat{a}\bar{p}} \\ \mathbf{C}_{\bar{a}\bar{p}}^{\text{T}} & \mathbf{C}_{\hat{a}\hat{p}}^{\text{T}} & \mathbf{M}_{\bar{p}\bar{p}} + \Delta t \mathbf{K}_{\bar{p}\bar{p}} & \mathbf{M}_{\hat{p}\hat{p}} + \Delta t \mathbf{K}_{\hat{p}\hat{p}} \\ \mathbf{C}_{\hat{a}\bar{p}}^{\text{T}} + \Delta t \mathbf{C}_{\bar{a}\hat{p}}^{\text{T}} & \mathbf{C}_{\hat{a}\bar{p}}^{\text{T}} + \Delta t \mathbf{C}_{\hat{a}\hat{p}}^{\text{T}} & \mathbf{M}_{\bar{p}\bar{p}} + \Delta t \mathbf{K}_{\bar{p}\bar{p}}^{\text{T}} & \mathbf{M}_{\hat{p}\hat{p}} + \Delta t \mathbf{K}_{\hat{p}\hat{p}}^{\text{T}} \end{bmatrix} \begin{pmatrix} d\bar{a} \\ d\hat{a} \\ d\bar{p} \\ d\hat{p} \end{pmatrix} = \times \begin{pmatrix} \mathbf{f}_{\bar{a}}^{\text{ext},t+\Delta t} - \mathbf{f}_{\bar{a},j-1}^{\text{int},t+\Delta t} \\ \mathbf{f}_{\hat{a}}^{\text{ext},t+\Delta t} - \mathbf{f}_{\hat{a},j-1}^{\text{int},t+\Delta t} \\ \Delta t \mathbf{f}_{\bar{p}}^{\text{ext},t+\Delta t} - \mathbf{f}_{\bar{p},j-1}^{\text{int},t+\Delta t} \\ \Delta t \mathbf{f}_{\hat{p}}^{\text{ext},t+\Delta t} - \mathbf{f}_{\hat{p},j-1}^{\text{int},t+\Delta t} \end{pmatrix} \quad (54) $$

with the stiffness matrices:

$$ \mathbf{K}_{\bar{a}\bar{a}} = \int_{\Omega} \nabla \mathbf{N}^{\text{T}} \mathbf{D} \nabla \mathbf{N} d\Omega $$

$$ \mathbf{K}_{\bar{a}\hat{a}} = \int_{\Omega} \mathcal{H}_{\Gamma_d} \nabla \mathbf{N}^{\text{T}} \mathbf{D} \nabla \mathbf{N} d\Omega $$

$$ \mathbf{K}_{\hat{a}\hat{a}} = \int_{\Omega} \mathcal{H}_{\Gamma_d}^2 \nabla \mathbf{N}^{\text{T}} \mathbf{D} \nabla \mathbf{N} d\Omega + \int_{\Gamma_d} \mathbf{N}^{\text{T}} \mathbf{T} \mathbf{N} d\Gamma $$

$$ \mathbf{K}_{\bar{p}\bar{p}} = - \int_{\Omega} k_f \nabla \mathbf{H}^{\text{T}} \nabla \mathbf{H} d\Omega $$
$$ 
\begin{aligned}
\mathbf{K}_{\bar{p}\bar{p}} &= -\int_{\Omega} k_f \nabla \mathbf{H}^{\mathrm{T}} \nabla (\mathcal{D}_{\Gamma_d} \mathbf{H}) \mathrm{d}\Omega \\
\mathbf{K}_{\hat{p}\hat{p}} &= -\int_{\Omega} k_f \nabla (\mathcal{D}_{\Gamma_d} \mathbf{H})^{\mathrm{T}} \nabla (\mathcal{D}_{\Gamma_d} \mathbf{H}) \mathrm{d}\Omega
\end{aligned}
\qquad (55) $$

the mass matrices:

$$ 
\begin{aligned}
\mathbf{M}_{\bar{p}\bar{p}} &= -\int_{\Omega} Q^{-1} \mathbf{H}^{\mathrm{T}} \mathbf{H} \mathrm{d}\Omega \\
\mathbf{M}_{\hat{p}\hat{p}} &= -\int_{\Omega} \mathcal{D}_{\Gamma_d} Q^{-1} \mathbf{H}^{\mathrm{T}} \mathbf{H} \mathrm{d}\Omega \\
\mathbf{M}_{\hat{p}\bar{p}} &= -\int_{\Omega} \mathcal{D}_{\Gamma_d}^2 Q^{-1} \mathbf{H}^{\mathrm{T}} \mathbf{H} \mathrm{d}\Omega
\end{aligned}
\qquad (56) $$

and the coupling matrices:

$$ 
\begin{aligned}
\mathbf{C}_{\bar{a}\bar{p}} &= -\int_{\Omega} \alpha \nabla \mathbf{N}^{\mathrm{T}} \mathbf{m} \mathbf{H} \mathrm{d}\Omega \\
\mathbf{C}_{\hat{a}\hat{p}} &= -\int_{\Omega} \alpha \mathcal{D}_{\Gamma_d} \nabla \mathbf{N}^{\mathrm{T}} \mathbf{m} \mathbf{H} \mathrm{d}\Omega \\
\mathbf{C}_{\hat{a}\bar{p}} &= -\int_{\Omega} \alpha \mathcal{H}_{\Gamma_d} \nabla \mathbf{N}^{\mathrm{T}} \mathbf{m} \mathbf{H} \mathrm{d}\Omega \\
\mathbf{C}^*_{\hat{a}\bar{p}} &= -\int_{\Gamma_d} \mathbf{N}^{\mathrm{T}} \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \mathbf{H} \mathrm{d}\Gamma \\
\mathbf{C}_{\hat{a}\hat{p}} &= -\int_{\Omega} \alpha \mathcal{H}_{\Gamma_d} \mathcal{D}_{\Gamma_d} \nabla \mathbf{N}^{\mathrm{T}} \mathbf{m} \mathbf{H} \mathrm{d}\Omega
\end{aligned}
\qquad (57) $$

The matrices $C_{\bar{a}\bar{p}}^{+}$, $C_{\hat{a}\hat{p}}^{+}$ and $K_{\bar{p}\bar{p}}^{+}$ result from differentiating the mass coupling term $\int_{\Gamma_d} \mathbf{H}^{\mathrm{T}} \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \mathbf{q}_d \mathrm{d}\Gamma$ with respect to $\bar{a}$, $\hat{a}$ and $\bar{p}$, respectively. The resulting expressions depend on the subscale model that has been chosen for the microflow, and are generally very complicated. In any case, these coupling terms and the coupling term $C_{\hat{a}\bar{p}}^{*}$ cause the tangent stiffness matrix to become unsymmetric. To restore symmetry, the non-symmetric contributions that arise from the mass and momentum couplings have been removed from the tangent stiffness matrix, and the iterations in the example calculations in the next section have been carried out with the resulting symmetric stiffness matrix.

## 7. Examples: stationary and propagating cracks

First, an open, stationary crack is considered in a linear-elastic medium and no tractions are transferred between both sides of the crack. Hence, the first subscale model for the fluid flow in the cavity is used and singularity functions that stem from linear-elastic fracture mechanics are added as enrichment functions at the crack tips [10]. The specimen is a square-shaped fractured block
<figure><img src="image_5.png" /><figcaption>Figure 4. Fine and coarse mesh composed of quadrilateral elements</figcaption></figure>

under plane-strain conditions. A normal fluid flux $q_o = 10^{-4} \text{ ms}^{-1}$ starting at $t = 0 \text{ s}$ is imposed at the bottom, 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 $10 \text{ m} \times 10 \text{ m}$ and consists of a porous material with a fluid volume 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's 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{Ns}$.

Initially, one fault is considered. This fault, at the centre of the specimen, is 2 m long, and is inclined 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 Fig. 4 has been used with quadrilateral elements equipped with quadratic shape functions.

The results obtained for an inclination angle of 30° are shown in Figs. 5–7. Because of the imposed fluid flux at the bottom, the pressure increases in the specimen and inside the fault, which subsequently opens, and fluid can flow inside, see Fig. 5. The discontinuity in the normal derivative of the pressure is illustrated in Fig. 6. 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 7 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
<figure><img src="image_6.png" /><figcaption>Figure 5. Deformed mesh for a crack angle of 30° (magnified ×1000)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. 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>

of the fault. In the cavity, the velocity of the fluid is very high compared with the velocity in the bulk because there is no resisting solid skeleton.

Next, a set of ten stationary faults is randomly generated. The length of the faults is between 1 and 3 m, while the fault angles vary from −10° to +30°. Figure 8 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 Fig. 8 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.

Finally, a cohesive crack is simulated which propagates in a plate under plane-strain conditions. Now, the second model for fluid flow inside the cavity has been used, in which progressive decohesion within the crack is assumed
<figure><img src="image_8.png" /><figcaption>Figure 7. Fluid velocity field for a crack angle of 30°</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. $L_2$ - Norm of the pressure gradient</figcaption></figure>

[11]. The plate has sides of 0.25 m and an initial notch which is located at the symmetry axis and is 0.05 m deep. A fixed vertical velocity $v = 2.35 \times 10^{-2}$ µm/s is prescribed in a opposite direction at the bottom and at the top of the plate (tensile loading). All other boundaries of the plate are assumed to be impervious. The mesh consists of 20 × 20 quadrilateral elements with bilinear
<figure><img src="image_10.png" /><figcaption>Figure 9. Pressure field in Pa. Left: The case with full coupling; Right: The case without coupling</figcaption></figure>

shape functions for the pressure and the displacement fields. The time step size is 1 s and the analysis is continued until the crack reaches the right side of the plate.

Obviously, the amount of fluid that flows into the cavity is closely related to the crack opening displacement and vanishes at the crack tip. As a consequence the water pressure decreases in the vicinity of the crack. Figure 9 gives the pressure fields for a simulation with a full coupling at the interface as derived in the previous section and for a simulation without a mass transfer coupling term (and also without pressure enrichment). In the latter case, the crack is not seen as a discontinuity in the pressure field and the fluid phase flows through the crack as it does in the bulk. Accordingly, the pressure and pressure gradient are continuous at the interface. Indeed, the results presented in the two graphs of Fig. 9 are very different. Because of the mass transfer coupling, the water is sucked into the crack and high negative pressures occur around the crack. As a consequence, the water saturation decreases in this zone and intense cavitation occurs. Because of the negative values of the water pressure, sucking tractions modify the global response of the plate. As shown in Fig. 10, the load-displacement curve obtained with the coupling term results in a higher load-carrying capacity. Indeed, the effect of the mass transfer coupling is strong and changes the fluid flow in the entire domain. Moreover, the fluid velocity is increased by an order of magnitude.

Finally, the profile of the tangential component of the Darcian velocity is shown in Fig. 11. The graphs are plotted such that the velocity is positive when it is oriented from the actual crack tip to the initial pre-notch. When the coupling is activated, the water flows from the actual tip to the initial notch with a
<figure><img src="image_11.png" /><figcaption>Figure 10. Load-displacement curves</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Water tangential Darcian velocity profile inside the crack in m/s</figcaption></figure>

value depending on the crack opening displacement whereas the orientation of the flow in the uncoupled case depends on the position on the interface.

## 8. 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 discretisation. 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 inde-
pendently 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 homogenised “constitutive” rela-
tions 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.

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*, Second Edition. 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* 1997; **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, Wien, 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] Réthoré J, de Borst R, Abellan MA. A discrete model for the dynamic propagation of shear bands in a fluid-saturated medium. *International Journal for Numerical and Analytical Methods in Geomechanics* 2007; **31**: 347–370.

[10] Réthoré J, de Borst R, Abellan MA. A two-scale approach for fluid flow in fractured porous media. *International Journal for Numerical Methods in Engineering* 2007; **71**: 780–800.

[11] Réthoré J, de Borst R, Abellan MA. A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks. *Computational Mechanics* 2008; **42**: 227–238.

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

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

[14] 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.
[15] 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.

[16] 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.

[17] 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.

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

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

[20] 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.

[21] 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.

[22] 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.

[23] de Borst R, Remmers JJC, Needleman A. Mesh-independent numerical representations of cohesive-zone models. *Engineering Fracture Mechanics* 2006; **173**: 160–177.