This is a repository copy of On viscoplastic regularisation of strain‐softening rocks and
soils.
White Rose Research Online URL for this paper:
https://eprints.whiterose.ac.uk/156517/
Version: Accepted Version
Article:
Borst, R. and Duretz, T. (2020) On viscoplastic regularisation of strain‐softening rocks and
soils. International Journal for Numerical and Analytical Methods in Geomechanics, 44 (6).
pp. 890-903. ISSN 0363-9061
https://doi.org/10.1002/nag.3046
This is the peer reviewed version of the following article: de Borst, R, Duretz, T. On
viscoplastic regularisation of strain-softening rocks and soils. Int J Numer Anal Methods
Geomech. 2020; 44: 890– 903, which has been published in final form at
https://doi.org/10.1002/nag.3046. This article may be used for non-commercial purposes in
accordance with Wiley Terms and Conditions for Use of Self-Archived Versions. This
article may not be enhanced, enriched or otherwise transformed into a derivative work,
without express permission from Wiley or by statutory rights under applicable legislation.
Copyright notices must not be removed, obscured or modified. The article must be linked
to Wiley’s version of record on Wiley Online Library and any embedding, framing or
otherwise making available the article or pages thereof by third parties from platforms,
services and websites other than Wiley Online Library must be prohibited.
eprints@whiterose.ac.uk
https://eprints.whiterose.ac.uk/
Reuse
Items deposited in White Rose Research Online are protected by copyright, with all rights reserved unless
indicated otherwise. They may be downloaded and/or printed for private study, or other acts as permitted by
national copyright laws. The publisher or other rights holders may allow further reproduction and re-use of
the full text version. This is indicated by the licence information on the White Rose Research Online record
for the item.
Takedown
If you consider content in White Rose Research Online to be in breach of UK law, please notify us by
emailing eprints@whiterose.ac.uk including the URL of the record and the reason for the withdrawal request.
Received: Added at production Revised: Added at production Accepted: Added at production
DOI: xxx/xxxx
ARTICLE TYPE
On viscoplastic regularisation of strain-softening rocks and soils
René de Borst*
1
| Thibault Duretz
2
1
Department of Civil and Structural
Engineering, University of Sheffield,
Sheffield, UK
2
Univ Rennes, CNRS, Géosciences Rennes
UMR 6118, F-35000 Rennes, France
Correspondence
*René de Borst, Department of Civil and
Structural Engineering, University of
Sheffield, Sheffield S1 3JD, UK. Email:
r.deborst@sheffield.ac.uk
Funding information
Horizon 2020 European Research Council
Grant 664734 "PoroFrac"
Summary
Constitutive models for rocks and soils typically incorporate some form of strain
softening. Moreover, many plasticity models for frictional materials use a non-
associated flow rule. Strain softening and non-associated flow rules can cause loss of
well-posedness of the initial-value problem, which can lead to a severe mesh depen-
dence in simulations and poor convergence of the iterative solution procedure. The
inclusion of viscosity, which is a common property of materials, seems a natural
way to restore well-posedness, but the mathematical properties of a rate-dependent
model, and therefore the effectiveness with respect to the removal of mesh depen-
dence, can depend strongly on how the viscous element is incorporated. Herein,
we show that rate-dependent models which are commonly applied to problems in
the Earth’s lithosphere, such as plate tectonics, is very different from the approach
typically adopted for more shallow geotechnical engineering problems. We analyse
the properties of these models under dynamic loadings, using dispersion analyses
and one-dimensional finite difference analyses, and complement them with two-
dimensional simulations of a typical strain localisation problem under quasi-static
loading conditions. Finally, we point out t hat a combined model, which features
two viscous elements, may be the best way forward for modelling time-dependent
failure processes in the deeper layers of the Earth, since it enables modelling of
the creep characteristics typical of long-term behaviour, but also regularises the
initial/boundary-value problem.
KEYWORDS:
viscoplasticity, strain softening, regularisation, failure, non-associated flow
1
INTRODUCTION
Strain localisation occurs in all geological materials, for a wide range of length scales, varying from faults which are hundreds of
kilometers long
1,2
, to decimeter-size laboratory samples
3,4,5
, and has been documented for a wide variety of materials, including
metals and marble, as early as 1931
6
. Strain localisation refers to the phenomenon of the concentration of strains in narrow zones
when the applied load exceeds a certain threshold level. It not only manifests itself at a broad range of length scales, but also t he
temporal scales can differ widely. While for most buildings time-dependent deformations and the ensuing damage occur at the
scale of several decades, geological processes within the lithosphere e.g., rock faults and shear zones, can stretch over millions
of years.
Strain localisation poses a challenge to one of the cornerstones of classical continuum mechanics, namely the assumption that
the mechanical and physical behaviour in a point is representative for a small, but finite volume sur rounding it. This assumption,
2 de Borst and Duretz
while mostly correct, and ubiquitous in geomechanics and geophysics, fails for highly localised deformations, like fault move-
ment or shear bands. A salient characteristic of the ensuing continuum is that it is devoid of an internal length scale, leading to
the prediction of zero-width localisation bands.
In fact, much of the observed structural softening when testing, for instance, laboratory-sized specimens under smooth, well-
controlled boundary conditions, is due to inhomogeneous deformations, in particular strain localisation, that develop in the
specimen at failure
4,5
. The challenge is therefore to separate the observed structural softening into a contribution due to strain
localisation and another by the ’true’ material behaviour. The scale of this difficulty is illustrated by the simple example of
non-associated plasticity. Even without the explicit introduction of strain softening at material level, i.e. for hardening or ideal
plasticity, numerical simulations can show structural softening
7,8,9,10,11
.
In practice, the introduction of non-associated flow or strain softening at a material level is unavoidable when carrying out
simulations based on the concept of a continuum. However, unless a continuum (plasticity) model which features strain soft-
ening, or non-associated flow, is equipped with an internal length scale and can thus predict a finitely-sized localisation zone,
it will lose ellipticity under (quasi)static loading conditions, and hyperbolicity under dynamic loading conditions, rendering
the initial/boundary-value problem ill-posed. This observation is the underlying reason of the reported mesh sensitivity and
the often erratic and unsatisfactory convergence behaviour of the equilibrium-searching iterative procedure
12,13
, which can be
aggravated upon mesh refinement
14
. A number of approaches have been proposed to repair this deficiency
15
, including Cosserat
plasticity
16,17,18
, non-local plasticity
19,20,21
and gradient plasticity
22,23,24
.
For a number of geomechanical applications, in particular those which involve very long-term processes in the Earth’s crust,
the consideration of a deformation-limiting viscosity represents an alternative to non-local models. It is emphasised though,
that not all rate-dependent rheologies solve the issue of mesh dependence, and that a pure series arrangement of the rheological
elements
25,26,27
, which is typically used for modelling the Earth’s mantle or lithosphere, see Figure 1(a) for instance, does not
introduce a length scale and therefore does not remove the mesh dependence. By contrast, a plasticity model which relies on the
introduction of a rate-limiting viscosity in a parallel ar rangement with a plastic slider
28,29
, Figure 1(b), does introduce a length
scale and can therefore provide mesh-independent numerical solutions
29,30
.
The aim of this contribution is to elucidate that different elasto-visco-plastic rheologies lead to continuum descriptions which
have very different mathematical properties. After a succinct recapitulation for a rate-independent elasto-plastic model, we
derive the governing differential equation for each of the rate-dependent approaches in a one-dimensional context. We apply a
dispersion analysis and verify this analysis with a finite difference computation. These analytical and numerical findings give
insight in the possible dispersive properties of the wave propagation and the presence, or not, of an internal length scale, which
is a condition for obtaining shear bands with a finite width. The analyses are augmented by two-dimensional finite difference
analyses of shear band localisation under quasi-static loading conditions, and confirm the statements regarding mesh dependence
for the two rate-dependent models. They also show an improved convergence behaviour of the (non-linear) iterative procedure
for a well-posed problem. To point out a possible way forward in modelling time-dependent localised deformations in the deeper
layers of the Earth’s crust, we have also examined a model which involves two viscous elements. Such a rheology can capture
the viscous creeping mode as well as brittle yielding of deforming rocks.
2
DISPERSION ANALYSIS AND INTERNAL LENGTH SCALE
2.1 Rate-independent elasto-plasticity
We consider a one-dimensional bar with a uniform cross-section 𝐴 and a constant mass density 𝜌. Assuming small deformation
gradients the balance of momentum is given by
𝜕𝜎
𝜕𝑥
= 𝜌
𝜕
2
𝑢
𝜕𝑡
2
(1)
while the kinematic relation reads
𝜀 =
𝜕𝑢
𝜕𝑥
(2)
with 𝜎 the axial stress and 𝑢 the axial displacement. As usual, 𝑥 and 𝑡 denote the spatial coordinate and time, respectively. As an
introductory step, we consider a rate-independent elasto-plastic solid, so that after the onset of plastic yielding, i.e. when 𝜎 = 𝜎
y
,
with 𝜎
y
the yield strength, the strain 𝜀 can be partitioned into an elastic component and a plastic component, as follows:
𝜀 = 𝜀
e
+ 𝜀
p
(3)
de Borst and Duretz 3
σ
ε
E
e
ε
σ
ε
η
v
v p
σ
ε
E
e
σ
η
ε
vp
vp
σ
σ
FIGURE 1 Rheological models: (a) Top: visco-elasto-plasticity; (b) Bottom: elasto-viscoplasticity
The stress is related to the elastic strain following Hooke’s relation, i.e.
𝜎 = 𝐸𝜀
e
(4)
with 𝐸 the time and spatially independent Young’s modulus. In a one-dimensional context and for monotonic loading, the plastic
strain rate can be derived by dividing the stress rate by the plastic hardening modulus ℎ:
𝜀
p
=
𝜎
ℎ
(5)
Combining Eqs (3)–(5) leads to:
𝜎 =
𝐸ℎ
𝐸 + ℎ
𝜀 (6)
and combination with the momentum balance, Eq. (1) results in:
𝜌
𝐸 + ℎ
𝐸
𝜕
2
𝑢
𝜕𝑡
2
= ℎ
𝜕
2
𝑢
𝜕𝑥
2
(7)
Next, we assume a harmonic wave,
𝑢(𝑥, 𝑡) = 𝐴𝑒
i(𝑘𝑥−𝜔𝑡)
(8)
with 𝑘 the wave number counting the number of wave lengths 𝓁 in the bar over 2𝜋, i.e. 𝑘 = 2𝜋∕𝓁, and 𝜔 the angular frequency.We
substitute t his expression into Eq. (7) to give:
𝜔 = 𝑘
𝐸ℎ
𝜌(𝐸 + ℎ)
(9)
or, using the definition for the phase velocity
𝑐
𝑓
=
𝜔
𝑘
=
𝐸ℎ
𝜌(𝐸 + ℎ)
(10)
A special case is the elastic bar velocity, which is obtained for a purely elastic material, i.e. ℎ → ∞, so that 𝑐
𝑓
→
𝐸∕𝜌. The
more important observation is that when ℎ < 0, t he phase velocity becomes imaginary, implying non-real wave speeds and a
wave propagation problem which is no longer hyperbolic, but elliptic, cf. Eq. (10)
2.2
Visco-elasto-plasticity
For sufficiently long time scales, all materials exhibit observable time-dependent deformations. For soft clays this holds for time
scales that equal the life time of structures, while for hard rocks viscous behaviour shows up at geological time scales. For the
4 de Borst and Duretz
latter type of applications it has therefore become customary to include a viscous component in the description of rocks for their
modelling at geological time scales
25,26
. The series arrangement of Eq. (3) is thus extended to include a viscous component:
𝜀
v
=
𝜎
𝜂
v
(11)
with 𝜂
v
the viscosity of the material, such that:
𝜀 = 𝜀
e
+ 𝜀
p
+ 𝜀
v
(12)
Using Eqs (4), (5), and (11) the strain rate reads:
𝜀 =
𝐸 + ℎ
𝐸ℎ
𝜎 +
𝜎
𝜂
v
(13)
Rearranging, differentiating with respect to 𝑥 and using the kinematic relation (2) gives:
𝜕
2
𝜎
𝜕𝑥𝜕𝑡
=
𝐸ℎ
𝐸 + ℎ
𝜕
3
𝑢
𝜕𝑥
2
𝜕𝑡
−
1
𝜂
v
𝜕𝜎
𝜕𝑥
(14)
Using the balance of momentum, Eq. (1) and integrating with respect to 𝑡 gives:
𝜕𝜎
𝜕𝑥
=
𝐸ℎ
𝐸 + ℎ
𝜕
2
𝑢
𝜕𝑥
2
−
𝜌
𝜂
v
𝜕𝑢
𝜕𝑡
(15)
while using the balance of momentum again results in:
𝜌
𝐸
𝐸 + ℎ
𝜕
2
𝑢
𝜕𝑡
2
− ℎ
𝜕
2
𝑢
𝜕𝑥
2
+
𝜌ℎ
𝜂
v
𝜕𝑢
𝜕𝑡
= 0 (16)
Evidently, in the rate-independent limit, i.e. when 𝜂
v
→ ∞ the last term cancels and we retrieve the governing equation for an
elasto-plastic solid, Eq. (7).
When comparing Eqs (7) and (16) it is clear that the viscous term in the rheological model has resulted in the addition of a
first-order term in the governing equation. Since the character of a partial differential equation is determined by the higher-order
terms
31
, no change is induced by this first-order term and it cannot be expected that the addition of the damper in the rheological
model will affect the change in character from hyperbolic to elliptic upon the introduction of strain softening (ℎ < 0).
To further investigate this, we substitute a single, linear harmonic wave, Eq. (8), into Eq. (16). This results in:
−𝜌(𝐸 + ℎ)𝜔
2
+ 𝐸ℎ𝑘
2
−
𝜌𝐸ℎ
𝜂
v
𝜔i = 0 (17)
No solution is possible when the wave number 𝑘 is assumed to be real. Hence, we assume 𝑘 = 𝑘
𝑟
+ 𝛼i to be complex, which
corresponds to the propagation of an exponentially attenuated harmonic wave:
𝑢(𝑥, 𝑡) = 𝐴𝑒
−𝛼𝑥
𝑒
i(𝑘
𝑟
𝑥−𝜔𝑡)
(18)
Substitution of the complex expression for 𝑘 into Eq. (17) and separating the resulting equation into real and imaginary parts
gives:
−𝜌(𝐸 + ℎ) 𝜔
2
+ 𝐸ℎ(𝑘
2
𝑟
− 𝛼
2
) = 0 (19a)
𝑘
𝑟
𝛼 −
𝜌𝜔
2𝜂
v
= 0 (19b)
Multiplying the first identity with 𝑘
2
𝑟
and substituting the second identity allow to solve for 𝑘
2
𝑟
:
𝑘
2
𝑟
=
𝜌𝜔
2
(𝐸 + ℎ)
2𝐸ℎ
1 +
1 +
1
(𝜂
v
𝜔)
2
𝐸ℎ
𝐸 + ℎ
2
(20)
From the latter identity we directly obtain the phase velocity:
𝑐
𝑓
=
𝜔
𝑘
𝑟
=
2𝐸ℎ
𝜌(𝐸 + ℎ)
⋅
1
1 +
1 +
1
(𝜂
v
𝜔)
2
𝐸ℎ
𝐸+ℎ
2
(21)
In the absence of plastic straining, or when ℎ is positive, Eq. (21) shows that waves travel at a lower speed due to the presence
of a damper because t hen
1
(𝜂
v
𝜔)
2
𝐸ℎ
𝐸+ℎ
2
> 0. Moreover, wave propagation is dispersive, since according to Eq. (21), the phase
velocity 𝑐
𝑓
is a function of 𝜔. Nevertheless, for strain softening, i.e. when ℎ < 0, the phase velocity still becomes imaginary.
