# A numerical approach for arbitrary cracks in a fluid-saturated medium

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

► To cite this version:

René de Borst, Julien Réthoré, Marie-Angèle Abellan. A numerical approach for arbitrary cracks in a fluid-saturated medium. Archive of Applied Mechanics, 2006, 75, pp.595-606. 10.1007/s00419-006-0023-y. hal-00938532

HAL Id: hal-00938532

https://hal.science/hal-00938532

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 Attribution 4.0 International License
René de Borst · Julien Réthoré · Marie-Angèle Abellan

# A numerical approach for arbitrary cracks in a fluid-saturated medium

**Abstract** A finite element method is proposed that can capture arbitrary discontinuities in a two-phase medium. The discontinuity is described in an exact manner 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, while at the discontinuity a discrete analogon of Darcy's relation is proposed. The results of this finite element model are independent of the original discretisation, as is demonstrated by an example of shear banding in a biaxial, plane-strain specimen.

**Keywords** Fracture · Shear bands · Multi-phase medium · Finite elements · Partition-of-unity method

## 1 Introduction

Broadly speaking, two approaches exist to 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 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 perhaps to model discontinuities in a discrete manner, thus reflecting the change in topology that actually takes place in the solid when a discontinuity propagates. In finite element analysis this approach was probably pioneered by Ngo and Scordelis [2]. In their seminal paper, nodes were simply split when the tensile force in this node exceeded a fracture criterion. Evidently, this approach suffers from the constraint that discontinuities can only propagate along the mesh lines, thus

R. de Borst (✉) · J. Réthoré
Delft University of Technology, Faculty of Aerospace Engineering, P.O. Box 5058,
2600 GB Delft, The Netherlands
E-mail: r.deborst@tudelft.nl
Tel.: +31-15-2785464
Fax: +31-15-2611465

M.-A. Abellan
LTDS-ENISE - UMR CNRS 5513, 42043 Saint-Etienne, France
biasing the possible crack directions by the discretisation. Advances in mesh generation and automatic reme-shing have largely allowed to overcome this restriction, e.g. [3,4]. More recently, finite element methods have been constructed that exploit the partition-of-unity property of finite element shape functions [5–15]. While preserving the original discretisation, the addition of extra degrees-of-freedom to elements which are crossed by a discontinuity allows to construct two continuous displacement fields that are separated by a Heaviside function at the discontinuity.

Increasingly, contemporary problems in mechanics 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. Application areas include durability problems, where moisture, ion and heat transfer are crucial, joining technologies such as welding, micro-devices, where thermal and electrical excitations are used to let the device perform certain mechanical functions, and biomaterials, where electrical, chemical, and hydro-mechanical interactions can play a role. It is the purpose of this manuscript to formulate a numerical model that is capable of describ-ing fracture in a porous medium, with a solid skeleton and an interstitial fluid as the constituent phases, in a mesh-independent manner. The model also 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 prevail.

The manuscript starts with a concise derivation of the governing equations for a fluid-saturated porous medium. Subsequently, the general methodology and the assumptions regarding the introduction of a dis-continuity 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. The versatility of the method is illustrated by an analysis of a biaxial, fluid-saturated, plane-strain specimen that is loaded in uniaxial compression.

## 2 Governing equations

We consider a two-phase medium subject to the restriction of small displacement gradients and small varia-tions in the concentrations [16]. 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(\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 = \text{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=\text{s,f}} \hat{\mathbf{p}}_{\pi} = \mathbf{0}. \quad (2) $$

We now neglect convective terms and confine attention to quasi-static processes, so that the momentum bal-ances reduce to:

$$ \nabla \cdot \boldsymbol{\sigma}_{\pi} + \hat{\mathbf{p}}_{\pi} = \mathbf{0}. \quad (3) $$

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

$$ \nabla \cdot \boldsymbol{\sigma} = \mathbf{0}, \quad (4) $$

where

$$ \boldsymbol{\sigma} = \boldsymbol{\sigma}_{\text{s}} + \boldsymbol{\sigma}_{\text{f}} = \boldsymbol{\sigma}_{\text{s}} - \alpha p \mathbf{I} \quad (5) $$
with p the fluid pressure, **I** the second-order identity tensor, and α the Biot coefficient, cf. [17].

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

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

In consideration of the assumption of small changes in the concentrations, the mass balances can be simplified to give:

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

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

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

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

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_f}{\rho_s} \frac{\partial \rho_s}{\partial t} \qquad (10) $$

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$, this equation can be rewritten as

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

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

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

with the compressibility modulus

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

where $K_f$ is the bulk modulus of the fluid. Inserting relations (11) and (12) into the balance of mass of the total medium, Eq. (9), 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 (14) $$

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

$$ \epsilon_s = \nabla^s \mathbf{u}_s \qquad (15) $$

with $\mathbf{u}_s$, $\epsilon_s$ the displacement and strain fields of the solid, respectively, the superscript s denoting the symmetric part of the gradient operator, and an incrementally linear stress–strain relation for the solid skeleton,

$$ d\boldsymbol{\sigma}_s = \mathbf{D}^{\mathrm{tan}} : d\boldsymbol{\epsilon}_s, \qquad (16) $$
where $D^{\text{tan}}$ is the fourth-order tangent stiffness tensor of the solid material and the d symbol denotes an infinitesimally small increment. For the pore fluid flow, Darcy's relation for isotropic media is assumed to hold,

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

with $k_f$ the permeability coefficient of the porous medium. The boundary conditions

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

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) = \mathbf{q}_p, \quad p = p_p \qquad (19)$$

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

### 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 [6,7], exploiting the partition-of-unity property of finite element shape functions [5]. 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{\bar{a}}_{jk} \right) \qquad (20)$$

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 Fig. 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 (21)$$

where $\mathcal{H}_{\Gamma_d}$ is the Heaviside step function centred at the discontinuity:

$$\mathcal{H}_{\Gamma_d} = \begin{cases} 0 & \text{if } \mathbf{x} \in \Omega^{-} \\ 1 & \text{if } \mathbf{x} \in \Omega^{+}. \end{cases} \qquad (22)$$

<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>
The decomposition in Eq. (21) has a structure similar to the interpolation in Eq. (20). This can be seen directly by rewriting and specialising Eq. (20) as:

$$ \mathbf{v} = \mathbf{N}(\dot{\mathbf{a}} + \mathcal{H}_{\Gamma_d} \dot{\mathbf{a}}) = \mathbf{N}\dot{\mathbf{a}} + \mathcal{H}_{\Gamma_d} \mathbf{N}\dot{\mathbf{a}} = \bar{\mathbf{v}} + \mathcal{H}_{\Gamma_d} \tilde{\mathbf{v}}, \quad (23) $$

where N contains the standard shape functions $\varphi_j$, and $\bar{\mathbf{a}}$ and $\tilde{\mathbf{a}}$ collect the conventional and the additional displacement nodal degrees-of-freedom, respectively. Accordingly, the partition-of-unity property of finite element shape functions can be used in a straightforward fashion to incorporate discontinuities, and thus discrete crack models, in a manner that preserves the discontinuous character of cracks.

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. (21), we have:

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

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. (24), we obtain:

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

## 4 Constitutive model at the discontinuity

The above kinematics and the balance equations for momentum and mass have to be complemented by constitutive equations. For the bulk material the governing equations for the solid skeleton and for the fluid flow have been outlined in one of the preceding sections. In the example, a linear-elastic constitutive relation for the skeleton has been adopted, so that $D^{\text{tan}}$ in Eq. (16) reduces to the linear-elastic stiffness tensor $D$.

Using a maximum shear stress criterion in the spirit of Tresca, $\tau = |\sigma_2 - \sigma_1|/2$, with $\sigma_1, \sigma_2$ the principal stresses, a shear-band like discontinuity is created when the criterion exceeds the critical value $\tau_c$. The orientation of the interface is chosen such that it maximises the shear stress. 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 have been added to the finite element model that describe this sliding mode, which is different from earlier shear-band simulations (for single-phase media) that have exploited the partition-of-unity property of finite element shape functions [8, 9, 12, 15].

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

$$ \mathbf{t}_d = \mathbf{t}_d(\mathbf{\delta}, \kappa) \quad (26) $$

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

$$ \dot{\mathbf{t}}_d = \mathbf{T}\dot{\mathbf{\delta}} \quad (27) $$

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

$$ \mathbf{T} = \frac{\partial \mathbf{t}_d}{\partial \mathbf{\delta}} + \frac{\partial \mathbf{t}_d}{\partial \kappa} \frac{\partial \kappa}{\partial \mathbf{\delta}}. \quad (28) $$

A key element is the presence of a mode-II fracture energy, $G_c^{\text{II}}$, which governs the shear band evolution and enters the interface constitutive relation (26) 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. [20]:

$$ G_c^{\text{II}} = \int_{\delta_t=0}^{\infty} \tau d\delta_t \quad (29) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2 Relation between relative sliding at the discontinuity and shear tractions</figcaption></figure>

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

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_{x \in \Gamma_d}, \qquad (30) $$

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. (30). 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. Another possible boundary condition for the fluid flow at an internal boundary $\Gamma_d$ is

$$ \mathbf{n}_{\Gamma_d} \cdot \mathbf{q}_d = q_d \mid_{x \in \Gamma_d}. \qquad (31) $$

Such an internal boundary condition represents the existence of a drain (or a line source) with a capacity $q_d$ per unit length.

## 5 Weak form and numerical elaboration

To arrive at the weak form of the balance equations, we multiply the momentum balance (4) and the mass balance (14) 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 = \bar{\eta} + \mathcal{H}_{\Gamma_d} \tilde{\eta} \qquad (32) $$

for the velocities, and

$$ \zeta = \bar{\zeta} + \mathcal{H}_{\Gamma_d} \tilde{\zeta} \qquad (33) $$

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

$$ \int_{\Omega} (\bar{\eta} + \mathcal{H}_{\Gamma_d} \tilde{\eta}) \cdot \nabla \cdot \boldsymbol{\sigma} \, d\Omega = 0 \qquad (34) $$
and

$$ \int_{\Omega} (\bar{\xi} + \mathcal{H}_{\Gamma_d} \tilde{\xi}) \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) \mathrm{d}\Omega = 0. \quad (35) $$

Using the standard procedure of applying the divergence theorem, using the external boundary conditions (18) and (19), 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:

$$ \int_{\Omega} (\nabla \cdot \bar{\eta}) \cdot \boldsymbol{\sigma} \mathrm{d}\Omega + \int_{\Omega^+} (\nabla \cdot \tilde{\eta}) \cdot \boldsymbol{\sigma} \mathrm{d}\Omega + \int_{\Gamma_d} \tilde{\eta} \cdot \mathbf{t}_d \mathrm{d}\Omega = \int_{\Gamma} (\bar{\eta} + \mathcal{H}_{\Gamma_d} \tilde{\eta}) \cdot \mathbf{t}_p \mathrm{d}\Omega \quad (36) $$

and

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

where for the latter equation also Darcy's relation (17) has been employed. We now switch to matrix-vector notation and discretise the trial functions $v_s$ and $p$ and the test functions $\eta$ and $\zeta$ as:

$$ \mathbf{v}_s = \mathbf{N}(\dot{\bar{\mathbf{a}}} + \mathcal{H}_{\Gamma_d} \dot{\tilde{\mathbf{a}}}), \quad (38) $$

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

$$ \boldsymbol{\eta} = \mathbf{N}(\bar{\mathbf{w}} + \mathcal{H}_{\Gamma_d} \tilde{\mathbf{w}}), \quad (40) $$

$$ \zeta = \mathbf{H}(\bar{\mathbf{z}} + \mathcal{H}_{\Gamma_d} \tilde{\mathbf{z}}). \quad (41) $$

Inserting Eqs. (40) and (41) into Eqs. (36) and (37) and requiring that the result holds for all admissible $\bar{\mathbf{w}}, \bar{\mathbf{z}}, \tilde{\mathbf{w}}$ and $\tilde{\mathbf{z}}$ gives:

$$ \int_{\Omega} \nabla \mathbf{N}^{\mathrm{T}} \boldsymbol{\sigma} \mathrm{d}\Omega = \int_{\Gamma} \mathbf{N}^{\mathrm{T}} \mathbf{t}_p \mathrm{d}\Gamma \quad (42) $$

$$ \int_{\Omega^+} \nabla \mathbf{N}^{\mathrm{T}} \boldsymbol{\sigma} \mathrm{d}\Omega + \int_{\Gamma_d} \mathbf{N}^{\mathrm{T}} \mathbf{t}_d \mathrm{d}\Gamma = \int_{\Gamma} \mathcal{H}_{\Gamma_d} \mathbf{N}^{\mathrm{T}} \mathbf{t}_p \mathrm{d}\Gamma \quad (43) $$

$$ - \int_{\Omega} \alpha \mathbf{H}^{\mathrm{T}} \mathbf{m}^{\mathrm{T}} \nabla \mathbf{v}_s \mathrm{d}\Omega + \int_{\Omega} k_f \nabla \mathbf{H}^{\mathrm{T}} \nabla p \mathrm{d}\Omega - \int_{\Omega} Q^{-1} \mathbf{H}^{\mathrm{T}} \dot{p} \mathrm{d}\Omega = \int_{\Gamma} \mathbf{H}^{\mathrm{T}} \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \mathbf{q}_p \mathrm{d}\Gamma \quad (44) $$

$$ \begin{aligned} & - \int_{\Omega^+} \alpha \mathbf{H}^{\mathrm{T}} \mathbf{m}^{\mathrm{T}} \nabla \mathbf{v}_s \mathrm{d}\Omega + \int_{\Omega^+} k_f \nabla \mathbf{H}^{\mathrm{T}} \nabla p \mathrm{d}\Omega - \int_{\Omega^+} Q^{-1} \mathbf{H}^{\mathrm{T}} \dot{p} \mathrm{d}\Omega + \int_{\Gamma_d} \mathbf{H}^{\mathrm{T}} \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \mathbf{q}_d \mathrm{d}\Gamma \\ & = \int_{\Gamma} \mathcal{H}_{\Gamma_d} \mathbf{H}^{\mathrm{T}} \mathbf{n}_{\Gamma_d}^{\mathrm{T}} \mathbf{q}_p \mathrm{d}\Gamma \end{aligned} \quad (45) $$

with, for two dimensions, $\mathbf{m} = [1, 1, 0]$
For use in a Newton-Raphson solution method, Eqs. (42)–(45) need to be linearised. To this end, the stress and the pressure are decomposed as

$$ \sigma_j = \sigma_{j-1} + d\sigma \qquad (46) $$

and

$$ p_j = p_{j-1} + dp \qquad (47) $$

with the subscripts $j-1$ and $j$ signifying the iteration numbers. Substituting these decompositions into the discrete set of equations (42)–(45), utilising Eq. (5), the kinematic relation (15), the stress–strain relation (16) for the bulk material and the traction–relative displacement relation (26) at the interface, using Darcy’s relation for the fluid flow in the porous medium (17), its discrete equivalent (30) at the interface and the expression (25) for the pressure gradient, and using the interpolations for the displacement and the pressure according to Eqs. (38) and (39), leads to the following set of equations linearised at iteration $j-1$:

$$ \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{K}_{\tilde{a}\tilde{p}}^{\mathrm{T}} & \mathbf{K}_{\tilde{a}\tilde{p}}^{\mathrm{T}} & \mathbf{K}_{\tilde{p}\tilde{p}}^{(1)} & \mathbf{K}_{\tilde{p}\tilde{p}}^{(1)} \\ \mathbf{K}_{\tilde{a}\tilde{p}}^{\mathrm{T}} & \mathbf{K}_{\tilde{a}\tilde{p}}^{\mathrm{T}} & \mathbf{K}_{\tilde{p}\tilde{p}}^{(1)} & \mathbf{K}_{\tilde{p}\tilde{p}}^{(1)} \end{bmatrix} \begin{pmatrix} \dot{\tilde{a}} \\ \dot{\tilde{a}} \\ \dot{\tilde{p}} \\ \dot{\tilde{p}} \end{pmatrix} + \begin{bmatrix} \mathbf{K}_{\tilde{a}\tilde{a}} & \mathbf{K}_{\tilde{a}\tilde{a}} & \mathbf{K}_{\tilde{a}\tilde{p}} & \mathbf{K}_{\tilde{a}\tilde{p}} \\ \mathbf{K}_{\tilde{a}\tilde{a}} & \mathbf{K}_{\tilde{a}\tilde{a}} & \mathbf{K}_{\tilde{a}\tilde{p}} & \mathbf{K}_{\tilde{a}\tilde{p}} \\ \mathbf{0} & \mathbf{0} & \mathbf{K}_{\tilde{p}\tilde{p}}^{(2)} & \mathbf{K}_{\tilde{p}\tilde{p}}^{(2)} \\ \mathbf{0} & \mathbf{0} & \mathbf{K}_{\tilde{p}\tilde{p}}^{(2)} & \mathbf{K}_{\tilde{p}\tilde{p}}^{(2)} \end{bmatrix} \begin{pmatrix} d\tilde{a}_j \\ d\tilde{a}_j \\ d\tilde{p}_j \\ d\tilde{p}_j \end{pmatrix} = \begin{pmatrix} \mathbf{f}_{\tilde{a}}^{\mathrm{ext}} & -\mathbf{f}_{\tilde{a},j-1}^{\mathrm{int}} \\ \mathbf{f}_{\tilde{a}}^{\mathrm{ext}} & -\mathbf{f}_{\tilde{a},j-1}^{\mathrm{int}} \\ \mathbf{f}_{\tilde{p}}^{\mathrm{ext}} & -\mathbf{f}_{\tilde{p},j-1}^{\mathrm{int}} \\ \mathbf{f}_{\tilde{p}}^{\mathrm{ext}} & -\mathbf{f}_{\tilde{p},j-1}^{\mathrm{int}} \end{pmatrix} \quad (48) $$

with $\mathbf{f}_{\tilde{a}}^{\mathrm{ext}} \cdots \mathbf{f}_{\tilde{p}}^{\mathrm{ext}}$ given by the right-hand sides of Eqs. (42)–(45), respectively, and with $\mathbf{f}_{\tilde{a}}^{\mathrm{int}} \cdots \mathbf{f}_{\tilde{p}}^{\mathrm{int}}$ given by the left-hand sides of Eqs (42)–(45). The stiffness matrices are defined as:

$$ \mathbf{K}_{\tilde{a}\tilde{a}} = \int_{\Omega} \mathbf{B}^{\mathrm{T}} \mathbf{DB} d\Omega, \qquad (49) $$

$$ \mathbf{K}_{\tilde{a}\tilde{a}} = \int_{\Omega^{+}} \mathbf{B}^{\mathrm{T}} \mathbf{DB} d\Omega, \qquad (50) $$

$$ \mathbf{K}_{\tilde{a}\tilde{a}} = \int_{\Omega^{+}} \mathbf{B}^{\mathrm{T}} \mathbf{DB} d\Omega + \int_{\Gamma_d} \mathbf{N}^{\mathrm{T}} \mathbf{T} \mathbf{N} d\Gamma, \qquad (51) $$

$$ \mathbf{K}_{\tilde{a}\tilde{p}} = - \int_{\Omega} \alpha \mathbf{B}^{\mathrm{T}} \mathbf{m} \mathbf{H} d\Omega, \qquad (52) $$

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

$$ \mathbf{K}_{\tilde{p}\tilde{p}}^{(1)} = - \int_{\Omega} Q^{-1} \mathbf{H}^{\mathrm{T}} \mathbf{H} d\Omega, \qquad (54) $$

$$ \mathbf{K}_{\tilde{p}\tilde{p}}^{(1)} = - \int_{\Omega^{+}} Q^{-1} \mathbf{H}^{\mathrm{T}} \mathbf{H} d\Omega, \qquad (55) $$

$$ \mathbf{K}_{\tilde{p}\tilde{p}}^{(2)} = - \int_{\Omega} k_{\mathrm{f}} \nabla \mathbf{H}^{\mathrm{T}} \nabla \mathbf{H} d\Omega, \qquad (56) $$

$$ \mathbf{K}_{\tilde{p}\tilde{p}}^{(2)} = - \int_{\Omega^{+}} k_{\mathrm{f}} \nabla \mathbf{H}^{\mathrm{T}} \nabla \mathbf{H} d\Omega, \qquad (57) $$

$$ \mathbf{K}_{\tilde{p}\tilde{p}}^{(2)} = - \int_{\Omega^{+}} k_{\mathrm{f}} \nabla \mathbf{H}^{\mathrm{T}} \nabla \mathbf{H} d\Omega - \int_{\Gamma_d} k_{\mathrm{d}} \mathbf{H}^{\mathrm{T}} \mathbf{H} d\Gamma. $$

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

$$ \left(\frac{d(\cdot)}{dt}\right)^{t+\Delta t} = \frac{(\cdot)^{t+\Delta t} - (\cdot)^t}{\Delta t}, \qquad (58) $$
where $\Delta t$ is the time increment, while $(\cdot)^t$ and $(\cdot)^{t+\Delta t}$ denote the unknowns at $t$ and $t + \Delta t$, respectively. Furthermore, a decomposition of the arrays that contain the nodal displacements and nodal pressures is adopted similar to that of the stresses and the pressures, cf. Eqs.(46) and (47), $(\cdot)_j = (\cdot)_{j-1} + d(\cdot)$. Combining this with Eq. (58) yields

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

Substituting this identity into Eq.(48) for all arrays finally gives

$$ \begin{bmatrix} \mathbf{K}_{\bar{a}\bar{a}} & \mathbf{K}_{\bar{a}\bar{a}} & \mathbf{K}_{\bar{a}\bar{p}} & \mathbf{K}_{\bar{a}\bar{p}} \\ \mathbf{K}_{\bar{a}\bar{a}}^{\mathrm{T}} & \mathbf{K}_{\bar{a}\bar{a}} & \mathbf{K}_{\bar{a}\bar{p}} & \mathbf{K}_{\bar{a}\bar{p}} \\ \mathbf{K}_{\bar{a}\bar{p}}^{\mathrm{T}} & \mathbf{K}_{\bar{a}\bar{p}}^{\mathrm{T}} & \mathbf{K}_{\bar{p}\bar{p}}^{(1)} + \Delta t \mathbf{K}_{\bar{p}\bar{p}}^{(2)} & \mathbf{K}_{\bar{p}\bar{p}}^{(1)} + \Delta t \mathbf{K}_{\bar{p}\bar{p}}^{(2)} \\ \mathbf{K}_{\bar{a}\bar{p}}^{\mathrm{T}} & \mathbf{K}_{\bar{a}\bar{p}}^{\mathrm{T}} & \mathbf{K}_{\bar{p}\bar{p}}^{(1)} + \Delta t \mathbf{K}_{\bar{p}\bar{p}}^{(2)} & \mathbf{K}_{\bar{p}\bar{p}}^{(1)} + \Delta t \mathbf{K}_{\bar{p}\bar{p}}^{(2)} \end{bmatrix} \begin{pmatrix} d\bar{a}_j \\ d\tilde{a}_j \\ d\bar{p}_j \\ d\tilde{p}_j \end{pmatrix} = \begin{pmatrix} \mathbf{f}_{\bar{a}}^{\mathrm{ext}} & - \mathbf{f}_{\bar{a},j-1}^{\mathrm{int}} \\ \mathbf{f}_{\bar{a}}^{\mathrm{ext}} & - \mathbf{f}_{\bar{a},j-1}^{\mathrm{int}} \\ \Delta t \mathbf{f}_{\bar{p}}^{\mathrm{ext}} & - \mathbf{f}_{\bar{p},j-1}^{\mathrm{int},t+\Delta t} \\ \Delta t \mathbf{f}_{\bar{p}}^{\mathrm{ext}} & - \mathbf{f}_{\bar{p},j-1}^{\mathrm{int},t+\Delta t} \end{pmatrix} \quad (60) $$

with

$$ \mathbf{f}_{\bar{p},j-1}^{\mathrm{int},t+\Delta t} = \mathbf{K}_{\bar{a}\bar{p}}^{\mathrm{T}} (\bar{a}_{j-1}^{t+\Delta t} - \bar{a}^t) + \mathbf{K}_{\bar{a}\bar{p}}^{\mathrm{T}} (\tilde{a}_{j-1}^{t+\Delta t} - \tilde{a}^t) + \mathbf{K}_{\bar{p}\bar{p}}^{(1)} (\bar{p}_{j-1}^{t+\Delta t} - \bar{p}^t) \\ + \mathbf{K}_{\bar{p}\bar{p}}^{(1)} (\tilde{p}_{j-1}^{t+\Delta t} - \tilde{p}^t) + \Delta t \mathbf{K}_{\bar{p}\bar{p}}^{(2)} \bar{p}_{j-1}^{t+\Delta t} + \Delta t \mathbf{K}_{\bar{p}\bar{p}}^{(2)} \tilde{p}_{j-1}^{t+\Delta t}, \quad (61) $$

$$ \mathbf{f}_{\bar{p},j-1}^{\mathrm{int},t+\Delta t} = \mathbf{K}_{\bar{a}\bar{p}}^{\mathrm{T}} (\bar{a}_{j-1}^{t+\Delta t} - \bar{a}^t) + \mathbf{K}_{\bar{a}\bar{p}}^{\mathrm{T}} (\tilde{a}_{j-1}^{t+\Delta t} - \tilde{a}^t) + \mathbf{K}_{\bar{p}\bar{p}}^{(1)} (\bar{p}_{j-1}^{t+\Delta t} - \bar{p}^t) \\ + \mathbf{K}_{\bar{p}\bar{p}}^{(1)} (\tilde{p}_{j-1}^{t+\Delta t} - \tilde{p}^t) + \Delta t \mathbf{K}_{\bar{p}\bar{p}}^{(2)} \bar{p}_{j-1}^{t+\Delta t} + \Delta t \mathbf{K}_{\bar{p}\bar{p}}^{(2)} \tilde{p}_{j-1}^{t+\Delta t}. \quad (62) $$

<figure><img src="image_4.png" /><figcaption>Fig. 3 Geometry and boundary conditions of two-dimensional example problem</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Evolution of the displacement field at $t = 6, 10, 14, 18, 22$ and $26$ s</figcaption></figure>

## 6 Example: shear band formation

We consider a two-dimensional specimen with a width $w = 0.04$ m and a height $H = 0.1$ m, Figure 3, 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^{-4}$ m/s. The pore pressure at the top of the specimen is the reference pressure, here zero, and undrained boundary conditions have been imposed on the other boundaries. 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 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. Shear-band formation was triggered by a small imperfection, see Fig. 3. The critical shear stress at which nucleation occurred, was taken as $\tau_c = 100$ MPa and after inception the shear-band evolution is controlled by $G_c^{II} = 500$ J/m$^2$. In the example calculations, the permeability of the diaphragm has been assigned a value, $k_d = 0.5 \times 10^{-14}$ m$^3$/Ns, which is half of that in the bulk material.

The specimen has been discretised with quadrilateral basis elements with a bilinear interpolation for the displacements as well as for the pressure. It is recognised that this scheme is not optimal, but the added complexity which an unequal interpolation would cause in the presence of a propagating discontinuity, and especially with respect to the integration of the different parts of the load vectors and stiffness matrices in elements which are crossed by the discontinuity, has led to a decision which favours simplicity over accuracy.
<figure><img src="image_6.png" /><figcaption>Fig. 5 Evolution of the pressure field at $t = 6, 10, 14, 18, 22$ and $26$ s</figcaption></figure>

In the simulation of Figs. 4 and 5, 24 elements have been used across the width of the specimen and 60 elements have been used over the height. The duration of the simulation, which has been carried out in 65 equally sized time steps, is 26 s.

Figures 4 and 5 show the evolution of the displacement field and the pressure field, respectively. At the end of the simulation, the discontinuity has traversed the entire specimen. At this point, the bulk has unloaded elastically and the external loading is transmitted entirely through the discontinuity, at which the deformations have localised. The global inclination of the interface is approximately $45^\circ$, which is in conformity with the behaviour in a single phase medium.

The pressure distribution is strongly influenced by the propagation of the interface even if $k_f$ and $k_d$ are of the same order of magnitude. Indeed, the pressure discontinuity is significant as observed in Fig. 5. Accordingly, the relatively lower permeability at the discontinuity has a major influence on the fluid flow, and therefore on the stress distribution inside the specimen.

## 7 Discussion and conclusions

A finite element method that exploits the partition-of-unity property of finite element shape functions has been used to analyse the stress evolution and fluid flow in a biaxial, plane-strain specimen with a propagating
discontinuity, e.g. a crack or a shear band. The approach appears to be accurate and powerful, and enables
analyses that are free from the underlying discretisation. The present formulation is restricted to fluid-saturated
(two-phase) media, but the methodology can be extended to a host of coupled problems in a rather straight-
forward manner.

Nevertheless, some limitations still exist. One of them, which specifically relates to the case studied, con-
cerns the velocity of the tip of the discontinuity. As the underlying two-phase medium appears to have a length
scale that vanishes in the short-wave limit [21], the tip velocity remains undetermined for quasi-static loading
conditions, and, in the simulations, is set by the space and time discretisation vtip = h/Δt = 0.004 m/s, where
h is the spacing of the discretisation and Δt is the time step. It is believed that the inclusion of inertia terms
in the formulation can remedy this deficiency. This will be the discussed in greater detail in a forthcoming
publication.

References

1. de Borst, R.:Damage, material instabilities, and failure. In: Stein, E., de Borst, R., Hughes, T.J.R. (eds.) Encyclopedia of computational mechanics, Vol. 2, Chapter 10. Wiley, Chichester (2004)

2. Ngo, D., Scordelis, A.C.: Finite element analysis of reinforced concrete beams. J Am Concrete Inst **64**, 152–163 (1967)

3. Ingraffea, A.R., Saouma, V.: Numerical modelling of discrete crack propagation in reinforced and plain concrete. In: Sih, G.C., de Tomaso, A. (eds.) Fracture Mechanics of Concrete, pp. 171–225. Martinus Nijhoff Publishers, Dordrecht (1985)

4. Camacho, G.T., Ortiz, M.: Computational modeling of impact damage in brittle materials. Int J Solids Struct **33**, 2899–2938 (1996)

5. Babuska, I., Melenk, J.M.: The partition of unity method. Int J Num Meth Eng **40**, 727–758 (1997)

6. Belytschko, T., Black, T.: Elastic crack growth in finite elements with minimal remeshing. Int J Num Meth Eng **45**, 601–620 (1999)

7. Moës, N., Dolbow, J., Belytschko, T.: A finite element method for crack growth without remeshing. Int J Num Meth Eng **46**, 131–150 (1999)

8. Wells, G.N., Sluys, L.J.: Discontinuous analysis of softening solids under impact loading. Int J Num Anal Meth Geomech **25**, 691–709 (2001)

9. Wells, G.N., Sluys, L.J., de Borst, R.: Simulating the propagation of displacement discontinuities in a regularized strain-softening medium. Int J Num Meth Eng **53**, 1235–1256 (2002)

10. Wells, G.N., de Borst, R., Sluys, L.J.: A consistent geometrically non-linear approach for delamination. Int J Num Meth Eng **54**, 1333–1355 (2002)

11. Remmers, J.J.C., de Borst, R., Needleman, A.: A cohesive segments method for the simulation of crack growth. Comput Mech **31**, 69–77 (2003)

12. Samaniego, E., Belytschko, T.: Continuum–discontinuum modelling of shear bands. Int J Num Meth Eng **62**, 1857–1872 (2005)

13. Réthoré, J., Gravouil, A., Combescure, A.: An energy-conserving scheme for dynamic crack growth using the extended finite element method. Int J Num Meth Eng **63**, 631–659 (2005)

14. Réthoré, J., Gravouil, A., Combescure, A.: A combined space–time extended finite element method. Int J Num Meth Eng **64**, 260–284 (2005)

15. Areias, P.M.A., Belytschko, T.: Twoscale shear band evolution by local partition of unity. Int J Num Meth Eng doi: 10.1002/nme.1589 (2006)

16. Jouanna, P., Abellan, M.A.: Generalized approach to heterogeneous media. In: Gens, A., Jouanna, P., Schrefler, B. (eds.) Modern Issues in non-saturated soils, pp. 1–128, Springer-Verlag, Wien, New York (1995)

17. Lewis, R.W., Schrefler, B.A.: The finite element method in the static and dynamic deformation and consolidation of porous media, 2nd Edn. Chichester: John Wiley & Sons (1998)

18. Armero, F., Callari, C.: An analysis of strong discontinuities in a saturated poro-plastic solid. Int J Num Meth Eng **46**, 1673–1698 (1999)

19. Larsson, J., Larsson, R.: Localization analysis of a fluidsaturated elastoplastic porous medium using regularized discontinuities. Mech Coh-frict Mat **5**, 565–582 (2000)

20. Feenstra, P.H., de Borst, R.: A composite plasticity model for concrete. Int J Solids Struct **33**, 707–730 (1996)

21. Abellan, M.A., de Borst, R.: Wave propagation and localisation in a softening two-phase medium. Comp Meth Appl Mech Eng doi: 10.1016/j.cma.2005.05.056 (2006)