# A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks

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 model for fluid flow in an unsaturated porous medium with cohesive cracks. Computational Mechanics, 2008, Challenges in Computational Mechanics, 42 (2), pp.227-238. 10.1007/s00466-007-0178-6 . hal-00938323

HAL Id: hal-00938323

https://hal.science/hal-00938323

Submitted on 27 May 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 Attribution 4.0 International License
# A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks

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

**Abstract** A two-scale model is developed for fluid flow in a deforming, unsaturated and progressively fracturing porous medium. At the microscale, the flow in the cohesive crack is modelled using Darcy's relation for fluid flow in a porous medium, taking into account changes in the permeability due to the progressive damage evolution inside the cohesive zone. 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 an unsaturated 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 are independent from the underlying discretization. The resulting discrete equations are nonlinear due to the cohesive crack model and the nonlinearity 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. The calculations indicate that the evolving cohesive cracks can have a significant influence on the fluid flow and vice versa.

**Keywords** Fracture · Unsaturated porous medium · Multiscale method · Multiphase medium · Fluid flow · Cohesive crack

## 1 Introduction

Since the pioneering work of Terzaghi [24] and Biot [4] the flow of fluids in deforming porous media has received considerable attention. Recently, Lewis and Schrefler [13] have given an account of a 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 [10, 14, 23]. A general framework for accommodating multi-field problems has been presented by Jouanna and Abellan [12].

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 [6, 22]. 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 as shown by Schrefler et al. [22]. 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

<footnote>J. Réthoré · R. de Borst (✉)
Faculty of Aerospace Engineering,
Delft University of Technology, Delft, Netherlands
e-mail: R.deBorst@TUDelft.nl</footnote>

<footnote>J. Réthoré
e-mail: J.Rethore@LR.TUDelft.nl</footnote>

<footnote>M.-A. Abellan
LTDS-ENISE-UMR CNRS, 5513 Saint-Etienne, France</footnote>
influence the flow and deformations in the entire body of interest.

In this contribution, we will develop a general numerical model for flow in progressively fracturing porous media. Building on earlier work in which we have constructed a numerical model for shear-band propagation in deforming, fluid-saturated porous media [8,20], we now extend the theory to include flow inside a propagating cohesive crack. Like in Réthoré et al. [21] the flow inside the evolving crack can also 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 cohesive crack is modelled 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 an unsaturated 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 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 discretizations, these requirements can be satisfied by exploiting the partition-of-unity property of finite element shape functions [2], as has been done successfully in applications to cracking in single-phase media: Black and Belytschko [5], Moës et al. [16], Belytschko et al. [3], Wells and Sluys [25,26], Wells et al. [27,28], Remmers et al. [17], Réthoré et al. [18,19], de Borst et al. [7].

To provide a proper setting, we will first briefly recapitulate the governing equations for a deforming unsaturated 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 the subject of the subsequent section, where it will be demonstrated how the momentum and mass couplings of the micromechanical flow model to the surrounding, unsaturated porous medium can be accomplished in the weak formulation. Time integration and a consistent linearization 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 a propagating cohesive crack. 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

### 2.1 Strong form

The bulk is considered as a three-phase medium subject to the restrictions of small displacement gradients and small variations in the concentrations [12]. The problem is formulated in terms of the displacement of the solid phase $\mathbf{u}_s$ and the water and air pressures $p_w$ and $p_a$, respectively. The voids of the solid skeleton (porosity $n$) are partially filled with water and partially with air. The partition is given by the degree of saturation for each fluid phase $S_\pi$ which sum to one:

$$S_w + S_a = 1 \qquad (1)$$

The degrees of saturation are described by means of experimentally determined functions of the capillary pressure $p_c = p_a - p_w$:

$$S_w = S_w(p_c) \qquad (2)$$

In the example, the Van Genuchten relation is used, e.g. Meschke and Grasberger [15], which reads:

$$S_w(p_c) = S_{irr} + (1 - S_{irr}) \left( 1 + \left( \frac{p_c}{p_{ref}} \right)^{\frac{1}{1-l}} \right)^l \qquad (3)$$

The degree of saturation for the water is not allowed to decrease to the irreducible saturation $S_{irr}$, and the reference pressure $p_{ref}$ is used as a scaling factor for the capillarity pressure $p_c$. $l$ is the porosity index which characterizes the micro-structure of the solid porous skeleton.

The solid phase is supposed to sustain an effective stress state given by the Cauchy stress tensor $\sigma_s$. Furthermore, the assumptions are made 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 balance of linear momentum for the mixture reads, e.g. Abellan and de Borst [1]:

$$\nabla \cdot \sigma = 0 \qquad (4)$$

where $\sigma$ is the total stress and the body force has been neglected. For future use we define $p$ as the average pressure in the fluid phases. Assuming immiscibility we have:

$$p = S_w p_w + S_a p_a \qquad (5)$$

The Biot coefficient $\alpha$, cf. Lewis and Schrefler [13], is given by the following relation:

$$\alpha = 1 - K_1/K_s \qquad (6)$$
with $K_s$ the bulk modulus of the solid material and $K_t$ the overall bulk modulus of the porous medium.

Under the same assumptions as for the balance of momentum, one can write the balance of mass for each phase. This is achieved by considering a control volume moving with the solid skeleton. Over this control volume the mass balance for each phase $\pi$ reads:

$$S_{\pi} \frac{\alpha - n}{K_s} \dot{p} + \alpha S_{\pi} \nabla \cdot \mathbf{v}_s + n \dot{S}_{\pi} + \frac{n S_{\pi}}{K_{\pi}} \dot{p}_{\pi} + n \nabla \cdot (\mathbf{v}_{\pi} - \mathbf{v}_s) = 0 \quad (7)$$

with a superimposed dot denoting time differentiation, $\rho_{\pi}$ the mass density and $\mathbf{v}_{\pi}$ the absolute velocity of constituent $\pi$. The first two terms account for the changes in the porous skeleton, the third term is due to the variation in the saturation and the fourth term represents the change in mass density due to the pressure change. The last term represents the outflow from the control volume. Taking into account the expression of the average pressure (5) and Eq. (2), the time derivative of the average pressure is:

$$\dot{p} = \left( S_w + (p_a - p_w) \frac{\partial S_w}{\partial p_c} \right) \dot{p}_w + \left( S_a - (p_a - p_w) \frac{\partial S_w}{\partial p_c} \right) \dot{p}_a \quad (8)$$

The governing equations, i.e. the balance of momentum of the saturated medium, Eq. (4), and the balances of mass, Eq. (7), are complemented by the kinematic relation,

$$\epsilon_s = \nabla^s \mathbf{u}_s \quad (9)$$

with $\mathbf{u}_s$, $\epsilon_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 the solid skeleton. The effective stress increment in the solid skeleton, $d\sigma'_s$ is related to the strain increment $d\epsilon_s$ by an incrementally linear stress-strain relation for the solid skeleton,

$$d\sigma'_{s} = \bar{\mathbf{D}}^{\text{tan}} : d\epsilon_{s} \quad (10)$$

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 $\sigma'_s = \sigma_s/(1-n)$, the above relation can be replaced by

$$d\sigma_s = \mathbf{D}^{\text{tan}} : d\epsilon_s \quad (11)$$

where the notation $\mathbf{D}^{\text{tan}} = (1-n)\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 for each of the fluid phases,

$$n(\mathbf{v}_{\pi} - \mathbf{v}_{s}) = -k_{\pi} \nabla p_{\pi} \quad (12)$$

with $k_{\pi}$ the permeability coefficient of the porous medium with respect to the fluid phase $\pi$. This parameter depends on the viscosity $\mu_{\pi}$ of the phase, and on the microstructure of the solid skeleton through the intrinsic permeability $k$. A dependency on the degree of saturation can be included via the relative permeability coefficient $k_{r\pi}$, such that:

$$k_{\pi} = \frac{k}{\mu_{\pi}} k_{r\pi} (S_{\pi}) \quad (13)$$

Following Van Genuchten, the relative permeability for the water is defined as

$$k_{rw} = \sqrt{S_e} \left( 1 - (1 - S_e^T)^l \right)^2 \quad (14)$$

and that for the air phase as

$$k_{ra} = (1 - S_e)^2 \left( 1 - S_e^{\frac{2+3l}{l}} \right) \quad (15)$$

where $S_e$ is a relative saturation

$$S_e = \frac{S_w - S_{irr}}{1 - S_{irr}} \quad (16)$$

The boundary conditions (see Fig. 1)

$$\mathbf{n}_{\Gamma} \cdot \boldsymbol{\sigma} = \mathbf{t}_p, \quad \mathbf{u} = \mathbf{u}_p \quad (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(\mathbf{v}_{\pi} - \mathbf{v}_{s}) \cdot \mathbf{n}_{\Gamma} = \mathbf{q}_{\pi p} \cdot \mathbf{n}_{\Gamma}, \quad p_{\pi} = p_{\pi p} \quad (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}$

<figure><img src="image_2.png" /><figcaption>Fig. 1 Boundary conditions and internal discontinuity $\Gamma_d$ for body $\Omega$</figcaption></figure>
and $p_p$ being the prescribed outflow of pore fluids and the prescribed pressures, respectively.

## 2.2 Weak form and coupling

To arrive at the weak form of the balance equations, we multiply the momentum balance (4) and the mass balances (7) by admissible test functions for the displacements $\eta$ of the skeleton and for the pressures $\zeta_\pi$. Substitution into Equations (4) and (7), using Darcy's relation (12), 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 S_w \zeta_w \nabla \cdot \mathbf{v}_s \, d\Omega + \int_{\Omega} k_w \nabla \zeta_w \cdot \nabla p_w \, d\Omega \\ & - \int_{\Omega} \zeta_w \frac{\dot{p}_w}{Q_{ww}} \, d\Omega - \int_{\Omega} \zeta_w \frac{\dot{p}_a}{Q_{wa}} \, d\Omega \\ & + \int_{\Gamma_d} \mathbf{n}_{\Gamma_d} \cdot [[\zeta_w n(\mathbf{v}_w - \mathbf{v}_s)]] \, d\Gamma = \int_{\Gamma} \zeta_w \mathbf{n}_{\Gamma} \cdot \mathbf{q}_{wp} \, d\Gamma \end{aligned} \quad (20) $$

for the water and

$$ \begin{aligned} & - \int_{\Omega} \alpha S_a \zeta_a \nabla \cdot \mathbf{v}_s \, d\Omega + \int_{\Omega} k_a \nabla \zeta_a \cdot \nabla p_a \, d\Omega \\ & - \int_{\Omega} \zeta_a \frac{\dot{p}_w}{Q_{aw}} \, d\Omega - \int_{\Omega} \zeta_a \frac{\dot{p}_a}{Q_{aa}} \, d\Omega \\ & + \int_{\Gamma_d} \mathbf{n}_{\Gamma_d} \cdot [[\zeta_a n(\mathbf{v}_a - \mathbf{v}_s)]] \, d\Gamma = \int_{\Gamma} \zeta_a \mathbf{n}_{\Gamma} \cdot \mathbf{q}_{ap} \, d\Gamma \end{aligned} \quad (21) $$

for the air. The coefficients $Q_{ij}$ are defined as:

$$ \begin{aligned} \frac{1}{Q_{ww}} &= S_w \frac{\alpha - n}{K_s} \left( S_w + \frac{\partial S_w}{\partial p_c} p_c \right) + \frac{n S_w}{K_w} - n \frac{\partial S_w}{\partial p_c} \\ \frac{1}{Q_{wa}} &= S_w \frac{\alpha - n}{K_s} \left( 1 - S_w - \frac{\partial S_w}{\partial p_c} p_c \right) + n \frac{\partial S_w}{\partial p_c} \\ \frac{1}{Q_{aw}} &= (1 - S_w) \frac{\alpha - n}{K_s} \left( S_w + \frac{\partial S_w}{\partial p_c} p_c \right) + n \frac{\partial S_w}{\partial p_c} \\ \frac{1}{Q_{aa}} &= (1 - S_w) \frac{\alpha - n}{K_s} \left( 1 - S_w - \frac{\partial S_w}{\partial p_c} p_c \right) \\ & \qquad + \frac{n(1 - S_w)}{K_a} - n \frac{\partial S_w}{\partial p_c} \end{aligned} $$

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 stems from the cohesive tractions due to the cracking of the solid phase and the tractions applied by the fluid in the discontinuity onto the faces of the discontinuity. We assume that the fluid 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} = \mathbf{t}_d - p \mathbf{n}_{\Gamma_d} \quad \text{on } \Gamma_d \quad (22) $$

with $\mathbf{t}_d$ the cohesive tractions. Now, 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 (23) $$

Since the cohesive 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_{\Gamma_d} \mathbf{n}_{\Gamma_d} \cdot [[\zeta_w n(\mathbf{v}_w - \mathbf{v}_s)]] \, d\Gamma = \\ & \int_{\Gamma_d} \zeta_w n \mathbf{n}_{\Gamma_d} \cdot [[\mathbf{v}_w - \mathbf{v}_s]] \, d\Gamma = \int_{\Gamma_d} \zeta_w \mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_{wd} \, d\Gamma \end{aligned} \quad (24) $$

where $\mathbf{q}_{wd}$ denotes the flux of the water 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 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. Next, the influence of the tangential 'micro'-flow inside the discontinuity on that at the 'macro'-scale will be quantified.

In order to derive a relation for the mass transfer between the cavity and the bulk, we consider the fluid flow inside the cohesive zone. The permeability inside this porous and damaged zone is denoted by $k_{\pi d}$. The geometry of the cavity is shown in Fig. 2. For simplicity, and because the contribution of the gas phase can be neglected, $p_a$ is assigned a zero value over the entire domain. The problem is now specialized to an unsaturated porous medium with a passive air phase.
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Cavity geometry</figcaption></figure>

Under the assumptions listed in the preceding section, the mass balance for the water inside the cohesive zone reads:

$$ \alpha \nabla \cdot \mathbf{v}_s + n \nabla \cdot (\mathbf{v}_w - \mathbf{v}_s) + \frac{1}{Q_{ww}} \frac{\partial p_w}{\partial t} = 0 \quad (25) $$

Because the width of the cavity $2h$ is negligible compared to its length, the mass balance is enforced in an average sense over the cross section:

$$ \int_{-h}^{h} [\alpha \nabla \cdot \mathbf{v}_s + n \nabla \cdot (\mathbf{v}_w - \mathbf{v}_s) + \frac{1}{Q_{ww}} \frac{\partial p_w}{\partial t}] dy = 0 \quad (26) $$

The variation of the pressure over the height of the cavity is neglected and the last term becomes:

$$ \int_{-h}^{h} \frac{1}{Q_{ww}} \frac{\partial p_w}{\partial t} dy = 2h \frac{1}{Q_{ww}} \frac{\partial p_w}{\partial t} \quad (27) $$

For the two first terms, the divergence operator leads to different contributions. 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 \quad (28) $$

where $v_s$ and $w_s$ are the component of the solid velocity tangential and normal to the crack, respectively. The contribution that involves derivatives with respect to $y$ is directly rewritten as:

$$ \int_{-h}^{h} \alpha \frac{\partial w_s}{\partial y} dy = \alpha[w_s(h) - w_s(-h)] = \alpha[\llbracket w_s \rrbracket] \quad (29) $$

The second term is elaborated as follows:

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

Assuming that $v_s$ varies linearly with $y$, the integral can be solved analytically:

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

Because $\alpha$ depends on the capillarity pressure only, it is supposed to be constant over a cross section. Rearranging Eq. (31) and defining $\langle v_s \rangle = \frac{v_s(h)+v_s(-h)}{2}$:

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

Applying the same operations to the second term of Eq. (26), the following expression can be derived:

$$ \int_{-h}^{h} n \nabla \cdot (\mathbf{v}_w - \mathbf{v}_s) dy = \llbracket n(w_w - w_s) \rrbracket \\ + \frac{\partial}{\partial x} \left( \int_{-h}^{h} n(v_w - v_s) dy \right) \\ - \left( n(v_w - v_s)(h) \frac{\partial h}{\partial x} - n(v_w - v_s)(-h) \frac{\partial (-h)}{\partial x} \right) \quad (33) $$

The term $\llbracket n(w_w - w_s) \rrbracket$ can be identified as the coupling term $n \Gamma_d \cdot q_d$ of the weak form of the mass balance. Introducing Darcy's relation projected onto the axis tangential to the crack, $n(v_w - v_s) = -k \frac{\partial p_w}{\partial x}$, with $k$ a generic permeability, the following relation is obtained:

$$ \int_{-h}^{h} n \nabla \cdot (\mathbf{v}_w - \mathbf{v}_s) dy = \llbracket n(w_w - w_s) \rrbracket - \frac{\partial}{\partial x} \left( \int_{-h}^{h} k \frac{\partial p_w}{\partial x} dy \right) \\ + \left( k \frac{\partial p_w}{\partial x}(h) \frac{\partial h}{\partial x} - k \frac{\partial p_w}{\partial x}(-h) \frac{\partial (-h)}{\partial x} \right) \quad (34) $$
At the faces of the cavity the permeability equals that of the bulk $k_w$, and for continuity reasons $\frac{\partial p_w}{\partial x}(h) = \frac{\partial p_w}{\partial x}(-h)$. Inside the crack, $k$ is assumed not to depend on $y$, but it is supposed to be affected by the decohesion of the solid skeleton. Consequently, the following expression is assigned to the permeability inside the cavity:

$$k = k_{wd}(S_w, \llbracket \mathbf{u}_s \rrbracket) \qquad (35)$$

which will be specialized in the next subsection. Finally, introducing this expression in Eq. (34) and inserting the result obtained in Eq. (26), the coupling term is written as:

$$\mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_{wd} = 2h \frac{1}{Q_{ww}} \frac{\partial p_w}{\partial t} + \alpha \llbracket w_s \rrbracket \\ +2\alpha \langle \frac{\partial v_s}{\partial x} \rangle h - \frac{\partial}{\partial x} \left( 2k_{wd} \frac{\partial p_w}{\partial x} h \right) \\ +2k_w \frac{\partial p_w}{\partial x} \frac{\partial h}{\partial x} \qquad (36)$$

## 2.3 Cohesive interface relation

The decohesion of the solid skeleton inside the interface is modelled using a discrete relation between the interface tractions $\mathbf{t}_d$ and the relative displacements $\llbracket \mathbf{u}_s \rrbracket$:

$$\mathbf{t}_d = \mathbf{t}_d(\llbracket \mathbf{u}_s \rrbracket, \kappa) \qquad (37)$$

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

$$\delta \mathbf{t}_d = \mathbf{T} \delta \llbracket \mathbf{u}_s \rrbracket \qquad (38)$$

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

$$\mathbf{T} = \frac{\partial \mathbf{t}_d}{\partial \llbracket \mathbf{u}_s \rrbracket} + \frac{\partial \mathbf{t}_d}{\partial \kappa} \frac{\partial \kappa}{\partial \llbracket \mathbf{u}_s \rrbracket} \qquad (39)$$

Whichever formalism is used, plasticity or damage, a key element in a cohesive law is the presence of a work of separation or fracture energy, $G_c$, which governs nucleation of voids and enters the interface constitutive relation (37) in addition to the tensile strength $f_t$. It is defined as the work needed to create a unit area of fully developed crack:

$$G_c = \int_0^\infty \mathbf{t}_d(\llbracket \mathbf{u}_s \rrbracket, \kappa) \, d\llbracket \mathbf{u}_s \rrbracket \qquad (40)$$

and thus equals the area under the decohesion curve. The initiation criterion is written in term of maximum traction at the tip

$$t_{max} = \max_{\theta \in [-\pi, \pi]} t_{eq} = t_{eq}(\theta_{max}) \ge t_c \qquad (41)$$

where $t_c$ is the cohesive strength and $\mathbf{t}$ gives the angle $\theta$ with the local tangent vector to the crack. The equivalent traction

$t_{eq}$ [9] is computed from the tractions of an averaged stress tensor at the crack tip:

$$t_{eq} = \sqrt{t_n^2 + t_r^2 / \beta} \qquad (42)$$

where

$$t_n = \mathbf{n} \cdot \boldsymbol{\sigma}_{tip} \cdot \mathbf{n}, \quad t_s = \mathbf{t} \cdot \boldsymbol{\sigma}_{tip} \cdot \mathbf{n} \qquad (43)$$

Typically, $\beta = 1.5$.

The stresses vary strongly in the vicinity of the tip and an accurate computation is less straightforward. Following Wells and Sluys [25] and Jirasek [11] we use a smoothing of the stresses around the tip and compute the stress at the tip by the following nonlocal-like procedure:

$$\boldsymbol{\sigma}_{tip} = \frac{\int_{\Omega} w \boldsymbol{\sigma} \, d\Omega}{\int_{\Omega} w \, d\Omega} \qquad (44)$$

where $w$ is a Gaussian weight function:

$$w = e^{-r^2 / 2l^2}$$

with $r$ the distance to the tip, and $l$ a characteristic length which defines the size of region of influence of the stress. Because of the variation of the stress, a small value of $l$ is desired, preferably in the same order of magnitude as the characteristic element length. This is accomplished in the following manner. By virtue of the linear behaviour of the solid phase in the bulk, a separate, independent integration domain can be defined, which follows the tip during propagation. This domain contains integration cells smaller than those of the mesh used in the discretization—typically their length is in the order of 15–20% of the element size. Moreover, a higher-order Gaussian quadrature is used over this domain, which results in a very accurate determination of the tip stress.

For the bulk material, it has been assumed that the current value of the porosity equals its initial value. Thus, the intrinsic permeability of the skeleton is not affected by the strains. This hypothesis is considered to be less applicable to the cohesive zone. At this point, discrete relation is introduced between the porosity inside the interface, $n_d$, and the displacement jump $\llbracket \mathbf{u}_s \rrbracket$:

$$n_d = n_d(\llbracket \mathbf{u}_s \rrbracket) \qquad (45)$$

The change in porosity affects the permeability in the process zone via a coefficient $k_n$, as $k_{wd}$ is defined by:

$$k_{wd} = \frac{k k_n k_{rw}}{\mu_w} \qquad (46)$$

where $k_n$ is given by the following heuristic relation [15]:

$$k_n(n_d) = 10^{\delta_n} \quad \text{with} \quad \delta_n = \frac{6(n_d - n)}{0.3 - 0.4n} \qquad (47)$$
# 3 Numerical elaboration

## 3.1 Discretization

The model presented in the previous section puts some restrictions 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 [2].

In the spirit of previous work on cracks in single-phase media [3, 17] the interpolation of each component of the displacement field of the solid phase is enriched with discontinuous 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 \quad (48)$$

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 cut by the discontinuity, Fig. 3. They hold additional degrees of freedom $\hat{u}_i$ corresponding to the discontinuous function

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

<figure><img src="image_4.png" /><figcaption>Fig. 3 Enrichment strategy: circles denote $\mathcal{N}_{\text{cut}}$</figcaption></figure>

Symbolically, Eq. (48) can be written as

$$\mathbf{u}_s = \mathbf{N}\mathbf{U} \quad (50)$$

where the matrix $\mathbf{N}$ contains the standard interpolation polynomials $N_i$ as well as the discontinuous function $\mathcal{H}_{\Gamma_d}$, and the array $\mathbf{U}$ contains the displacement degrees of freedom $\bar{u}_i$ and $\hat{u}_i$. Furthermore, we define the matrix $\llbracket \mathbf{N} \rrbracket$ which contains the jumps in the interpolation polynomials, and the matrix $\mathbf{B}$, which contains the bounded part of spatial derivatives of $\mathbf{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} \quad (51)$$

All nodes whose support is cut by the discontinuity, $\mathcal{N}_{\text{cut}}$, hold additional pressure degrees of freedom:

$$p_w = \sum_{i \in \mathcal{N}} H_i \bar{p}_i + \sum_{i \in \mathcal{N}_{\text{cut}}} H_i \mathcal{D}_{\Gamma_d} \hat{p}_i \quad (52)$$

$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 Eq. (52) symbolically as:

$$p_w = \mathbf{H}\mathbf{P} \quad (53)$$

where $\mathbf{H}$ contains the shape functions $H_i$ as well as the distance function $\mathcal{D}_{\Gamma_d}$, and the array $\mathbf{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 Eq. (36). 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 the order of $H_i$ for consistency in the discrete balance of momentum equation. In practice, we will use quadrilateral elements with bilinear shape functions for the displacement as well as for the pressure discretization.

## 3.2 Discrete equations and resolution

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

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

$$ -\int_{\Omega} \alpha \mathbf{H}^{\mathrm{T}} \mathbf{m}^{\mathrm{T}} \nabla \hat{\mathbf{u}}_{s} \mathrm{d}\Omega + \int_{\Omega} k_{w} \nabla \mathbf{H}^{\mathrm{T}} \nabla p_{w} \mathrm{d}\Omega $$

$$ -\int_{\Omega} Q_{ww}^{-1} \mathbf{H}^{\mathrm{T}} \dot{\mathbf{p}}_{w} \mathrm{d}\Omega + \mathbf{Q}_{\mathrm{inter}} = \mathbf{Q}_{\mathrm{ext}} \quad (55) $$

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}_{p} \mathrm{d}\Gamma \quad (56) $$

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

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

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 (55) leads to:

$$ -\int_{\Omega} \alpha \mathbf{H}^{\mathrm{T}} \mathbf{m}^{\mathrm{T}} \Delta \nabla \mathbf{u}_{s} \mathrm{d}\Omega + \Delta t \int_{\Omega} k_{w} \nabla \mathbf{H}^{\mathrm{T}} \nabla p_{w} \mathrm{d}\Omega \\ -\int_{\Omega} Q_{ww}^{-1} \mathbf{H}^{\mathrm{T}} \Delta p_{w} \mathrm{d}\Omega + \Delta t \mathbf{Q}_{\mathrm{inter}} = \Delta t \mathbf{Q}_{\mathrm{ext}} \quad (59) $$

The interfacial force vector $\mathbf{F}_{\mathrm{inter}}$ is derived from Eq. (23) by integrating $\llbracket \boldsymbol{\eta} \rrbracket \cdot (\mathbf{t}_d - p_w \mathbf{n}_{\Gamma_d})$ along $\Gamma_d$:

$$ \mathbf{F}_{\mathrm{inter}} = \int_{\Gamma_d} \llbracket \mathbf{N} \rrbracket^{\mathrm{T}} \mathbf{t}_d \mathrm{d}\Gamma - \left( \int_{\Gamma_d} \llbracket \mathbf{N} \rrbracket^{\mathrm{T}} \mathbf{n}_{\Gamma_d} \mathbf{H} \mathrm{d}\Gamma \right) \mathbf{P} \quad (60) $$

In a similar fashion, Eq. (24) is integrated along $\Gamma_d$ to give:

$$ \mathbf{Q}_{\mathrm{inter}} = \int_{\Gamma_d} \mathbf{H}^{\mathrm{T}} \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \mathbf{q}_{wd} \mathrm{d}\Gamma \quad (61) $$

The interfacial flux vector $\mathbf{Q}_{\mathrm{inter}}$ is then elaborated by incorporating the discrete form of Eq. (36), yielding:

$$ \begin{aligned} \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \mathbf{q}_{wd} = \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \llbracket \mathbf{N} \rrbracket \mathbf{U} \frac{1}{Q_{ww}} \mathbf{H} \dot{\mathbf{P}} + \alpha \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \llbracket \mathbf{N} \rrbracket \dot{\mathbf{P}} \\ + \alpha \left( \mathbf{t}_{\Gamma_d}^{\mathrm{T}} \langle \mathbf{N} \rangle \dot{\mathbf{U}} \right) \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \llbracket \mathbf{N} \rrbracket \mathbf{U} \\ - \mathbf{t}_{\Gamma_d}^{\mathrm{T}} \nabla \left( k_{wd} \left( \mathbf{t}_{\Gamma_d}^{\mathrm{T}} \nabla \mathbf{H} \mathbf{P} \right) \left( \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \llbracket \mathbf{N} \rrbracket \mathbf{U} \right) \right) \\ + k_w \left( \mathbf{t}_{\Gamma_d}^{\mathrm{T}} \nabla \mathbf{H} \mathbf{P} \right) \left( \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \llbracket \nabla \mathbf{N} \rrbracket \mathbf{t}_{\Gamma_d} \mathbf{U} \right) \end{aligned} \quad (62) $$

Because of the cohesive-zone model in the interface, and since $\mathbf{Q}_{\mathrm{inter}}$ and $\mathbf{F}_{\mathrm{inter}}$ are nonlinear, 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}^{\mathrm{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(\frac{\mathbf{U}}{\mathbf{P}}\right)^i + \left( \frac{\mathbf{F}_{\mathrm{inter}}}{\Delta t \mathbf{Q}_{\mathrm{inter}}} \right)^i - \left( \frac{\mathbf{F}_{\mathrm{ext}}}{\Delta t \mathbf{Q}_{\mathrm{ext}}} \right) \quad (63) $$

where the stiffness matrices are defined as:

$$ \mathbf{K}_{uu} = \int_{\Omega} \mathbf{B}^{\mathrm{T}} \mathbf{D}^{\mathrm{tan}} \mathbf{B} \mathrm{d}\Omega \quad (64) $$

$$ \mathbf{K}_{up} = - \int_{\Omega} \alpha \mathbf{B}^{\mathrm{T}} \mathbf{m} \mathbf{H} \mathrm{d}\Omega \quad (65) $$

$$ \mathbf{K}_{pp}^{(1)} = - \int_{\Omega} Q_{ww}^{-1} \mathbf{H}^{\mathrm{T}} \mathbf{H} \mathrm{d}\Omega \quad (66) $$

$$ \mathbf{K}_{pp}^{(2)} = - \int_{\Omega} k_w \nabla \mathbf{H}^{\mathrm{T}} \nabla \mathbf{H} \mathrm{d}\Omega \quad (67) $$

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} + \frac{\partial \mathbf{F}_{\mathrm{inter}}}{\partial \mathbf{U}} & \mathbf{K}_{up} + \frac{\partial \mathbf{F}_{\mathrm{inter}}}{\partial \mathbf{P}} \\ \mathbf{K}_{up}^{\mathrm{T}} + \Delta t \frac{\partial \mathbf{Q}_{\mathrm{inter}}}{\partial \mathbf{U}} & \mathbf{K}_{pp}^{(1)} + \Delta t \mathbf{K}_{pp}^{(2)} + \Delta t \frac{\partial \mathbf{Q}_{\mathrm{inter}}}{\partial \mathbf{P}} \end{bmatrix} \quad (68) $$

where all quantities are evaluated at iteration $i$.

The coupling terms, $\mathbf{F}_{\mathrm{inter}}$ and $\mathbf{Q}_{\mathrm{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} + \frac{\partial \mathbf{F}_{\mathrm{inter}}}{\partial \mathbf{U}} & \mathbf{K}_{up} \\ \mathbf{K}_{up}^{\mathrm{T}} & \mathbf{K}_{pp}^{(1)} + \Delta t \mathbf{K}_{pp}^{(2)} \end{bmatrix} \quad (69) $$

where

$$ \frac{\partial \mathbf{F}_{\mathrm{inter}}}{\partial \mathbf{U}} = \int_{\Gamma_d} \mathbf{N}^{\mathrm{T}} \mathbf{T} \mathbf{N} \mathrm{d}\Gamma \quad (70) $$

This restoration of symmetry probably decreases the convergence rate as this leads to a modified Newton-Raphson algorithm. Nevertheless, the symmetric format of the iterative matrix allows for more flexibility in the implementation as well as a better conditioning of the matrix.
## 4 Example calculation

To illustrate the effects of the fluid flow inside the crack at the ‘macro’-scale as well as at the ‘micro’-scale, we simulate the rupture of a plate under plane-strain conditions. As shown in Fig. 4, a pre-notch is located at the symmetry axis of the plate. The plate has sides of 0.25 m and the pre-notch 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 boundaries of the plate are assumed to be impervious. The analysis is continued until the crack reaches the right side of the plate.

The material properties are summarized in Table 1. The mesh consists of 20 × 20 quadrilateral elements with bilinear shape functions for the pressure and the displacement fields. This degree is high enough to obtain non-zero second derivatives of the pressure field and the convergence of the discretization scheme has been checked in Réthoré et al. [21]. The time step size is 1 s.

Figure 5 gives the amount of fluid attracted into the cavity of the crack. Obviously, this amount is closely related the crack opening displacement and vanishes at the crack tip. As a consequence the water pressure decreases in the vicinity of the crack.

Figure 6 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 recognized as a discontinuity in the pressure field and the fluid phase flows through the crack

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Geometry and loading for the pre-notched plate</figcaption></figure>

<table><caption>**Table 1** Material properties</caption><thead><tr><th>Name</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Young's modulus</td><td><i>E</i></td><td>25.85 GPa</td></tr><tr><td>Poisson's ratio</td><td><i>ν</i></td><td>0.18</td></tr><tr><td>Tensile strength</td><td><i>t</i><sub>c</sub></td><td>2.7 MPa</td></tr><tr><td>Cohesive fracture energy</td><td><i>G</i><sub>c</sub></td><td>95.0 J/m<sup>2</sup></td></tr><tr><td>Initial porosity</td><td><i>n</i></td><td>0.2</td></tr><tr><td>Intrinsic permeability</td><td><i>k</i></td><td>2.78e-21 m<sup>2</sup></td></tr><tr><td>Biot's coefficient</td><td><i>α</i></td><td>1.0</td></tr><tr><td>Solid bulk modulus</td><td><i>K</i><sub>s</sub></td><td>13.46 GPa</td></tr><tr><td>Water density</td><td><i>ρ</i><sub>w</sub></td><td>1000.0 kg/m<sup>3</sup></td></tr><tr><td>Water bulk modulus</td><td><i>K</i><sub>w</sub></td><td>0.2 GPa</td></tr><tr><td>Water dynamic viscosity</td><td><i>μ</i><sub>w</sub></td><td>5.0e-4 N/m<sup>2</sup>s</td></tr><tr><td>Irreducible saturation</td><td><i>S</i><sub>irr</sub></td><td>0.0</td></tr><tr><td>Reference pressure</td><td><i>p</i><sub>ref</sub></td><td>18.6 MPa</td></tr><tr><td>Porous index</td><td><i>l</i></td><td>0.4396</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Profile of the coupling term $n_{\Gamma_d} \cdot q_{wd}$ along the crack [m$^2$/s]</figcaption></figure>

as it does in the bulk. Accordingly, the pressure and pressure gradient are continuous at the interface. The results presented in the two graphs of Fig. 6 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, Fig. 7 (left).

Because of the negative values of the water pressure, sucking tractions modify the global response of the plate. As shown in Fig. 8, the load-displacement curve obtained with the coupling term result in a higher load-carrying capacity.

This effect can also be observed on Fig. 9 which shows the norm of the pressure gradient distribution. The coupling induces high gradients in the elements surrounding the crack.
<figure><img src="image_7.png" /><figcaption>**Fig. 7** Water saturation. Left
The case with full coupling;
Right The case without coupling</figcaption></figure>

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

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Load-displacement curves</figcaption></figure>

The water velocity fields shown in Fig. 10 also illustrate the
difference between both simulations. 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 the tangential component of the Dar-
cian velocity is plotted 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 flow from the actual tip to the initial
pre-notch with a value depending on the crack opening dis-
placement whereas the orientation of the flow in the uncou-
pled case depends on the position on the interface (the peak
close to the actual crack is due to a numerical effect). Again,
the velocity is increased by the mass transfer coupling.

**5 Concluding remarks**

A methodology has been proposed to insert discontinuities
such as cracks, faults, or shear bands, in an unsaturated porous
medium. The discontinuities can be located arbitrarily, not
related to the underlying discretization. For the fluid flow
<figure><img src="image_10.png" /><figcaption>**Fig. 10** Water velocity in m/s. *Above* The case with full coupling; *Below* The case without coupling</figcaption></figure>

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

in the progressively fracturing 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
<figure><img src="image_12.png" /><figcaption>**Fig. 11** Water tangential Darcian velocity profile inside the crack in m/s</figcaption></figure>

property of the finite element polynomial shape functions. The examples confirm the efficiency of the method and the strong effect of the coupling between the flow induced by damage inside the crack and the flow in the bulk.

References

1. Abellan MA, de Borst R (2006) Wave propagation and localisation in a softening two-phase medium. Comput Methods Appl Mech Eng 195(37-40):5011-5019

2. Babuska I, Melenk JM (1997) The partition of unity method. Int J Numer Methods Eng 40:727-758

3. Belytschko T, Moës N, Usui S, Parimi C (2001) Arbitrary discontinuities in finite elements. Int J Numer Methods Eng 50(4):993-1013

4. Biot MA (1965) Mechanics of incremental deformations. Wiley, Chichester

5. Black T, Belytschko T (1999) Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng 45:601-620

6. Boone T, Ingraffea AR (1990) A numerical procedure for simulation of hydraulically-driven fracture propagation in poroelastic media. Int J Numer Anal Methods Geomech 14(2):27-47

7. de Borst R, Remmers JJC, Needleman A (2006) Mesh-independent numerical representations of cohesive-zone models. Eng Fracture Mech 173(2):160-177

8. de Borst R, Réthoré J, Abellan MA (2006) A numerical approach for arbitrary cracks in a fluid-saturated medium. Arch Appl Mech 75:595-606

9. Camacho GT, Ortiz M (1996) Computational modeling of impact damage in brittle materials. Int J Solids Struct 33:1267-1282

10. Huyghe HM, Janssen JD (1997) Quadriphasic mechanics of swelling incompressible media. Int J Eng Sci 35:793-802

11. Jirasek M (1998) Computational modelling of concrete structures, Balkema, Rotterdam, chap Embedded crack models for concrete fracture, pp 291-300

12. Jouanna P, Abellan MA (1995) Modern issues in non-saturated soils, Springer, Heidelberg, chap Generalized approach to heterogeneous media

13. Lewis RW, Schrefler BA (1998) The finite element method in the static and dynamic deformation and consolidation of porous media, 2nd edn. Wiley, Chichester

14. van Loon R, Huyghe JM, Wijlaars MW, Baaijens FPT (2003) 3d fe implementation of an incompressible quadriphasic mixture model. Int J Numer Methods Eng 57:1243-1258

15. Meschke G, Grasberger S (2003) Numerical modeling of coupled hygromechanical degradation of cementitious materials. J Eng Mech 129(4):383-392

16. Moës N, Dolbow J, Belytschko T (1999) A finite element method for crack growth without remeshing. Int J Numer Methods Eng 46(1):133-150

17. Remmers JJC, de Borst R, Needleman A (2003) A cohesive segments method for the simulation of crack growth. Comput Mech 31:69-77

18. Réthoré J, Gravouil A, Combescure A (2005) A combined space time extended finite element method. Int J Numer Methods Eng 64:260-284

19. Réthoré J, Gravouil A, Combescure A (2005) An energy conserving scheme for dynamic crack growth with the extended finite element method. Int J Numer Methods Eng 63:631-659

20. Réthoré J, de Borst R, Abellan MA (2007a) A discrete model for the dynamic propagation of shear bands in a fluid-saturated medium. Int J Numer Anal Methods Geomech 31:347-370

21. Réthoré J, de Borst R, Abellan MA (2007b) A two-scale approach for fluid flow in fractured porous media. Int J Numer Methods Eng. doi: 10.1002/nme.1962

22. Schrefler BA, Secchi S, Simoni L (2006) On adaptive refinement techniques in multi-field problems including cohesive fracture. Comput Methods Appl Mech Eng 195(4-6):444-461

23. Snijders H, Huyghe JM, Janssen JD (1995) Triphasic finite element model for swelling porous media. Int J Numer Methods Fluids 20:1039-1046

24. Terzaghi K (1943) Theoretical soil mechanics. Wiley, New York

25. Wells GN, Sluys LJ (2001) Discontinuous analysis of softening solids under impact loading. Int J Numer Anal Methods Geomech 25:691-709

26. Wells GN, Sluys LJ (2001) A new method for modeling cohesive cracks using finite elements. Int J Numer Methods Eng 50:2667-2682

27. Wells GN, de Borst R, Sluys LJ (2002) A consistent geometrically non-linear approach for delamination. Int J Numer Methods Eng 54:1333-1355

28. Wells GN, Sluys LJ, de Borst R (2002) Simulating the propagation of displacement discontinuities in a regularized strain-softening medium. Int J Numer Methods Eng 53:1235-1256