# A discrete model for the dynamic propagation of shear bands 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 dynamic propagation of shear bands in a fluid-saturated medium. International Journal for Numerical and Analytical Methods in Geomechanics, 2007, Advanced Models for Fracture in Quasi-Brittle Materials, 31 (2), pp.347-370. 10.1002/nag.575 . hal-00380798

HAL Id: hal-00380798

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

Submitted on 15 Jul 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 dynamic propagation of shear bands
in a fluid-saturated medium

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

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

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

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

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 is 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. Finally, a limited quantitative study is made on the importance of including a so-called dynamic seepage term in Darcy’s relation when considering dynamic shear banding.

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

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

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

<footnote>Contract/grant sponsor: European Commission; contract/grant number: HPRN-CT-2002-00198</footnote>
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 character, from elliptic to hyperbolic for quasi-static loadings, and from hyperbolic to
elliptic for dynamic loadings. Unless a regularization is applied, the resulting initial/boundary
value problem becomes ill-posed, resulting in numerical solutions that depend severely on the
discretization [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–12]. While preserving the
original discretization, 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 discretization.

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 recognized 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 specialized 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 localization on cavitation in a fluid-saturated porous medium is highlighted. Finally, a quantitative study is carried out regarding the role of the so-called dynamic seepage term in Darcy's law.

## 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 supposing 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{\boldsymbol{p}}_{\pi} + \rho_{\pi} \boldsymbol{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, $\boldsymbol{g}$ is the gravity acceleration and $\hat{\boldsymbol{p}}_{\pi}$ is the source of momentum for constituent $\pi$ from the other constituent, which takes into account the possible local drag interaction between the solid and the fluid. Evidently, the latter source terms must satisfy the momentum production constraint:

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

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

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

Adding both momentum balances, and taking into account Equation (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: $\partial \mathbf{v}_s / \partial t \approx \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. At the end of the paper, we will assess its accuracy, also in the presence of a discontinuity.

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

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

$$ \sum_{\pi=s,f} n_{\pi} = 1 \quad (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 \quad (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 n_s \partial \rho_s}{K_t \rho_s \partial t} \quad (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, $\alpha = 1 - K_t/K_s$ [16], this equation can be rewritten as

$$ (\alpha - 1)\nabla \cdot \mathbf{v}_s = \frac{n_s \partial \rho_s}{\rho_s \partial t} \quad (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} \mathrm{d}p = \frac{n_f}{\rho_f} \mathrm{d}\rho_f \quad (13) $$

with the overall compressibility or Biot modulus

$$ \frac{1}{Q} = \frac{\alpha - n_f}{K_s} + \frac{n_f}{K_f} \quad (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, Equation (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 \quad (15) $$

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

$$ \mathbf{n}_{\Gamma} \cdot \mathbf{\sigma} = \mathbf{t}_{\mathrm{p}}, \quad \mathbf{v} = \mathbf{v}_{\mathrm{p}} \quad (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 \quad (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} \quad (18) $$

close the initial value problem.
### 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 Black [3] and Moës et al. [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 $v$ can be interpolated as

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

with $\tilde{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 kth 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 $v$ can be written as the sum of two continuous velocity fields $\tilde{v}$ and $\tilde{v}$:

$$v = \tilde{v} + \mathcal{H}_{\Gamma_d} \tilde{v} \quad (20)$$

where $\mathcal{H}_{\Gamma_d}$ is the Heaviside step function centred at the discontinuity. The decomposition in Equation (20) has a structure similar to the interpolation in Equation (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 \tilde{v}_s + \mathcal{H}_{\Gamma_d} \nabla^s \tilde{v}_s + \delta_{\Gamma_d} (\tilde{v}_s \otimes 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 Equation (20), we have

$$p = \tilde{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 = \tilde{p}$) and is also different from that of

<figure><img src="image_2.png" /><figcaption>Figure 1. Body composed of continuous displacement fields at each side of the discontinuity $\Gamma_d$.</figcaption></figure>
Larsson and Larsson [19], who assume that a regularized 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 Equation (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) $$

# 4. CONSTITUTIVE EQUATIONS

## 4.1. Models for the bulk

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

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

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

$$ d\boldsymbol{\sigma}_s = \mathbf{D}^{\text{tan}} : d\boldsymbol{\varepsilon}_s \quad (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 \quad (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) \quad (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 $\partial \mathbf{v}_s / \partial t \approx \partial \mathbf{v}_f / \partial t$ for relatively slow dynamic loadings, Equation (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) \quad (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) \quad (29) $$

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

$$ \dot{\mathbf{t}}_d = \mathbf{T} \dot{\boldsymbol{\delta}} \quad (30) $$
with T the material tangent stiffness matrix of the discrete traction-separation law:

$$ T = \frac{\partial t_d}{\partial \delta} + \frac{\partial t_d}{\partial \kappa} \frac{\partial \kappa}{\partial \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 maximizes the shear stress. In this orientation, the shear stress $\tau = |\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, $\mathcal{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]:

$$ \mathcal{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. $\mathcal{G}_c^{II}$ equals the area under the decohesion curves shown in Figure 2.

<figure><img src="image_3.png" /><figcaption>Figure 2. Relation between relative sliding at the discontinuity and shear tractions.</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_{\text{coh}} + \tau_{\text{fr}} \qquad (33) $$

with

$$ \tau_{\text{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. Equation (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_{\text{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}|_{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 Equation (34). Conversely, ideal permeability requires that $k_d \to \infty$, so that $\mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d$ can only be

<figure><img src="image_4.png" /><figcaption>Figure 3. Relation between relative sliding at the discontinuity and friction coefficient.</figcaption></figure>
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 = \tilde{\zeta} + \mathcal{H}_{\Gamma_d} \tilde{\zeta} \qquad (36) $$

for the pressures. Substitution into Equations (4) and (15), and integrating over the domain $\Omega$ leads to the corresponding weak forms:

$$ \int_{\Omega} (\tilde{\eta} + \mathcal{H}_{\Gamma_d} \tilde{\eta}) \cdot (\nabla \cdot \sigma - \rho \frac{\partial v_s}{\partial t}) d\Omega = 0 \qquad (37) $$

and

$$ \int_{\Omega} (\tilde{\zeta} + \mathcal{H}_{\Gamma_d} \tilde{\zeta}) (\alpha \nabla \cdot \mathbf{v}_s + n_f \nabla \cdot (\mathbf{v}_f - \mathbf{v}_s) + \frac{1}{Q} \frac{\partial p}{\partial t}) d\Omega = 0 \qquad (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 \tilde{\eta} \cdot \dot{\mathbf{v}}_s \, d\Omega + \int_{\Omega^+} \rho \tilde{\eta} \cdot \dot{\mathbf{v}}_s \, d\Omega + \int_{\Omega} (\nabla \cdot \tilde{\eta}) \cdot \sigma \, d\Omega + \int_{\Omega^+} (\nabla \cdot \tilde{\eta}) \cdot \sigma \, d\Omega + \int_{\Gamma_d} \tilde{\eta} \cdot \mathbf{t}_d \, d\Omega \\ & = \int_{\Gamma} (\tilde{\eta} + \mathcal{H}_{\Gamma_d} \tilde{\eta}) \cdot \mathbf{t}_p \, d\Omega \end{aligned} \qquad (39) $$

and

$$ \begin{aligned} & - \int_{\Omega} k_f \rho_f \nabla \tilde{\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_{\tilde{\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 \tilde{\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} \tilde{\zeta} Q^{-1} \dot{p} \, d\Omega - \int_{\Omega^+} \tilde{\zeta} Q^{-1} \dot{p} \, d\Omega \\ & = \int_{\Gamma} (\tilde{\zeta} + \mathcal{H}_{\Gamma_d} \tilde{\zeta}) q_p \, d\Gamma \end{aligned} \qquad (40) $$

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

We now switch to matrix-vector notation and discretize the trial functions $v_s$ and $p$ and the test functions $\eta$ and $\zeta$ as

$$ 
\begin{aligned}
v_s &= N(\dot{\tilde{a}} + \mathcal{H}_{\Gamma_d} \dot{\tilde{a}}) \\
p &= H(\tilde{p} + \mathcal{H}_{\Gamma_d} \tilde{p}) \\
\eta &= N(\tilde{w} + \mathcal{H}_{\Gamma_d} \tilde{w}) \\
\zeta &= H(\tilde{z} + \mathcal{H}_{\Gamma_d} \tilde{z})
\end{aligned} 
\quad (41) $$

Inserting Equations (41) into Equations (39) and (40) and requiring that the result holds for all admissible $\tilde{w}, \tilde{z}, \tilde{w}$ and $\tilde{z}$ gives:

$$ 
\begin{aligned}
\mathbf{M}_{\tilde{a}\tilde{a}}\ddot{\tilde{a}} + \mathbf{M}_{\tilde{a}\tilde{a}}\dot{\ddot{a}} + \mathbf{f}_{\tilde{a}}^{\text{int}} &= \mathbf{f}_{\tilde{a}}^{\text{ext}} \\
\mathbf{M}_{\tilde{a}\tilde{a}}\ddot{\tilde{a}} + \mathbf{M}_{\tilde{a}\tilde{a}}\dot{\ddot{a}} + \mathbf{f}_{\tilde{a}}^{\text{int}} &= \mathbf{f}_{\tilde{a}}^{\text{ext}} \\
\mathbf{M}_{\tilde{p}\tilde{a}}\ddot{\tilde{a}} + \mathbf{M}_{\tilde{p}\tilde{a}}\dot{\ddot{a}} + \mathbf{M}_{\tilde{p}\tilde{p}}\dot{\tilde{p}} + \mathbf{M}_{\tilde{p}\tilde{p}}\dot{\dot{p}} + \mathbf{K}_{\tilde{a}\tilde{p}}^{\text{T}}\dot{\tilde{a}} + \mathbf{K}_{\tilde{a}\tilde{p}}^{\text{T}}\dot{\dot{a}} + \mathbf{K}_{\tilde{p}\tilde{p}}\tilde{p} + \mathbf{K}_{\tilde{p}\tilde{p}}\tilde{\tilde{p}} &= \mathbf{f}_{\tilde{p}}^{\text{ext}} \\
\mathbf{M}_{\tilde{p}\tilde{a}}\ddot{\tilde{a}} + \mathbf{M}_{\tilde{p}\tilde{a}}\dot{\ddot{a}} + \mathbf{M}_{\tilde{p}\tilde{p}}\dot{\tilde{p}} + \mathbf{M}_{\tilde{p}\tilde{p}}\dot{\dot{p}} + \mathbf{K}_{\tilde{a}\tilde{p}}^{\text{T}}\dot{\tilde{a}} + \mathbf{K}_{\tilde{a}\tilde{p}}^{\text{T}}\dot{\dot{a}} + \mathbf{K}_{\tilde{p}\tilde{p}}\tilde{p} + \mathbf{K}_{\tilde{p}\tilde{p}}\tilde{\tilde{p}} &= \mathbf{f}_{\tilde{p}}^{\text{ext}}
\end{aligned} 
\quad (42) $$

with the mass matrices:

$$ 
\begin{aligned}
\mathbf{M}_{\tilde{a}\tilde{a}} &= \int_{\Omega} \rho \mathbf{N}^{\text{T}} \mathbf{N} \, \mathrm{d}\Omega, & \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, & \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, & \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
\end{aligned} 
\quad (43) $$

the stiffness matrices:

$$ 
\begin{aligned}
\mathbf{K}_{\tilde{a}\tilde{p}} &= -\int_{\Omega^{+}} \alpha \mathbf{B}^{\text{T}} \mathbf{m} \mathbf{H} \, \mathrm{d}\Omega, & \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} \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, & \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
\end{aligned} 
\quad (44) $$

the external force vectors:

$$ 
\begin{aligned}
\mathbf{f}_{\tilde{a}}^{\text{ext}} &= \int_{\Gamma} \mathbf{N}^{\text{T}} \mathbf{t}_{\mathrm{p}} \, \mathrm{d}\Gamma, & \mathbf{f}_{\tilde{a}}^{\text{ext}} &= \int_{\Gamma} \mathcal{H}_{\Gamma_d} \mathbf{N}^{\text{T}} \mathbf{t}_{\mathrm{p}} \, \mathrm{d}\Gamma \\
\mathbf{f}_{\tilde{p}}^{\text{ext}} &= \int_{\Gamma} \mathbf{H}^{\text{T}} \mathbf{q}_{\mathrm{p}} \, \mathrm{d}\Gamma, & \mathbf{f}_{\tilde{p}}^{\text{ext}} &= \int_{\Gamma} \mathcal{H}_{\Gamma_d} \mathbf{H}^{\text{T}} \mathbf{q}_{\mathrm{p}} \, \mathrm{d}\Gamma
\end{aligned} 
\quad (45) $$
and the internal force vectors:

$$ \mathbf{f}_{\bar{a}}^{\text{int}} = \int_{\Omega} \mathbf{B}^{\text{T}} \boldsymbol{\sigma} \, \mathrm{d}\Omega \qquad (46) $$

$$ \mathbf{f}_{\bar{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 N$ and, for two dimensions, $\mathbf{m} = [1, 1, 0]$.

### 5.3. Temporal integration and linearization

The semi-discrete initial value problem (42) 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 + (\frac{1}{2} - \beta) \Delta t^2 \ddot{\mathbf{y}}_n + \beta \Delta t^2 \ddot{\mathbf{y}}_{n+1} \quad (47) $$

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

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

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 &= \left(\frac{\gamma}{2\beta} - 1\right) \Delta t \end{aligned} $$

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

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

Application of the time integration scheme (49) to Equations (42) 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 Equations (50), one obtains a sequence of linearized problems, which for implementation purposes are conveniently be cast in a matrix-vector format:

$$ \left[ \begin{array}{llll} \alpha_0 \mathbf{M}_{\bar{a}\bar{a}} + \mathbf{K}_{\bar{a}\bar{a}} & \alpha_0 \mathbf{M}_{\bar{a}\bar{a}} + \mathbf{K}_{\bar{a}\bar{a}} & \mathbf{K}_{\bar{a}\bar{p}} & \mathbf{K}_{\bar{a}\bar{p}} \\ \alpha_0 \mathbf{M}_{\bar{a}\bar{a}} + \mathbf{K}_{\bar{a}\bar{a}} & \alpha_0 \mathbf{M}_{\bar{a}\bar{a}} + \mathbf{K}_{\bar{a}\bar{a}} & \mathbf{K}_{\bar{a}\bar{p}} & \mathbf{K}_{\bar{a}\bar{p}} \\ \alpha_0 \mathbf{M}_{\bar{p}\bar{a}} + \alpha_1 \mathbf{K}_{\bar{a}\bar{p}}^{\text{T}} & \alpha_0 \mathbf{M}_{\bar{p}\bar{a}} + \alpha_1 \mathbf{K}_{\bar{a}\bar{p}}^{\text{T}} & \alpha_1 \mathbf{M}_{\bar{p}\bar{p}} + \mathbf{K}_{\bar{p}\bar{p}} & \alpha_1 \mathbf{M}_{\bar{p}\bar{p}} + \mathbf{K}_{\bar{p}\bar{p}} \\ \alpha_0 \mathbf{M}_{\bar{p}\bar{a}} + \alpha_1 \mathbf{K}_{\bar{a}\bar{p}}^{\text{T}} & \alpha_0 \mathbf{M}_{\bar{p}\bar{a}} + \alpha_1 \mathbf{K}_{\bar{a}\bar{p}}^{\text{T}} & \alpha_1 \mathbf{M}_{\bar{p}\bar{p}} + \mathbf{K}_{\bar{p}\bar{p}} & \alpha_1 \mathbf{M}_{\bar{p}\bar{p}} + \mathbf{K}_{\bar{p}\bar{p}} \end{array} \right] \begin{pmatrix} \mathrm{d}\bar{a} \\ \mathrm{d}\bar{a} \\ \mathrm{d}\bar{p} \\ \mathrm{d}\bar{p} \end{pmatrix} = \begin{pmatrix} \mathbf{f}_{\bar{a}}^{*} \\ \mathbf{f}_{\bar{a}}^{*} \\ \mathbf{f}_{\bar{p}}^{*} \\ \mathbf{f}_{\bar{p}}^{*} \end{pmatrix} \quad (51) $$
with the stiffness matrices

$$ 
\begin{align} 
\mathbf{K}_{\tilde{a}\tilde{a}} &= \int_{\Omega} \mathbf{B}^{\mathrm{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}^{\mathrm{T}} \mathbf{D}\mathbf{B} \, \mathrm{d}\Omega \tag{52} \\ 
\mathbf{K}_{\tilde{a}\tilde{a}} &= \int_{\Omega^{+}} \mathbf{B}^{\mathrm{T}} \mathbf{D}\mathbf{B} \, \mathrm{d}\Omega + \int_{\Gamma_{\mathrm{d}}} \mathbf{N}^{\mathrm{T}} \mathbf{T}\mathbf{N} \, \mathrm{d}\Gamma \nonumber
\end{align} 
$$

and the arrays at the right-hand side

$$ 
\begin{align} 
\mathbf{f}_{\tilde{a}}^{*} &= \mathbf{f}_{\tilde{a}}^{\mathrm{ext}} - (\mathbf{f}_{\tilde{a}}^{\mathrm{int}})^{i} - \alpha_{0} \mathbf{M}_{\tilde{a}\tilde{a}} \ddot{\tilde{a}}_{n+1}^{i} - \alpha_{0} \mathbf{M}_{\tilde{a}\tilde{a}} \ddot{\tilde{a}}_{n+1}^{i} \nonumber \\ 
\mathbf{f}_{\tilde{a}}^{*} &= \mathbf{f}_{\tilde{a}}^{\mathrm{ext}} - (\mathbf{f}_{\tilde{a}}^{\mathrm{int}})^{i} - \alpha_{0} \mathbf{M}_{\tilde{a}\tilde{a}} \ddot{\tilde{a}}_{n+1}^{i} - \alpha_{0} \mathbf{M}_{\tilde{a}\tilde{a}} \ddot{\tilde{a}}_{n+1}^{i} \nonumber \\ 
\mathbf{f}_{\tilde{p}}^{*} &= \mathbf{f}_{\tilde{p}}^{\mathrm{ext}} - \alpha_{0} \mathbf{M}_{\tilde{p}\tilde{a}} \ddot{\tilde{a}}_{n+1}^{i} - \alpha_{0} \mathbf{M}_{\tilde{p}\tilde{a}} \ddot{\tilde{a}}_{n+1}^{i} - \alpha_{1} \mathbf{M}_{\tilde{p}\tilde{p}} \dot{\tilde{p}}_{n+1}^{i} - \alpha_{1} \mathbf{M}_{\tilde{p}\tilde{p}} \dot{\tilde{p}}_{n+1}^{i} \nonumber \\ 
&\quad - \alpha_{1} \mathbf{K}_{\tilde{a}\tilde{p}}^{\mathrm{T}} \dot{\tilde{a}}_{n+1}^{i} - \alpha_{1} \mathbf{K}_{\tilde{a}\tilde{p}}^{\mathrm{T}} \dot{\tilde{a}}_{n+1}^{i} \nonumber \\ 
\mathbf{f}_{\tilde{p}}^{*} &= \mathbf{f}_{\tilde{p}}^{\mathrm{ext}} - \alpha_{0} \mathbf{M}_{\tilde{p}\tilde{a}} \ddot{\tilde{a}}_{n+1}^{i} - \alpha_{0} \mathbf{M}_{\tilde{p}\tilde{a}} \ddot{\tilde{a}}_{n+1}^{i} - \alpha_{1} \mathbf{M}_{\tilde{p}\tilde{p}} \dot{\tilde{p}}_{n+1}^{i} - \alpha_{1} \mathbf{M}_{\tilde{p}\tilde{p}} \dot{\tilde{p}}_{n+1}^{i} \nonumber \\ 
&\quad - \alpha_{1} \mathbf{K}_{\tilde{a}\tilde{p}}^{\mathrm{T}} \dot{\tilde{a}}_{n+1}^{i} - \alpha_{1} \mathbf{K}_{\tilde{a}\tilde{p}}^{\mathrm{T}} \dot{\tilde{a}}_{n+1}^{i} \tag{53}
\end{align} 
$$

where the superscript signifies that the corresponding quantity has to be evaluated at iteration i. The quantities $\ddot{\tilde{a}}_{n+1}^{i}$, $\dot{\tilde{a}}_{n+1}^{i}$, $\dot{\tilde{a}}_{n+1}^{i}$, $\ddot{\tilde{p}}_{n+1}^{i}$, $\dot{\tilde{p}}_{n+1}^{i}$ are evaluated using Equations (49).

The stiffness matrix of Equation (51) is not symmetric. Symmetry can be restored by multiplying the third and the fourth rows 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 $\mathbf{M}_{\tilde{p}\tilde{a}}$, $\mathbf{M}_{\tilde{p}\tilde{a}}$, $\mathbf{M}_{\tilde{p}\tilde{a}}$ and $\mathbf{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 is affected.

### 5.4. 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 [22], we use a smoothing of the stresses around the tip and compute the stress at the tip by the following non-local-like procedure:

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

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 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.

## 6. EXAMPLE CALCULATIONS

All results of the computations are based on 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 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$/N s. The permeability of the diaphragm was assigned a value $k_d = 0.5 \times 10^{-14}$ m$^2$/N s. Shear-band formation was triggered by a small imperfection, see Figure 4.

A structured mesh that consists of 5841 four-noded elements with equal (bilinear) interpolations has been used for the displacements and the pressure. The simulation is started

<figure><img src="image_5.png" /><figcaption>Figure 4. Geometry and boundary conditions.</figcaption></figure>
using a time step of 0.4 s, which is small enough to accurately follow the pressure evolution and the near-quasi-static behaviour 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 reference 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². The dynamic seepage term has not been taken into account in the reference 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 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 utilize 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.

## 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 frictional softening law is $\mu_c = 0.2 \tan \phi$. A parametric study has been carried out with respect to the influence of $\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 µs.

A first simulation has been carried obtained for a fracture energy $G_c^{II} = 15$ J/m². In this case $\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 value of $\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.
<figure><img src="image_6.png" /><figcaption>Figure 5. Evolution of the pressure field for the Tresca criterion.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. Pressure field near the process zone and tractions at the discontinuity. The scale of the pressures is the same as in Figure 5. The magnitudes of the tractions are proportional to the lengths of the bars.</figcaption></figure>

Increasing the value of $\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 localization inside the shear band and the softening of the interface at the end of propagation. As a consequence, the pressure evolution in Figure 8 has no local cavitation zone. Indeed, cavitation is first obtained around the initiation locus and, subsequently, in the entire specimen. Instability patterns are observed in the last two figures that plot the pressure distribution.

Such instabilities, which are due to friction, are observed at a larger scale when $\delta_c/\Delta U_{ini}$ is increased further up to 0.24 (Figure 9). 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 localization is obtained, but only frictional instabilities are observed. No cavitation is induced because of the absence of strain localization.

The differences between the results using the Tresca-like criterion and those using the Coulomb initiation criterion are also illustrated in Figure 11, which compares the global load-displacement responses of the specimen for both criteria. In either case, three parts in the curve can be distinguished. First, the load remains constant. This corresponds to the wave propagation from the initial defect to the top of the specimen. In the second part, the load decreases because of the propagation of the shear band. As in the Coulomb interface, there is a frictional contribution in the load-carrying capacity, the load decrease is less pronounced as for the Tresca-like criterion. In the third part of the load-displacement curve, the global strength of the specimen is more and more due to transmission of tractions across the shear band. For the Tresca-like criterion, the load now decreases rapidly until the process zone is completely damaged. For the Coulomb criterion, the rate at which the load decreases depends on $\delta_c/\Delta U_{ini}$. When increasing $\delta_c/\Delta U_{ini}$, the rate becomes lower because the interface is less damaged and the frictional contribution remains active along the entire shear band.
<figure><img src="image_8.png" /><figcaption>Figure 7. Evolution of the pressure field for the Coulomb criterion with $\delta_c/\Delta U_{ini} = 0.024$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8. Evolution of the pressure field for the Coulomb criterion with $\delta_c/\Delta U_{ini} = 0.16$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9. Evolution of the pressure field for the Coulomb criterion with $\delta_c/\Delta U_{ini} = 0.24$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10. Shear traction distribution at the end of the shear-band propagation for different values of $\delta_c/\Delta U_{ini}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Load-displacement curves after initiation: comparison between Tresca and Coulomb criteria.</figcaption></figure>

### 6.3. Propagation of shear-band tip

As a consequence of the time discretization and of the rather homogeneous stress field in the specimen, the initiation criterion can be violated over a region that extends over several
elements. Hence, the nucleation criterion is checked in all elements ahead of the actual tip, until the initiation criterion is no longer violated. Figure 12 illustrates this procedure. Figure 13 shows the results for both initiation criteria in terms of the computed shear-band length as a function of time.

For both criteria, the length that has been plotted is that of the topological shear band. Indeed, a topological tip is defined by the enhanced basis terms, but it is not the actual position

<figure><img src="image_13.png" /><figcaption>Figure 12. Illustration of the shear-band tip propagation procedure. The stress tip is estimated at the square and the tip is propagated as long as the criterion is exceeded.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. Shear-band evolution for Tresca and Coulomb criteria.</figcaption></figure>
of the shear-band tip. For the Tresca-like criterion, the actual tip is located at the front of the cohesive zone. Taking this into account, the average speed of the shear-band propagation can be estimated to be in the same order as the dilatational wave speed. This value seems rather high, and a proper explanation for this result has not yet been found. For the Coulomb model, the determination of the actual tip is complicated. For this reason, only an estimate has been done of the propagation speed of the topological tip, which is around 80% of the dilatational wave speed.

### 6.4. Discontinuity in the pressure field

One of the constitutive assumptions in the interface model concerns the introduction of a discontinuity in the pressure field. Since this assumption is less intuitive than the discontinuity in the displacement field, a computation has also been carried out without this assumption. Figure 14 compares the pressure distributions around the process zone obtained in the reference simulation and in a simulation with a continuous pressure field. In the latter case, strong pressure gradients are observed in the elements that are crossed by the discontinuity. Moreover, the cavitation domain is translated to just below these elements, whereas it starts at the interface in the case of a discontinuous pressure field. Hence, the assumption of a discontinuity in the pressure field at the interface seems physically justified as the underlying process tends to create it anyway. Moreover, the computation with a pressure discontinuity seems more robust and permits the definition of different constitutive relations for the fluid flow inside the interface and in the bulk. The latter advantage is not to be underestimated, since the microstructure of the material inside a shear band is usually modified extensively during strain localization, resulting in a marked variation of material properties such as the permeability.

### 6.5. Influence of the relative acceleration

In the beginning of this article, a framework for fluid flow in deforming porous media was set up starting from the balance equations of the two constituent phases, thus including separate inertia terms for both. However, a model with two separate inertia terms requires the independent interpolation of three fields. To simplify the ensuing numerical model, the reduction to one inertia term was subsequently made. Nevertheless, a similar effect as having two separate inertia terms can be obtained by including a so-called dynamic seepage term in

<figure><img src="image_15.png" /><figcaption>Figure 14. Pressure distribution around the process zone for a: (a) discontinuous; and (b) continuous modelling.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15. Influence of the dynamic seepage term (zoom on the load-displacement curves after initiation).</figcaption></figure>

Darcy's law, Equation (28), at the expense of losing symmetry of the system, Equation (51). To assess the effect of including the dynamic seepage term in Darcy's relation, simulations have been carried out in which this effect has been incorporated. Figure 15 presents the beginning of the load-displacement curves for the reference simulation with the Tresca criterion and for a simulation in which the dynamic seepage term has been included. No clear distinction between both cases is possible. This observation also holds for local quantities, such as the pressures, where the differences are also negligible (about $10^{-6}$), and for a simulation in which both inertia terms were taken into account. Very similar results were obtained when including both inertia terms in the analysis.

## 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 discretization. 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 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 localization and strongly depends on the constitutive model for the discontinuity.

## REFERENCES

1. de Borst R. Damage, material instabilities, and failure. In *Encyclopedia of Computational Mechanics*, Stein E, de Borst R, Hughes TJR (eds), vol. 2, Chapter 10. Wiley: Chichester, 2004.
2. Babuska I, Melenk JM. The partition of unity method. *International Journal for Numerical Methods in Engineering* 1997; **40**:727–758.

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

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

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

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

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

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

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

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

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

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

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

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

15. Jouanna P, Abellan MA. Generalized approach to heterogeneous media. In *Modern Issues in Non-Saturated Soils*, Gens A, Jouanna P, Schrefler B (eds). Springer: Wien, New York, 1995; 1–128.

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

17. de Borst R, Réthoré J, Abellan MA. A numerical approach for arbitrary cracks in a fluid-saturated medium. *Archive of Applied Mechanics* 2006; **75**:595–606.

18. Armero F, Callari C. An analysis of strong discontinuities in a saturated poro-plastic solid. *International Journal for Numerical Methods in Engineering* 1999; **46**:1673–1698.

19. Larsson J, Larsson R. Localization analysis of a fluid-saturated elastoplastic porous medium using regularized discontinuities. *Mechanics of Cohesive-frictional Materials* 2000; **5**:565–582.

20. Schrefler, B, Scotta R. A fully coupled dynamic model for two-phase fluid flow in deformable porous media. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**:3223–3246.

21. Feenstra PH, de Borst R. A composite plasticity model for concrete. *International Journal of Solids and Structures* 1996; **33**:707–730.

22. Jirasek M. Embedded crack models for concrete fracture. In *Computational Modelling of Concrete Structures*, de Borst R, Bicanic N, Mang H, Meschke G (eds). Balkema: Rotterdam, 1998; 291–300.