# A discrete model for the propagation of discontinuities in a fluid-saturated medium

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 discrete model for the propagation of discontinuities in a fluid-saturated medium. IUTAM Symposium on Discretization Methods for Evolving Discontinuities, Sep 2006, Lyon, France. 10.1007/978-1-4020-6530-9_19. hal-00506330

HAL Id: hal-00506330

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

Submitted on 26 Mar 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 discrete model for the propagation of discontinuities in a fluid-saturated medium

Julien Réthoré<fnref n="1" />, René de Borst¹², and Marie-Angèle Abellan<fnref n="3" />

**Summary.** The first part of this manuscript discusses a finite element method that captures arbitrary discontinuities in a two-phase medium by exploiting the partition-of-unity property of finite element shape functions. The fluid flow away from the discontinuity is modelled in a standard fashion using Darcy's relation, and at the discontinuity a discrete analogy of Darcy's relation is used. Subsequently, dynamic shear banding is studied numerically for a biaxial, plane-strain specimen. A Tresca-like as well as a Coulomb criterion are used as nucleation criterion. Decohesion is controlled by a mode-II fracture energy, while for the Coulomb criterion, frictional forces are transmitted across the interface in addition to the cohesive shear tractions. The effect of the different interface relations on the onset of cavitation is studied.

**Key words:** shear band, dynamic fracture, two-phase medium, partition-of-unity method

## 1 Introduction

Broadly speaking, two approaches exist for the numerical analysis of the nucleation and propagation of discontinuities in solids, such as cracks, shear bands and faults. Within the classical theory of continuum mechanics, the approach in which discontinuities are distributed over a finite volume, so that relative displacements across the faces of a discontinuity are transformed into strains, is perhaps the most natural. Also in an engineering sense, it offers advantages, since there is no need to keep track of all individual microcracks that arise in a solid. Indeed, for computations of large structures, any attempt to model each individual crack would exceed even the currently available computing power. However, from a theoretical point of view, the modelling of discontinuities in a distributed or smeared sense has a limitation, namely that at a certain level of accumulated damage, the set of governing equations locally changes

<footnote n="1">¹ Faculty of Aerospace Engineering, Delft University of Technology, Delft, Netherlands r.deborst@tudelft.nl</footnote>

<footnote n="2">² LaMCoS - UMR CNRS 5514, INSA de Lyon, Villeurbanne, France</footnote>

<footnote n="3">³ LTDS-ENISE - UMR CNRS 5513, Saint-Etienne, France</footnote>
character, from elliptic to hyperbolic for quasi-static loadings, and from hyperbolic to elliptic for dynamic loadings. Unless a regularisation is applied, the resulting initial/boundary value problem becomes ill-posed, resulting in numerical solutions that depend severely on the discretisation [1].

Intuitively the most appealing approach is to model discontinuities in a discrete manner, thus reflecting the change in topology that actually takes place in the solid when a discontinuity propagates. Recently, such finite element methods have been constructed that exploit the partition-of-unity property of finite element shape functions [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. While preserving the original discretisation, the addition of extra degrees-of-freedom to nodes whose support is crossed by a discontinuity allows to construct two continuous displacement fields that are separated by a Heaviside function at the discontinuity. As a consequence, discontinuities can propagate, not biased by the original discretisation.

Many problems in geomechanics involve the coupling of the set of equations that describe the stress evolution and those which describe diffusion-type processes, e.g. water or ion transport. Indeed, hydro-mechanical interactions have been recognised to play a crucial role in geotechnical, petroleum and mining engineering since the pioneering works by Terzaghi [13] and Biot [14]. It is the purpose of this manuscript to formulate a numerical model that is capable of describing dynamic shear band propagation in a porous medium, with a solid skeleton and an interstitial fluid as the constituent phases, in a discrete, mesh-independent manner. The model exploits the partition-of-unity property of finite element shape functions, and can therefore be considered to be an extension to earlier works on fracture for single-phase media. On the other hand, the present methodology can be extended in a fairly straightforward manner to introduce discontinuities, including cracks, in initial value problems where several diffusion-type problems play a role.

The manuscript starts with a concise derivation of the balance equations for a fluid-saturated porous medium. Subsequently, the general methodology and the assumptions regarding the introduction of a discontinuity in a finite element model are discussed. The approach is specialised to a medium where the (discrete) failure mode is caused by exhaustion of the shear stress capacity on a critical plane. Next, studies are carried out for a plane-strain, biaxial specimen. Two different nucleation criteria for shear band propagation are investigated, and the role of localisation on cavitation in a fluid-saturated porous medium is highlighted.

## 2 Balance equations

We consider a two-phase medium subject to the restriction of small displacement gradients and small variations in the concentrations [15]. 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} = \mathbf{0} \quad (2) $$

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

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

Adding both momentum balances, and taking into account eq. (2), one obtains the momentum balance for the mixture:

$$ \nabla \cdot \boldsymbol{\sigma} - \rho_s \frac{\partial \mathbf{v}_s}{\partial t} - \rho_f \frac{\partial \mathbf{v}_f}{\partial t} = \mathbf{0} \quad (4) $$

where the stress is, as usual, composed of a solid and a fluid part,

$$ \boldsymbol{\sigma} = \boldsymbol{\sigma}_s + \boldsymbol{\sigma}_f \quad (5) $$

For relatively slow dynamic loadings, the assumption is often made that the accelerations of the solid and of the fluid are equal: $\frac{\partial \mathbf{v}_s}{\partial t} \approx \frac{\partial \mathbf{v}_f}{\partial t}$. With the mass density of the mixture, $\rho = \rho_s + \rho_f$, the balance of momentum (4) reduces to:

$$ \nabla \cdot \boldsymbol{\sigma} - \rho \frac{\partial \mathbf{v}_s}{\partial t} = \mathbf{0} \quad (6) $$

Numerical analyses are usually conducted with the latter equation as balance of momentum, cf [16], but the accuracy of this assumption is seldom quantified.

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 \quad (7) $$

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

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

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

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

with $K_s$ the bulk modulus of the solid material and $K_t$ the overall bulk modulus of the porous medium. Using the definition of the Biot coefficient, $1 - \alpha = K_t/K_s$ [16], 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 (12) $$

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

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

with the overall compressibility, or Biot modulus

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

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

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

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

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

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 \qquad (17)$$

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 which 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 \qquad (18)$$

close the initial value problem.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Body composed of continuous displacement fields at each side of the discontinuity $\Gamma_d$</figcaption></figure>

## 3 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 [3, 4], exploiting the partition-of-unity property of finite element shape functions [2]. 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{\bar{a}}_j + \sum_{k=1}^{m} \psi_k \dot{\tilde{a}}_{jk} \right) \qquad (19)$$

with $\bar{a}_j$ the ‘regular’ nodal degrees-of-freedom for the displacements, $\psi_k$ the enhanced basis terms, and $\tilde{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 Figure 1). The velocity field $\mathbf{v}$ can be written as the sum of two continuous velocity fields $\bar{\mathbf{v}}$ and $\tilde{\mathbf{v}}$:

$$\mathbf{v} = \bar{\mathbf{v}} + \mathcal{H}_{\Gamma_d} \tilde{\mathbf{v}} \qquad (20)$$

where $\mathcal{H}_{\Gamma_d}$ is the Heaviside step function centred at the discontinuity. The decomposition in eq. (20) has a structure similar to the interpolation in eq. (19),
e.g. [17]. 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 and shear bands is preserved. With the standard small-strain assumption that the strain-rate field of the solid, $\epsilon_s$, is derived from the symmetric part of the gradient of the velocity field, we obtain:

$$ \dot{\epsilon}_s = \nabla^s \bar{\nu}_s + \mathcal{H}_{\Gamma_d} \nabla^s \tilde{\nu}_s + \delta_{\Gamma_d} (\tilde{\nu}_s \otimes \mathbf{n}_{\Gamma_d})^s \quad (21) $$

with the superscript s denoting the symmetric part of the gradient operator.

With respect to the pore fluid, we consider the case that a diaphragm with a permeability $k_d$ is placed at the discontinuity in the displacement. As a consequence, the fluid pressure can be discontinuous across $\Gamma_d$ and, similar to eq. (20), we have:

$$ p = \bar{p} + \mathcal{H}_{\Gamma_d} \tilde{p} \quad (22) $$

It is noted that this assumption is different from that of Armero and Callari [18], who adopt a smooth pressure field (and therefore $p = \bar{p}$) and is also different from that of Larsson and Larsson [19], who assume that a regularised Dirac distribution is added to the continuous pressure field at the location of the discontinuity in the displacement field. For the fluid flow, gradients of the pressure need to be computed. Differentiating eq. (22), we obtain:

$$ \nabla p = \nabla \bar{p} + \mathcal{H}_{\Gamma_d} \nabla \tilde{p} + \delta_{\Gamma_d} \tilde{p} \mathbf{n}_{\Gamma_d} \quad (23) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Relation between relative sliding at the discontinuity and shear tractions</figcaption></figure>

# 4 Constitutive equations

## 4.1 Models for the bulk

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\boldsymbol{\sigma}'_s = \bar{\mathbf{D}}^{\text{tan}} : d\boldsymbol{\epsilon}_s \qquad (24)$$

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

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

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

For the flow of the pore fluid, Darcy's relation for isotropic media is assumed to hold,

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

with $k_f$ the permeability coefficient of the porous medium. For loading situations in which high strain rates play a significant role, Darcy's relation can be extended with a so-called dynamic seepage term [16, 20], which results in:

$$n_f(\mathbf{v}_f - \mathbf{v}_s) = -k_f \left( \nabla p + \rho_f \frac{\partial \mathbf{v}_f}{\partial t} \right) \qquad (27)$$

In line with the earlier assumption to neglect the gravity acceleration, this term has also been omitted here. In practical situations, following the assumption $\frac{\partial \mathbf{v}_s}{\partial t} \approx \frac{\partial \mathbf{v}_f}{\partial t}$ for relatively slow dynamic loadings, eq. (27) is often approximated by

$$n_f(\mathbf{v}_f - \mathbf{v}_s) = -k_f \left( \nabla p + \rho_f \frac{\partial \mathbf{v}_s}{\partial t} \right) \qquad (28)$$

## 4.2 Interface behaviour

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) \qquad (29)$$

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{\boldsymbol{\delta}} \qquad (30)$$

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

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

A first possibility that has been used in the example calculations for shear band initiation is the use of a maximum shear stress criterion in the spirit
of Tresca. With the resolved shear stress $\tau$, a shear-band like discontinuity is created when the criterion equals the critical value $\tau_c$: $\tau = \tau_c$. The orientation of the interface is such that it maximises the shear stress. In this orientation, the shear stress $\tau = \frac{|\sigma_1 - \sigma_2|}{2}$, $\sigma_1$, $\sigma_2$ being the principal stresses. A maximum shear stress nucleation criterion is primarily applicable when compressive stress states around the discontinuity prevail, such as in rocks and soils. Then, the failure mode will only involve sliding at the discontinuity, but no crack opening. For this reason, in the example calculations only degrees-of-freedom that describe this sliding mode have been added to the finite element model, which is different from earlier shear-band simulations (for single-phase media) [5, 6, 12] that have exploited the partition-of-unity property of finite element shape functions, but is similar to [9]. Dilatancy in the shear band can be incorporated when, in addition to the tangential degrees-of-freedom, during propagation extra degrees-of-freedom are activated which are normal to the shear band.

A key element is the presence of a mode-II fracture energy, $G_c^{II}$, which governs the shear band evolution and enters the interface constitutive relation (29) in addition to the shear strength $\tau_c$. It is defined as the work needed to create a unit area of fully developed shear band, e.g. [21]:

$$G_c^{II} = \int_{\delta_t=0}^{\infty} \tau d\delta_t \qquad (32)$$

with $\tau$ the shear stress across the shear band, and $\delta_t$ the relative sliding between both faces of the shear band. $G_c^{II}$ equals the area under the decohesion curves shown in Figure 2.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Relation between relative sliding at the discontinuity and friction coefficient</figcaption></figure>

Alternatively, a Coulomb criterion for local inception of the shear band has been used in the examples at the end of this paper. In this criterion nucleation starts when

$$\tau = \tau_{coh} + \tau_{fr} \qquad (33)$$

with

$$\tau_{coh} = c_0$$

the cohesive contribution and $c_0$ the cohesion. Decohesion is governed by the fracture energy $G_c^{II}$, similar to the Tresca-like criterion, cf. eq. (32). $\tau_{fr}$ is the
frictional contribution, which is defined as a function of the traction normal to the discontinuity and the effective friction coefficient $\mu$:

$$ \tau_{fr} = \mu \mathbf{n}_{\Gamma_d} \cdot \boldsymbol{\sigma} \cdot \mathbf{n}_{\Gamma_d} $$

The effective friction coefficient has a virgin value $\mu = \tan \phi$, with $\phi$ the friction angle. The vector $\mathbf{n}_{\Gamma_d}$ is such that it is normal to the critical plane where Coulomb's criterion for incipient shear failure is satisfied. A frictional softening relation models the microstructure evolution of the solid grains in the interface. The particular relation used in the example calculations is shown in Figure 3. In it, $\mu_c$ is the threshold value and $\delta_c$ the relative sliding at this value. It is assumed that $\delta_c = 2G_c^{II}/c_0$, which equals the value defined by the cohesive softening relation.

As with the Tresca-like criterion, it is assumed that the failure mode only involves sliding. Possible dilatancy effects are not included in the kinematics of the discontinuity. For this reason, one can also now suffice by adding only degrees-of-freedom to the finite element model that describe the discrete sliding mode. It is interesting to note that, unlike in non-associated plasticity, the resulting stiffness matrix remains symmetric.

A discrete equivalent of Darcy's relation is now defined for the fluid flow $\mathbf{q}_d$ at the discontinuity as:

$$ \mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d = -k_d(p^+ - p^-) = -k_d \tilde{p} \mid_{\mathbf{x} \in \Gamma_d} \qquad (34) $$

where $k_d$ is the permeability of the diaphragm that has been assumed to coincide with the displacement discontinuity $\Gamma_d$ and $p^+$ and $p^-$ are the pressures in the $\Omega^+$ and $\Omega^-$ domains, respectively. For an impervious boundary, $k_d = 0$, which implies that $\mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d = 0$ according to eq. (34). Conversely, ideal permeability requires that $k_d \to \infty$, so that $\mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d$ can only be bounded if $p^+ - p^- = 0$, which implies that no discontinuity can exist in the pressure field and the formulation of Armero and Callari [18] is retrieved.

# 5 Numerical elaboration

## 5.1 Weak forms

To arrive at the weak form of the balance equations, we multiply the momentum balance (4) and the mass balance (15) by test functions for the velocities of the skeleton and for the pressures. In the spirit of a standard Bubnov-Galerkin approach, they are assumed to be of the following format:

$$ \eta = \tilde{\eta} + \mathcal{H}_{\Gamma_d} \tilde{\eta} \qquad (35) $$

for the velocities, and

$$ \zeta = \bar{\zeta} + \mathcal{H}_{\Gamma_d} \tilde{\zeta} \qquad (36) $$
for the pressures. Substitution into eqs (4) and (15), and integrating over the domain $\Omega$ leads to the corresponding weak forms:

$$ \int_{\Omega} (\bar{\eta} + \mathcal{H}_{\Gamma_d} \tilde{\eta}) \cdot \left( \nabla \cdot \boldsymbol{\sigma} - \rho \frac{\partial \mathbf{v}_s}{\partial t} \right) d\Omega = 0 \quad (37) $$

and

$$ \int_{\Omega} (\bar{\zeta} + \mathcal{H}_{\Gamma_d} \tilde{\zeta}) \left( \alpha \nabla \cdot \mathbf{v}_s + n_f \nabla \cdot (\mathbf{v}_f - \mathbf{v}_s) + \frac{1}{Q} \frac{\partial p}{\partial t} \right) d\Omega = 0 \quad (38) $$

Using the standard procedure of applying the divergence theorem, using the external boundary conditions (16) and (17), eliminating the Heaviside functions by changing the integration domain from $\Omega$ to $\Omega^+$, eliminating the Dirac delta functions by transforming the volume integral into a surface integral, and introducing the shorter notation of a superimposed dot for $\partial/\partial t$, the balance equations take the form:

$$ \begin{aligned} & \int_{\Omega} \rho \bar{\eta} \cdot \dot{\mathbf{v}}_s d\Omega + \int_{\Omega^+} \rho \tilde{\eta} \cdot \dot{\mathbf{v}}_s d\Omega + \\ & \int_{\Omega} (\nabla \cdot \bar{\eta}) \cdot \boldsymbol{\sigma} d\Omega + \int_{\Omega^+} (\nabla \cdot \tilde{\eta}) \cdot \boldsymbol{\sigma} d\Omega + \\ & \int_{\Gamma_d} \tilde{\eta} \cdot \mathbf{t}_d d\Omega = \int_{\Gamma} (\bar{\eta} + \mathcal{H}_{\Gamma_d} \tilde{\eta}) \cdot \mathbf{t}_p d\Omega \end{aligned} \quad (39) $$

and

$$ \begin{aligned} & - \int_{\Omega} k_f \rho_f \nabla \bar{\zeta} \cdot \dot{\mathbf{v}}_s d\Omega - \int_{\Omega^+} k_f \rho_f \nabla \tilde{\zeta} \cdot \dot{\mathbf{v}}_s d\Omega \\ & - \int_{\Omega} \alpha \bar{\zeta} \nabla \cdot \mathbf{v}_s d\Omega - \int_{\Omega^+} \alpha \tilde{\zeta} \nabla \cdot \mathbf{v}_s d\Omega \\ & - \int_{\Omega} k_f \nabla \bar{\zeta} \cdot \nabla p d\Omega - \int_{\Omega^+} k_f \nabla \tilde{\zeta} \cdot \nabla p d\Omega - \int_{\Gamma_d} \tilde{\zeta} \mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d d\Gamma \\ & - \int_{\Omega} \bar{\zeta} Q^{-1} \dot{p} d\Omega - \int_{\Omega^+} \tilde{\zeta} Q^{-1} \dot{p} d\Omega = \int_{\Gamma} (\bar{\zeta} + \mathcal{H}_{\Gamma_d} \tilde{\zeta}) q_p d\Gamma \end{aligned} \quad (40) $$

where for the derivation of the latter equation also Darcy's relation (28) has been employed.

## 5.2 Discretisations

We now switch to matrix-vector notation and discretise the trial functions $\mathbf{v}_s$ and $p$ and the test functions $\eta$ and $\zeta$ as:

$$ \begin{aligned} \mathbf{v}_s &= \mathbf{N}(\dot{\bar{\mathbf{a}}} + \mathcal{H}_{\Gamma_d} \dot{\tilde{\mathbf{a}}}) \\ p &= \mathbf{H}(\bar{\mathbf{p}} + \mathcal{H}_{\Gamma_d} \tilde{\mathbf{p}}) \\ \eta &= \mathbf{N}(\bar{\mathbf{w}} + \mathcal{H}_{\Gamma_d} \tilde{\mathbf{w}}) \\ \zeta &= \mathbf{H}(\bar{\mathbf{z}} + \mathcal{H}_{\Gamma_d} \tilde{\mathbf{z}}) \end{aligned} \quad (41) $$
Eqs (41) can be inserted into eqs (39) and (40) to obtain the semi-discrete form under the requirement that the result holds for all admissible $\bar{\mathbf{w}}$, $\bar{\mathbf{z}}$, $\bar{\mathbf{w}}$ and $\tilde{\mathbf{z}}$.

The semi-discrete initial value problem is second order in time with respect to the displacement variables and first order for the fluid pore pressure. Yet, the same integration scheme will be used for both variables, in particular the Newmark method commonly used in structural dynamics. Let $\mathbf{y}$ denote an array which is a function of time. At the discrete time instant $t_n$ its value is $\mathbf{y}_n$. Under the assumption of time continuity, the updating equations of the Newmark method are:

$$ \mathbf{y}_{n+1} = \mathbf{y}_n + \Delta t \dot{\mathbf{y}}_n + \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{\mathbf{y}}_n + \beta \Delta t^2 \ddot{\mathbf{y}}_{n+1} \quad (42) $$

$$ \dot{\mathbf{y}}_{n+1} = \dot{\mathbf{y}}_n + (1 - \gamma) \Delta t \ddot{\mathbf{y}}_n + \gamma \Delta t \ddot{\mathbf{y}}_{n+1} \quad (43) $$

with $\beta, \gamma$ the parameters of the time integration scheme. The equations can be recast as:

$$ \begin{aligned} \ddot{\mathbf{y}}_{n+1} &= \alpha_0 (\mathbf{y}_{n+1} - \mathbf{y}_n) - \alpha_2 \dot{\mathbf{y}}_n - \alpha_4 \ddot{\mathbf{y}}_n \\ \dot{\mathbf{y}}_{n+1} &= \alpha_1 (\mathbf{y}_{n+1} - \mathbf{y}_n) - \alpha_3 \dot{\mathbf{y}}_n - \alpha_5 \ddot{\mathbf{y}}_n \end{aligned} \quad (44) $$

with

$$ \begin{aligned} \alpha_0 &= \frac{1}{\beta \Delta t^2}, & \alpha_2 &= \frac{1}{\beta \Delta t}, & \alpha_4 &= \frac{1}{2\beta} - 1 \\ \alpha_1 &= \frac{\gamma}{\beta \Delta t}, & \alpha_3 &= \frac{\gamma}{\beta} - 1, & \alpha_5 &= (\frac{\gamma}{2\beta} - 1) \Delta t \end{aligned} $$

For future use we also list the expressions for the variations that can be derived from expressions (44):

$$ \delta \ddot{\mathbf{y}}_{n+1} = \alpha_0 \delta \mathbf{y}_{n+1}, \quad \delta \dot{\mathbf{y}}_{n+1} = \alpha_1 \delta \mathbf{y}_{n+1} \quad (45) $$

Application of the time integration scheme (44) to semi-discrete balance equations results in a set of coupled, discrete equations, which is nonlinear. Therefore, an iterative solution procedure has to be applied within each time step $\Delta t$. When using the Newton-Raphson method, as has been done in the ensuing examples, and exploiting the variations defined in eqs (45), one obtains a sequence of linearised problems, which for implementation purposes are conveniently be cast in a matrix-vector format:

$$ \begin{bmatrix} \alpha_0 \mathbf{M}_{\bar{a}\bar{a}} + \mathbf{K}_{\bar{a}\bar{a}} & \alpha_0 \mathbf{M}_{\bar{a}\tilde{a}} + \mathbf{K}_{\bar{a}\tilde{a}} & \mathbf{K}_{\bar{a}\bar{p}} & \mathbf{K}_{\bar{a}\tilde{p}} \\ \alpha_0 \mathbf{M}_{\tilde{a}\bar{a}} + \mathbf{K}_{\tilde{a}\bar{a}} & \alpha_0 \mathbf{M}_{\tilde{a}\tilde{a}} + \mathbf{K}_{\tilde{a}\tilde{a}} & \mathbf{K}_{\tilde{a}\bar{p}} & \mathbf{K}_{\tilde{a}\tilde{p}} \\ \alpha_0 \mathbf{M}_{\bar{p}\bar{a}} + \alpha_1 \mathbf{K}_{\bar{a}\bar{p}}^T & \alpha_0 \mathbf{M}_{\bar{p}\tilde{a}} + \alpha_1 \mathbf{K}_{\bar{a}\tilde{p}}^T & \alpha_1 \mathbf{M}_{\bar{p}\bar{p}} + \mathbf{K}_{\bar{p}\bar{p}} & \alpha_1 \mathbf{M}_{\bar{p}\tilde{p}} + \mathbf{K}_{\bar{p}\tilde{p}} \\ \alpha_0 \mathbf{M}_{\tilde{p}\bar{a}} + \alpha_1 \mathbf{K}_{\tilde{a}\bar{p}}^T & \alpha_0 \mathbf{M}_{\tilde{p}\tilde{a}} + \alpha_1 \mathbf{K}_{\tilde{a}\tilde{p}}^T & \alpha_1 \mathbf{M}_{\tilde{p}\bar{p}} + \mathbf{K}_{\tilde{p}\bar{p}} & \alpha_1 \mathbf{M}_{\tilde{p}\tilde{p}} + \mathbf{K}_{\tilde{p}\tilde{p}} \end{bmatrix} \begin{bmatrix} d\bar{\mathbf{a}} \\ d\tilde{\mathbf{a}} \\ d\bar{\mathbf{p}} \\ d\tilde{\mathbf{p}} \end{bmatrix} = \begin{pmatrix} \mathbf{f}_{\bar{a}}^* \\ \mathbf{f}_{\tilde{a}}^* \\ \mathbf{f}_{\bar{p}}^* \\ \mathbf{f}_{\tilde{p}}^* \end{pmatrix} \quad (46) $$

with the external force vectors:
$$ \mathbf{f}_{\tilde{a}}^{\text{ext}} = \int_{\Gamma} \mathbf{N}^{\text{T}} \mathbf{t}_{p} \mathrm{d}\Gamma \ , \quad \mathbf{f}_{\tilde{a}}^{\text{ext}} = \int_{\Gamma} \mathcal{H}_{\Gamma_d} \mathbf{N}^{\text{T}} \mathbf{t}_{p} \mathrm{d}\Gamma $$

$$ \mathbf{f}_{\tilde{p}}^{\text{ext}} = \int_{\Gamma} \mathbf{H}^{\text{T}} q_{p} \mathrm{d}\Gamma \ , \quad \mathbf{f}_{\tilde{p}}^{\text{ext}} = \int_{\Gamma} \mathcal{H}_{\Gamma_d} \mathbf{H}^{\text{T}} q_{p} \mathrm{d}\Gamma $$

the internal force vectors:

$$ \mathbf{f}_{\tilde{a}}^{\text{int}} = \int_{\Omega} \mathbf{B}^{\text{T}} \boldsymbol{\sigma} \mathrm{d}\Omega $$

$$ \mathbf{f}_{\tilde{a}}^{\text{int}} = \int_{\Omega^{+}} \mathbf{B}^{\text{T}} \boldsymbol{\sigma} \mathrm{d}\Omega + \int_{\Gamma_d} \mathbf{N}^{\text{T}} \mathbf{t}_{d} \mathrm{d}\Gamma $$

with $\mathbf{B} = \nabla \mathbf{N}$ and, for two dimensions, $\mathbf{m} = [1, 1, 0]$. The mass matrices:

$$ \mathbf{M}_{\tilde{a}\tilde{a}} = \int_{\Omega} \rho \mathbf{N}^{\text{T}} \mathbf{N} \mathrm{d}\Omega \ , \quad \mathbf{M}_{\tilde{a}\tilde{a}} = \mathbf{M}_{\tilde{a}\tilde{a}} = \mathbf{M}_{\tilde{a}\tilde{a}} = \int_{\Omega^{+}} \rho \mathbf{N}^{\text{T}} \mathbf{N} \mathrm{d}\Omega $$

$$ \mathbf{M}_{\tilde{p}\tilde{a}} = - \int_{\Omega} k_{f} \rho_{f} \nabla \mathbf{H}^{\text{T}} \mathbf{N} \mathrm{d}\Omega \ , \quad \mathbf{M}_{\tilde{p}\tilde{a}} = \mathbf{M}_{\tilde{p}\tilde{a}} = \mathbf{M}_{\tilde{p}\tilde{a}} = - \int_{\Omega^{+}} k_{f} \rho_{f} \nabla \mathbf{H}^{\text{T}} \mathbf{N} \mathrm{d}\Omega $$

$$ \mathbf{M}_{\tilde{p}\tilde{p}} = - \int_{\Omega} Q^{-1} \mathbf{H}^{\text{T}} \mathbf{H} \mathrm{d}\Omega \ , \quad \mathbf{M}_{\tilde{p}\tilde{p}} = \mathbf{M}_{\tilde{p}\tilde{p}} = \mathbf{M}_{\tilde{p}\tilde{p}} = - \int_{\Omega^{+}} Q^{-1} \mathbf{H}^{\text{T}} \mathbf{H} \mathrm{d}\Omega $$

the stiffness matrices:

$$ \mathbf{K}_{\tilde{a}\tilde{p}} = - \int_{\Omega} \alpha \mathbf{B}^{\text{T}} \mathbf{m} \mathrm{H} \mathrm{d}\Omega \ , \quad \mathbf{K}_{\tilde{a}\tilde{p}} = \mathbf{K}_{\tilde{a}\tilde{p}} = \mathbf{K}_{\tilde{a}\tilde{p}} = - \int_{\Omega^{+}} \alpha \mathbf{B}^{\text{T}} \mathbf{m} \mathrm{H} \mathrm{d}\Omega $$

$$ \mathbf{K}_{\tilde{p}\tilde{p}} = - \int_{\Omega} k_{f} \nabla \mathbf{H}^{\text{T}} \nabla \mathbf{H} \mathrm{d}\Omega \ , \quad \mathbf{K}_{\tilde{p}\tilde{p}} = \mathbf{K}_{\tilde{p}\tilde{p}} = - \int_{\Omega^{+}} k_{f} \nabla \mathbf{H}^{\text{T}} \nabla \mathbf{H} \mathrm{d}\Omega $$

$$ \mathbf{K}_{\tilde{p}\tilde{p}} = - \int_{\Omega^{+}} k_{f} \nabla \mathbf{H}^{\text{T}} \nabla \mathbf{H} \mathrm{d}\Omega - \int_{\Gamma_{d}} k_{d} \mathbf{H}^{\text{T}} \mathbf{H} \mathrm{d}\Gamma $$

$$ \mathbf{K}_{\tilde{a}\tilde{a}} = \int_{\Omega} \mathbf{B}^{\text{T}} \mathbf{D} \mathbf{B} \mathrm{d}\Omega \ , \quad \mathbf{K}_{\tilde{a}\tilde{a}} = \mathbf{K}_{\tilde{a}\tilde{a}} = \int_{\Omega^{+}} \mathbf{B}^{\text{T}} \mathbf{D} \mathbf{B} \mathrm{d}\Omega $$

$$ \mathbf{K}_{\tilde{a}\tilde{a}} = \int_{\Omega^{+}} \mathbf{B}^{\text{T}} \mathbf{D} \mathbf{B} \mathrm{d}\Omega + \int_{\Gamma_{d}} \mathbf{N}^{\text{T}} \mathbf{T} \mathbf{N} \mathrm{d}\Gamma $$

and the arrays at the right-hand side:

$$ \mathbf{f}_{\tilde{a}}^{*} = \mathbf{f}_{\tilde{a}}^{\text{ext}} - (\mathbf{f}_{\tilde{a}}^{\text{int}})^{i} - \alpha_{0} \mathbf{M}_{\tilde{a}\tilde{a}} \ddot{\mathbf{a}}_{n+1}^{i} - \alpha_{0} \mathbf{M}_{\tilde{a}\tilde{a}} \ddot{\mathbf{a}}_{n+1}^{i} $$

$$ \mathbf{f}_{\tilde{a}}^{*} = \mathbf{f}_{\tilde{a}}^{\text{ext}} - (\mathbf{f}_{\tilde{a}}^{\text{int}})^{i} - \alpha_{0} \mathbf{M}_{\tilde{a}\tilde{a}} \ddot{\mathbf{a}}_{n+1}^{i} - \alpha_{0} \mathbf{M}_{\tilde{a}\tilde{a}} \ddot{\mathbf{a}}_{n+1}^{i} $$

$$ \mathbf{f}_{\tilde{p}}^{*} = \mathbf{f}_{\tilde{p}}^{\text{ext}} - \alpha_{0} \mathbf{M}_{\tilde{p}\tilde{a}} \ddot{\mathbf{a}}_{n+1}^{i} - \alpha_{0} \mathbf{M}_{\tilde{p}\tilde{a}} \ddot{\mathbf{a}}_{n+1}^{i} - \alpha_{1} \mathbf{M}_{\tilde{p}\tilde{p}} \dot{\mathbf{p}}_{n+1}^{i} - \alpha_{1} \mathbf{M}_{\tilde{p}\tilde{p}} \dot{\mathbf{p}}_{n+1}^{i} \\ - \alpha_{1} \mathbf{K}_{\tilde{a}\tilde{p}}^{\text{T}} \dot{\mathbf{a}}_{n+1}^{i} - \alpha_{1} \mathbf{K}_{\tilde{a}\tilde{p}}^{\text{T}} \dot{\mathbf{a}}_{n+1}^{i} $$

$$ \mathbf{f}_{\tilde{p}}^{*} = \mathbf{f}_{\tilde{p}}^{\text{ext}} + -\alpha_{0} \mathbf{M}_{\tilde{p}\tilde{a}} \ddot{\mathbf{a}}_{n+1}^{i} - \alpha_{0} \mathbf{M}_{\tilde{p}\tilde{a}} \dot{\mathbf{a}}_{n+1}^{i} - \alpha_{1} \mathbf{M}_{\tilde{p}\tilde{p}} \dot{\mathbf{p}}_{n+1}^{i} - \alpha_{1} \mathbf{M}_{\tilde{p}\tilde{p}} \dot{\mathbf{p}}_{n+1}^{i} \\ - \alpha_{1} \mathbf{K}_{\tilde{a}\tilde{p}}^{\text{T}} \dot{\mathbf{a}}_{n+1}^{i} - \alpha_{1} \mathbf{K}_{\tilde{a}\tilde{p}}^{\text{T}} \dot{\mathbf{a}}_{n+1}^{i} $$

where the superscript signifies that the corresponding quantity has to be evaluated at iteration $i$. The quantities $\ddot{\mathbf{a}}_{n+1}^{i}$, $\dot{\mathbf{a}}_{n+1}^{i}$, $\ddot{\mathbf{a}}_{n+1}^{i}$, $\dot{\mathbf{a}}_{n+1}^{i}$, $\dot{\mathbf{p}}_{n+1}^{i}$, $\dot{\mathbf{p}}_{n+1}^{i}$ are evaluated using eqs (44).
The stiffness matrix of eq. (46) is not symmetric. Symmetry can be restored by multiplying the third and the fourth row of submatrices by $\alpha_1^{-1}$ and omitting the contributions in the tangent stiffness matrix that are due to the dynamic seepage term – the submatrices $M_{\bar{p}\bar{a}}$, $M_{\bar{p}\tilde{a}}$, $M_{\tilde{p}\bar{a}}$ and $M_{\tilde{p}\tilde{a}}$. Since the corresponding terms are retained in the right-hand side, the results are not affected, only the convergence speed of the iterative procedure.

### 5.3 Stress computation at the tip

The nucleation criterion requires the determination of the stresses at the tip of the discontinuity. Unfortunately, the stresses vary strongly in the vicinity of the tip and an accurate estimate of them is difficult to obtain. In the present case, the stress in the bulk of the specimen is almost homogeneous except for a small area around the tip, which exacerbates the problem. Following Wells [5] and Jirasek [23] we use a smoothing of the stresses around the tip and compute the stress at the tip by the following nonlocal-like procedure:

$$ \sigma_{tip} = \frac{\int_{\Omega} w \sigma \, d\Omega}{\int_{\Omega} w \, d\Omega} \qquad (47) $$

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 nearly homogeneous stress state in the specimen, 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 discretisation – 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.

## 6 Example calculations

All results of the computations are based for the same two-dimensional specimen with a width $w = 0.04$ m and a height $H = 0.1$ m, see also Figure 4, which is loaded under plane-strain conditions. The sides are traction free and the external loading is applied via an imposed constant velocity $V_0 = -10^{-3}$ m/s. Undrained conditions have been imposed on the entire boundary of the specimen, because fast transient phenomena have been considered. The solid constituent is assumed to behave in a linear elastic manner
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Geometry and boundary conditions</figcaption></figure>

with a Young's modulus $E = 20$ GPa and a Poisson's ratio $\nu = 0.35$. The absolute mass densities are $\rho'_s = \rho_s/n_s = 2000$ kg/m$^3$ for the solid phase and $\rho'_f = \rho_f/n_f = 1000$ kg/m$^3$ for the fluid phase, while the fluid fraction $n_f = 0.3$. The Biot coefficient $\alpha$ has been set equal to 1, the Biot modulus has been assigned a value $Q = 5.0$ GPa, while the bulk material was assumed to have a permeability $k_f = 10^{-14}$ m$^3$/Ns. The permeability of the diaphragm was assigned a value $k_d = 0.5 \cdot 10^{-14}$ m$^2$/Ns. Shear-band formation was triggered by a small imperfection, see Figure 4.

A structured mesh has been used and consists of 5841 four-noded elements with equal (bilinear) interpolations for the displacements and the pressure. The simulation is started using a time step of 0.4 s, which is small enough to accurately follow the pressure evolution and the near quasi-static behavior of the solid skeleton before the onset of the shear band. When the shear band starts to propagate, the phenomenon becomes dynamic and the time step size is reduced severely in order to properly capture the propagation of the stress wave. The parameters of the Newmark scheme are $\gamma = 0.5$ and $\beta = 0.25$.

## 6.1 Tresca-like initiation criterion

The simulation for the Tresca-like nucleation criterion from which most of the results derive, has been obtained with the following parameters: the time step size during the shear-band formation equals 0.2 µs, nucleation traction $\tau_c = 50$ MPa and mode-II fracture energy $G_c^{II} = 750$ J/m$^2$. The dynamic seepage term has not been taken into account in the analysis.

The evolution of the pore pressure field following the time $t_0$ at which the shear band starts to propagate, is shown in Figure 5. The scale has been chosen such that the white regions on the picture have a pore pressure below the cavitation pressure (here: $-10^{+5}\text{Pa}$). One observes that, initially, cavitation occurs only in the close vicinity of the discontinuity. When the shear band tip reaches the centre of the specimen, the level of the pore pressure above the
<figure><img src="image_6.png" /><figcaption>Fig. 5. Evolution of the pressure field for the Tresca criterion</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Pressure field near the process zone and tractions at the discontinuity. The scale of the pressures is equal to that in Figure 5. The magnitudes of the tractions are proportional to the lengths of the bars</figcaption></figure>
discontinuity increases and pore pressures that exceed the cavitation pressure develop over a larger region below the interface. Subsequently, this region follows the tip of the shear band, and when it reaches the right boundary of the specimen, the cavitation phenomena extend over the entire specimen. Because of the limitations of the model – a gas phase has not been modelled separately – the physical interpretation of the numerical results at this advanced stage of shear band propagation becomes questionable.

As illustrated by Figure 6, it seems that in the present simulations that utilise the Tresca-like criterion, cavitation is a consequence of local elastic unloading behind the process zone, where cohesive softening takes place. Indeed, Figure 6 shows the pressure field as well as the values of the shear tractions in the cohesive interface (with bars orthogonal to the discontinuity). The cavitation front appears to coincide with the transition zone between the damaged and the intact parts of the interface.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Evolution of the pressure field for the Coulomb criterion with $\frac{\delta_c}{\Delta U_{ini}} = 0.024$</figcaption></figure>

## 6.2 Coulomb initiation criterion

For the Coulomb criterion, the following model parameters have been used: cohesion $c_0 = 10.0$ MPa, friction angle $\phi = 30^{\circ}$. The threshold value in the
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Evolution of the pressure field for the Coulomb criterion with $\frac{\delta_c}{\Delta U_{ini}} = 0.16$</figcaption></figure>

frictional softening law is $\mu_c = 0.2 \tan \phi$. A parametric study has been carried out with respect to the influence of $\frac{\delta_c}{\Delta U_{ini}}$, where $\delta_c$ is the tangential displacement jump when the cohesive part of the traction has vanished, and $\Delta U_{ini}$ the value of the prescribed displacement at the top of the specimen at shear-band initiation. The simulations have been carried out using a time step size of 0.25 $\mu$s.

A first simulation has been carried obtained for a fracture energy $G_c^{II} = 15$ J/m$^2$. In this case $\frac{\delta_c}{\Delta U_{ini}}$ equals 0.024. The evolution of the pressure field is shown in Figure 7. The results are quite similar to those obtained with the Tresca criterion since the values of $\frac{\delta_c}{\Delta U_{ini}}$ are close. Figure 10 shows that for this case the process zone, i.e. where the tractions do not vanish, is rather small. As a consequence, local unloadings are obtained behind this process zone and cavitation occurs only in a small zone behind the shear-band tip.

Increasing the value of $\frac{\delta_c}{\Delta U_{ini}}$ to 0.16, the length of the process zone becomes approximately equal to that of the shear band, see Figure 10. Now, local unloadings are not observed, but a global unloading occurs that can be associated with the strain localisation inside the shear band and the softening of the interface at the end of propagation. As a consequence, the pressure evolution on Figure 8 has no local cavitation zone. Indeed, cavitation is first obtained around the initiation locus and, subsequently, in the entire specimen.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Evolution of the pressure field for the Coulomb criterion with $\frac{\delta_c}{\Delta U_{ini}} = 0.24$</figcaption></figure>

Instability patterns are observed on the last two figures that plot the pressure distribution.

Such instabilities, which are due to friction, are observed at a larger scale when $\frac{\delta_c}{\Delta U_{ini}}$ is increased further up to 0.24. Now, the shear band propagates because the stress is higher than the material strength, but the mechanical energy is not sufficient to damage the interface. Consequently, no localisation is obtained, but only frictional instabilities are observed. No cavitation is induced because of the absence of strain localisation.

# 7 Concluding remarks

In this contribution a numerical model has been elaborated which can capture discontinuities, e.g. cracks or shear bands, in a fluid-saturated medium. The representation of the discontinuity is truly discrete and unbiased by the discretisation. Moreover, the constitutive relations for the bulk and for the discontinuity can be specified independently, for the solid phase as well as for the fluid phase. Example calculations of dynamic shear band propagation have been presented with a Tresca-like and a Coulomb criterion for shear band initiation. The results show that the propagation of the shear band is strongly
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Shear traction distribution at the end of the shear-band propagation for different values of $\frac{\delta_c}{\Delta U_{ini}}$</figcaption></figure>

influenced by the constitutive assumptions in the discontinuity. Indeed, as highlighted by the results of the calculations, the cavitation phenomenon is triggered by unloading of the solid skeleton, which is a direct consequence of strain localisation and strongly depends on the constitutive model for the discontinuity.

## 8 Acknowledgements

This work was supported by the European Commission under contract HPRN-CT-2002-00198.

## References

1. de Borst R (2004) Damage, material instabilities, and failure. In: Encyclopedia of Computational Mechanics. Wiley, Chichester

2. Babuska I, Melenk JM (1997) International Journal for Numerical Methods in Engineering 40:727–758

3. Belytschko T, Black T (1999) International Journal for Numerical Methods in Engineering 45:601–620
4. Moës N, Dolbow J, Belytschko T (1999) International Journal for Numerical Methods in Engineering 46:131–150

5. Wells GN, Sluys LJ (2001) International Journal for Numerical and Analytical Methods in Geomechanics 25:691–709

6. Wells GN, Sluys LJ, de Borst R (2002) International Journal for Numerical Methods in Engineering 53:1235–1256

7. Wells GN, de Borst R, Sluys LJ (2002) International Journal for Numerical Methods in Engineering 54:1333–1355

8. Remmers JJC, de Borst R, Needleman A (2003) Computational Mechanics 31:69–77

9. Samaniego E, Belytschko T (2005) International Journal for Numerical Methods in Engineering 62:1857–1872

10. Réthoré J, Gravouil A, Combescure A (2005) International Journal for Numerical Methods in Engineering 63:631–659

11. Réthoré J, Gravouil A, Combescure A (2005) International Journal for Numerical Methods in Engineering 64:260–284

12. Areias PMA, Belytschko T (2006) International Journal for Numerical Methods in Engineering 66:878–910

13. Terzaghi K (1943) Theoretical Soil Mechanics John Wiley & Sons, New York

14. Biot MA (1965) Mechanics of Incremental Deformations John Wiley & Sons, Chichester

15. Jouanna P, Abellan MA (1995) Generalized approach to heterogeneous media In: Modern Issues in Non-Saturated Soils. Springer-Verlag, Wien - New York, 1–128

16. Lewis RW, Schrefler BA (1998) The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, Second Edition John Wiley & Sons, Chichester

17. de Borst R, Réthoré J, Abellan MA (2006) Archive of Applied Mechanics 75:595–606

18. Armero F, Callari C (1999) International Journal for Numerical Methods in Engineering 46:1673–1698

19. Larsson J, Larsson R (2000) Mechanics of Cohesive-frictional Materials 5:565–582

20. Schrefler, B, Scotta R (2001) Computer Methods in Applied Mechanics and Engineering 190:3223–3246

21. Feenstra PH, de Borst R (1996) International Journal of Solids and Structures 33:707–730

22. Abellan MA, de Borst R (2006) Computer Methods in Applied Mechanics and Engineering 195:5011–5019

23. Jirasek M (1998) Embedded crack models for concrete fracture In: Computational Modelling of Concrete Structures. Balkema, Rotterdam, 291–300