A simple model for tertiary creep in geomaterials Open Access

A

total cross-sectional area of a surface within the unit cell

$An+1(k)$

algorithmic modulus

A_s

solid matrix area within A

C

friction angle

c₁

material parameter related to cohesion

c₂

material parameter related to the friction angle

c₃

material parameter related to the dilation angle

D

dissipation rate

D

viscoplastic-damage consistent tangent modulus

D^e

standard isotropic elasticity tensor

E

Young’s modulus

$e˙p$

plastic deviator strain rate

F

Helmholtz free-energy function

f

yield function

G

shear modulus

I

second-order unit tensor

K

bulk modulus

k

iteration counter

m

order of the dissipation function D

n

material constant

p

effective mean stress

q

deviatoric stress invariant

r_d

material constant that governs the ratio of damage

r_p

material constant that governs the ratio of viscoplasticity

t

time

α_d

material damage occurring with a representative continuum volume element

α_p

micro-plastic strain

$α˙v$

plastic volumetric strain rate

β

softening/hardening parameter

ϵ

total strain tensor

ϵ^e

elastic strain tensor

$ϵ˙p$

plastic strain rate

ε_f

tolerance

ε_r

relative tolerance

η

viscous coefficient that controls the extent of plastic strain

Λ

arbitrary Lagrangian multiplier

ν

Poisson’s ratio

Π(α_d)

hardening/softening function

σ

true stress tensor

φ

angle of friction

χ

generalised stress tensor

$χ¯$

dissipative stress tensor

ψ

dilatancy angle

Creep may be defined as continued deformation without a stress change. Creep has been studied since about 1905, although such behaviour was documented as early as 1833 (Griggs, 1939). Most early studies focused on the creep rupture of metals under tensile stress. However, later studies have been carried out on rocks, in particular salt rocks, as these soft rocks creep under temperature and stress conditions, as evidenced from laboratory data (Hayano et al., 2001). Determining the creep characteristics of rock is an important stage in developing a tool that can predict the time-dependent deformation of an underground cavity. Creep tests performed in the laboratory are very significant in mining and improved design of underground structures (Cristescu, 1989). Creep testing of rock in the laboratory has been carried out by a number of researchers (Heap et al., 2010; Langer, 1982; Le Comte, 1965; Li and Xia, 2010; Obert, 1965; Phueakphum et al., 2010; Scott-Duncan and Lajtai, 1993; Singh, 1975; Vouille et al., 1984; Yang and Jiang, 2010). The simplest creep tests are those during which the rock specimen is uniaxially loaded in compression. The testing procedure involves an increment of load applied quickly to the rock specimen, and the stress is held constant while the gradually increasing strain is recorded regularly (Goodman, 1989). Triaxial tests have also been carried out in which the sample is confined by an all-around pressure, more closely simulating in situ conditions. The duration of an individual creep test is generally several weeks or months. Tests lasting a number of years have also been reported (Cristescu, 1989).

Laboratory data from creep tests are mostly depicted in the form of strain–time curves of which the general form is displayed in Figure 1 (Goodman, 1989; Jaeger et al., 2007; Jeremic, 1994). An instantaneous elastic strain, ϵ_e, is followed by primary or transient creep (region 1) in which strain occurs at an ever-decreasing rate. Secondary creep (region 2) follows if the constant stress overcomes a given limit and is characterised by a constant strain rate. For higher constant stress levels, tertiary creep (region 3) is also observed, which is characterised by a strain rate increase with time and leads eventually to failure. As illustrated by Jeremic (1994), laboratory investigations have shown that removal of the applied load in region 1 at point P of Figure 1 causes the strain to decrease rapidly to point Q (change in strain equal to ϵ) and then asymptotically back to zero at point R. Thus, region 1 can be classified as viscoelastic. Removal of stress in region 2 at point T will result in permanent deformation (VO). Thus, region 2 can be classified as viscoplastic. Despite the classification of these regions as viscoelastic and viscoplastic being an idealisation, it is reasonable to think that at low levels of stress the material behaviour is roughly viscoelastic and at high levels the behaviour is viscoplastic. The tertiary stage of creep behaviour (region 3) appears due to progressive microcracking of the material and would result in a loss of strength and stiffness, which may eventually lead to failure and a complete loss of the load-carrying capacity of the material.

Creep behaviour can also be classified according to the acting stress; it is possible to divide creep behaviour into volumetric and deviatoric (or shear) creep (Tavenas et al., 1978). Volumetric creep is caused by constant volumetric stress, and deviatoric creep is caused by constant deviatoric stress. Generally, volumetric creep consists only of the primary phase of the creep deformation – that is, it tends to stabilise. Deviatoric creep may or may not consist of all three phases, depending on the shear mobilisation. If the deviatoric stress is low, then only the primary creep phase will appear, but after crossing some level of the shear mobilisation, the primary phase will be followed by the secondary phase, which can lead to the tertiary phase and creep rupture. Tertiary creep is often observed in soft rocks, and it represents a problem in a mining environment.

There are a great number of constitutive models developed to account for the time-dependent creep behaviour of geomaterials. Extensive reviews for constitutive modelling of creep can be found in the papers by Liingaard et al. (2004) and Grimstad et al. (2017). Existing constitutive models for creep behaviour can generally be classified into three categories (Liingaard et al., 2004): (a) elementary phenomenological models, (b) rheological models and (c) general stress–strain–time models.

In elementary phenomenological models, empirical relations are obtained by directly fitting the observed test data with simple mathematical functions. Logarithmic and exponential laws are often adopted to link creep strain to time. Examples of these models can be found in the papers by Singh and Mitchell (1968), Campanella and Vaid (1974), Leroueil (1987), Yin (1999) and Bi et al. (2019). One of the basic limitations of elementary phenomenological models is that they are strictly valid only for conditions that are identical to those of the test from which they were derived (e.g. one-dimensional (1D) condition). However, they can provide practical solutions to engineering problems, as far as the boundary conditions are consistent with the laboratory experiments.

Rheological models consist of a combination of different components, such as springs, dashpots and sliders, to describe the time-dependent creep behaviour of a material. The structure of these models is not related to a particular creep test, and therefore, only the model parameters change between tests in order to provide a fit for the strain–time data. Linear viscoelastic models consist of various combinations of two states of deformation, elastic behaviour (represented by a spring) and viscous behaviour (represented by a dashpot). Many different linear viscoelastic models have been proposed, such as the Maxwell, the Kelvin and the Burgers model. The simplest model that can be used to trace strain up to the onset of secondary creep is the Burgers model, which is composed of a Maxwell model and a Kelvin model connected in series (Figure 2(a)). An overview of these viscoelastic models can be found in the book by Jaeger et al. (2007). These viscoelastic models are typically formulated for 1D problems. However, they can be extended and applied to three-dimensional boundary-value problems. Birchall and Osman (2012), for example, presented analytical solutions for predicting the creep displacements of deep tunnels in three dimensions using the Burgers model. Linear viscoelastic models cannot simulate the failure of geomaterials. Slider elements (St Venant elements) are added to the elastic and viscous components of viscoelasticity in order to include the failure process. Typically, a slider element is placed in parallel with a dashpot element, known collectively as a Bingham unit (Figure 2(b)). The model by Perzyna (1966) is a well-known example of a viscoplasticity model. A viscoplastic material exhibits a time-dependent behaviour only in the plastic region. This behaviour corresponds to the secondary stage of creep. However, these models cannot replicate tertiary creep behaviour.

Generalised stress–strain models can describe both viscous effects and rate-independent behaviour of soils under general loading conditions. Several theoretical approaches have been used to derive these models. Borja and Kavazanjian (1985), for example, used the concept by Bjerrum (1967) of total strain decomposition into an immediate (time-independent) part and a delayed (time-dependent) part to model the stress–strain–time behaviour of wet clays. A modified Cam-Clay yield surface (Roscoe and Burland, 1968) is employed to characterise the time-independent stress–strain behaviour. Kaliakin and Dafalias (1990) presented a bounding-surface model for the time-dependent response of overconsolidated soils. The basic characteristic of the bounding-surface approach is that the material exhibits a memory of the loading history and the plastic deformation depends on the distance between the present stress point and an image stress point on the bounding surface defined by a mapping rule for stress states within the surface. Hashiguchi and Okayasu (2000) adopted the subsurface modelling approach to simulate the time-dependent behaviour. This model assumes a bounding surface and a subloading surface that always passes through a current stress point. The model can describe the continuous stress rate–strain rate relation in the loading process with smooth elastic–plastic transition. Despite the complexity of the generalised models mentioned earlier, numerical investigations have focused on volumetric creep under undrained conditions, despite the fact that creep is an inherently drained phenomenon, as it takes a long time to develop (Jardine, 2014, Kavvadas and Kalos, 2019). At large shear stress approaching strength, deviatoric creep strain rates accelerate and lead to tertiary creep and failure (Bishop, 1966).

The literature on constitutive formulations for predicting tertiary deviatoric creep is sparse. Dragon and Mroz (1979), Al-Shamrani and Sture (1998), Pellet et al. (2005) and Shao et al. (2003) modelled tertiary creep using a strain-softening/hardening approach based on the shear-stress-driven damage mechanism. Gioda and Gividini (1996) predicted tertiary creep by introducing viscous damage on a Bingham model. Kavvadas and Kalos (2019) presented a time-dependent constitutive model for cohesive structured soils (TMS model). The TMS model is constructed from three characteristic surfaces: (a) an internal plastic yield envelope (PYE), which is a small surface bounding the elastic domain; (b) an external structure envelope (SE) bounding the PYE and all accessible states; and (c) an intrinsic envelope, corresponding to an equivalent intrinsic state. The size, shape and location of the SE give the magnitude of the structure (and strength), which increases with viscous volumetric strains and degrade with plastic strains and viscous shear strains. The TMS model is capable of simulating tertiary creep; however, it requires identification of 21 material constants.

The main aim of this work is to model tertiary creep using a sound theoretical framework. The new simple model, presented here, makes use of damage mechanics concepts. The damage mechanics framework is a powerful tool that can be used to simulate many features in geomechanics, as illustrated by Xiao and Liu (2017) and Xiao and Desai (2019a, 2019b). This paper extends the hyperplasticity approach of Collins and Houlsby (1997), Houlsby and Puzrin (2000) and Puzrin and Houlsby (2001) to model coupled viscoplastic-damage frictional materials. This paper focuses on deviatoric creep and presents a model that obeys the laws of thermodynamics. The entire constitutive behaviour can be derived from two scalar potentials: a free-energy potential that provides the elasticity law and a dissipation potential that provides the yield function, the direction of plastic flow and the evolution of a damage variable. No additional assumptions are required. The motivation here is to derive a simple model. Thus, the new model has only a few parameters that have physical meanings and is capable of capturing the tertiary creep observed in soft rocks.

For isothermal deformations of a continuum undergoing small strains, the first and second laws of thermodynamics can be written as

where F represents the Helmholtz free-energy function; D is the dissipation rate, which is strictly non-negative; and σ denotes the true stress tensor and ϵ denotes the total strain rate tensor. Collins and Houlsby (1997) and Houlsby and Puzrin (2000) introduced the hyperplasticity theoretical framework, in which they demonstrated that the material behaviour can be determined from the knowledge of the two thermodynamic potential functions F and D. The local state of the material is assumed to be defined completely by the knowledge of (a) the strain tensor ϵ, (b) a set of internal tensor variables and (c) the entropy. For the isothermal case, the entropy does not enter into the formulation.

For a damaged material, the internal variables can, most usefully, be thought of as the micro-plastic strains α_p and the material damage α_d occurring within a representative continuum volume element. Therefore, free energy is chosen to be of the form F(ϵ, α_p, α_d), while the dissipation rate function is taken to be $D(ϵ˙,αp,α˙p,α˙p,α˙d)$⁠. If the dissipation function is homogenous of the first order, then Equation 1 can be rewritten as

It follows that

and

where the generalised stress tensors χ_p and χ_d are given by

and the dissipative stress tensors $χ¯p$ and $χ¯d$ are given by

Following the orthogonality condition by Ziegler (1983), which was adopted in the hyperplasticity framework of Houlsby and Puzrin (2000), a wide range of classes of materials can be described by enforcement of the stronger condition in Equation 4

To develop rate-independent thermomechanical models, Puzrin and Houlsby (2001) suggested a decoupled form of dissipation function that is appropriate for multiple surface kinematic hardening plasticity models. For a rate-independent damage material, the dissipation function takes the form of

As noted by Einav et al. (2007), adopting a decoupled dissipation for rate-independent plasticity-damage models may result in damage prior to plastic straining or vice versa. Furthermore, the uncoupled dissipation potentials imply multi-yield surfaces, each one associated with the evaluation of each internal variable. Therefore, Einav et al. (2007) proposed a coupled dissipation function for rate-independent materials, which leads to a single yield function, of the form

where c_i(ϵ, α) is a positive definite function and $Φi(α˙i)$ is a homogenous first-order function operator returning a positive scalar. It should be emphasised that for a rate-independent material, the dissipation function is chosen to be a homogeneous first-order function of the internal variable rate. Thus, the dissipation function can be expressed mathematically through the Euler equation

In rate-dependent materials, the dissipation function is not a first-order function. Thus, Equation 9 and the formulations mentioned earlier cannot be directly applied. If D is of order m, then Z = D/m can be defined so that

and the dissipative stress tensor can be derived as

For more details about the mathematical formulation of the hyperplasticity framework, the reader is referred to the book by Houlsby and Puzrin (2006).

This paper presents a new viscoplastic-damage model that is capable of simulating secondary creep and tertiary creep. The viscoplastic component of the model describes the secondary, steady-state creep, while the damage component of the model can describe the tertiary creep during which the strain rate increases with time due to a reduction in material stiffness. Here the theory of continuum damage mechanics is used, which was developed by Kachanov (1958), who introduced an internal scalar variable to model the creep failure of metals under uniaxial loads. Other significant contributions to this theory were made by Lemaitre and Chaboche (1990), among others. The authors combined continuum damage mechanics within the framework of hyperplasticity, thus encompassing viscoplasticity and damage within a single theory.

In this paper, the attention is focused on the case of isotropic damage, so that α_d is simply a scalar internal variable starting from 0 and increasing to a maximum value of 1. The damage variable can be defined as

where A is the total cross-sectional area of a surface within the unit cell in one of the three perpendicular directions and A_S is the solid matrix within A.

The free-energy potential can also be written as follows

where ϵ^e is the elastic strain tensor defined as the total strain minus the plastic strain and D^e is the standard isotropic elasticity tensor. From this potential, the elasticity law is given by

The stresses associated with the damage can be derived from the free-energy potential as follows

where G is the shear modulus; K is the bulk modulus; $p=tr[σ]/3$ is the effective mean stress; $tr[⋅]$ is the trace operator of $[⋅]$⁠; and $q=(3J2)1/2$ is the deviatoric stress invariant, with $J2=[(1/2)(s:s)]1/2$⁠. $s=σ−pI$⁠, with I as the second-order unit tensor.

For a coupled viscoplastic-damage material, the authors propose a dissipation function of the form

The term inside the first brackets under the square root, in the expression of the dissipation function, D, represents the dissipation due to plasticity, while the term inside the second brackets represents the dissipation due to damage. The additive term under the square root is the dissipation due to viscosity. $c1=c1*Π(αd)$ and $c2=c2*Π(αd)$⁠, where $c1*$ and $c2*$ are material parameters related to the cohesion and the friction angle, respectively. Π(α_d) is a hardening/softening function, $α˙p1$ and $α˙p2$ are internal variables associated with the micro-plastic strain, $α˙d$ is the material damage, η is the viscous coefficient that controls the extent of plastic strain, n is a material constant, while r_p and r_d are two constants governing the ratio of viscoplasticity and damage, and ρ is defined as

The internal variables $α˙p1$ and $α˙p2$ can be taken as equal to the plastic deviatoric shear strain rate $α˙γ$⁠, which is defined as

where

and $e˙p$ is the plastic deviator strain rate; $ϵ˙p$ is the plastic strain rate; and $α˙v=tr[ϵ˙p]$ is the plastic volumetric strain rate.

The dissipation function (Equation 17) is consistent with the coupled form given by Equation 9 and therefore yields a single yield surface. This coupled dissipation function, therefore, defines a model that introduces damage whenever plasticity occurs and vice versa. As will be illustrated later, the definition of ρ given by Equation 18 implies that the damage starts with the plasticity behaviour. In the absence of damage in a rate-independent material, this dissipation function collapses to the standard form of dissipation in a Drucker–Prager-type material (i.e. $D=p(c1+c2p)α˙γ$⁠).

A linear relationship between the volumetric and shear strain rates can be introduced as follows

where $c3=c3*Π(αd)$ is a material parameter related to the dilation angle. Therefore, the modified dissipation function can be written as

where Λ is an arbitrary Lagrangian multiplier.

Since the dissipation given by Equation 22 is the summation of homogeneous functions of different orders, Equations 11 and 12 can be applied to obtain the dissipative stresses. It follows that the dissipative deviatoric stress $χ¯q$ is given by

and the dissipative stress due to damage is given by

and the dissipative mean stress is equal to the Lagrangian multiplier

Manipulating the preceding stress expressions (Equations 23–25) leads to the following equation

Since there is no kinematic hardening, the true stresses and the dissipative stresses are identical; from Ziegler’s orthogonality, the dissipative stresses and the generalised stresses are equal – that is

An expression for the damage internal variable in the true stress space can be found by substituting Equations 23 and 24 into Equation 18 and making use of Equation 27

Equation 26 can be rewritten in the true stress space as

Rearranging the preceding equation, the following is obtained

The plastic deviatoric strain rate $α˙γ$ is a non-negative quantity; therefore, the expression inside the braces in Equation 30 must be greater than or equal to zero. As a result, the yield criterion can be easily derived and is shown to be given by

where

For Drucker–Prager-type materials, the constants r_p and r_d are related by the following expression

This is consistent with the definitions of r_p and r_d in a rate-independent material as proposed by Einav et al. (2007). It follows that r_p ≥ 1 and r_d ≥ 1. In the limiting case when r_p = 1, r_d → ∞ and the constitutive model simplifies to a viscoplastic model. It should be noted that Equation 17 implies that

Thus, comparing this with the definition of ρ given by Equation 18, the expression for the yield surface (Equation 31) can be retrieved. Thus, damage can start only with the start of the plastic behaviour.

By applying the dilatancy constraint of Equation 21 and using Equation 20, it can be shown that for a Drucker–Prager material, the components of the plastic strain rate are given by

The preceding relationship can be written as

where g = q − c₃p and $γ˙=α˙γ$⁠. For convenience, the set of constitutive viscoplastic-damage equations is summarised in Table 1.

In order to use the proposed damage viscoplasticity model in a displacement-based finite-element formulation, a stress integration procedure of the constitutive model, which takes place at the elemental Gauss points, is required. Here the iterative implicit backward Euler stress integration scheme based on the operator split methodology is used (Simo and Taylor, 1985). For a strain increment Δϵ over a time step Δt and the state variables at t_n, the updated stress vector σ_n+1 and the damage parameter $αdn+1$ can be obtained. This involves two basic sequences over a time interval [t_n, t_n+1] – namely, an elastic predictor and a plastic corrector, depending on whether the elastic trial stress falls inside or outside the yield function f.

During the elastic predictor, the plastic internal variables are kept constant and equal to their respective values at time t. For the initial iteration count k = 0, these variables are

where the subscript (n + 1) indicates that all values are obtained at the end of the increment. Furthermore, the elastic trial stress is expressed as

At this stage, if the trial state is admissible (i.e. $f(n+1(0),αdn+1(0))≤ϵf$⁠, where ϵ _f is a tolerance factor) then $(⋅)n+1=(⋅)n+1trial$⁠. Otherwise, the yield criterion is violated and plastic and damage strains are expected to occur.

The Newton–Raphson scheme, which can lead to efficient return mapping procedures (de Souza Neto et al., 2009), was adopted to evaluate the plastic corrector sequence. To this end, a system of non-linear implicit equations has to be solved simultaneously for the unknowns $an+1T={σn+1,αdn+1}T$ iteratively using the Newton–Raphson method. The discrete residuals associated with the unknowns are written as follows

where k is an iteration counter. In the linearisation process of Equation 39, the directional derivative of the system variables ${σn+1,αdn+1}$ in the direction of Δa _n +1, with Δϵ _n +1 being kept constant during the return mapping, implies that

where the components of $[An+1(k)]$ are the partial derivatives of the residual vector with respect to the non-linear system ${σn+1,αdn+1}$⁠. If the residuals are greater than a specified tolerance – for instance, $‖rn+1(k)‖/‖rn+1(0)‖<εr$⁠, where ε _r is a relative tolerance – Equation 39 is solved again to update the unknowns and iteration is continued until the residuals fall within the specified tolerance. For presentational convenience, Table 2 summarises more clearly the operations needed for the iterative implicit stress-update algorithm for the viscoplastic-damage model.

The next task consists of the derivation of the elasto-plastic tangent operator consistent with the algorithm employed to update the stress for use in finite-element computations. The viscoplastic-damage consistent tangent modulus is derived as

The computation of Equation 41 relies on the total derivatives $∂rn+1/∂ϵn+1$ of the implicit non-linear residual equations in Equation 39 and the inverse of [A], which can be extracted from the local iteration procedure of the return mapping algorithm (Rouainia and Muir Wood, 2001; Rouainia and Peric, 1998).

Many advanced constitutive models have been developed to reproduce some aspects of viscoplastic behaviour identified in the literature. However, these models are often characterised by their use of numerous parameters, many of which have no clear physical meaning and require extensive laboratory testing to be calibrated and be practical as a design tool. In addition to being thermodynamically consistent and derived entirely from two potentials without any additional assumptions, the proposed model requires only eight parameters to capture both the elastic and the coupled viscoplastic-damage behaviour: two for linear elasticity, three for the yield surface, two for the viscous behaviour and one for defining damage. These parameters have a physical meaning and can be obtained from relatively straightforward testing procedures.

The pre-yield behaviour is modelled linear elastically using Young’s modulus E and Poisson’s ratio ν, which can be obtained from the approximately linear part of a stress–strain curve obtained from laboratory testing (following a triaxial test, for example). The model requires the identification of three conventional parameters of the Drucker–Prager model (⁠$c1*$⁠, $c2*$ and $c3*$⁠), which are related to cohesion, C; friction angle, φ; and dilatancy angle, ψ. The cohesion, friction angle and dilatancy angle can be identified from conventional triaxial tests. In triaxial tests, the Drucker–Prager yield surface can be taken to circumscribe the Mohr–Coulomb yield surface as shown in Figure 3 in the deviatoric stress space. Therefore

In this model, the post-yield behaviour is described by the softening/hardening function Π(α_d) equal to (1 − βα_d)², where β is the softening/hardening parameter. This function fits well with experimental data, as demonstrated by Shao et al. (2003). The parameter β can be obtained from the post-yield stress–strain curve in triaxial tests; if the geomaterial softens, β > 0, and if it hardens, then β < 0.

The parameters r_p and r_d are related by Equation 33 and govern the damage as shown in Figure 4. When r_d/r_p = 0, the model performs as a classical viscoplastic model and the damage mechanism is deactivated. However, when 0 < r_d/r_p ≤ 1·0, the model is allowed to undergo damage and therefore softening. The ratio r_d/r_p can be estimated by measuring the damage in the material. The variations of the elasticity modulus can be found in a displacement control triaxial test by carrying out a post-yield unloading/reloading cycle and comparing the values of the material stiffness with that of its original value before the yield stress is reached. There are other techniques for measuring damage in materials reported throughout the literature (see e.g. the paper by Lemaitre and Dufailly (1987)). If a rapid test is carried out and the damage α_d is plotted against the accumulative deviatoric plastic strain α_γ, then the ratio r_d/r_p can be estimated from the initial slope of the α_γ − α_d curve. From Equation 28

Two parameters are used to describe the viscoplasticity response; viscosity is defined by the parameter η, and the rate-sensitivity parameter is described by the non-dimensional parameter n. Both parameters are positive. These parameters are analogous to the parameters of the model by Perzyna (1966), which are widely used in computational applications of viscoplasticity and can be obtained using triaxial tests from the relationships given by Equation 30.

The model has been specifically tested for its ability to reproduce a series of triaxial compression test results on three different geomaterials reported in the literature: (a) sandstone (Yang and Jiang, 2010), (b) Meuse–Haute/Marne claystone (Hu et al., 2014) and (c) Inada granite (Fujii et al., 1999).

Yang and Jiang (2010) carried out drained triaxial tests on sandstone core samples extracted at the Xiangjiaba Hydropower Project in south-western China. The samples were taken at a depth of between 61·2 and 62·5 m, where the Poisson’s ratio and density are approximately ν = 0·22 and 2·600 Mg/m³, respectively. Triaxial loading tests under drained conditions performed on isotropically compressed samples were used to illustrate the short-term and the long-term behaviour.

The short-term loading tests were conducted under three different confining pressures (σ ₃ = 3, 5 and 7 MPa), and axial stresses were imposed at a constant rate of 0·127 MPa/s until peak strength was reached. The measured angle of friction φ was found to be about 58·4°. Table 3 shows the mechanical behaviour observed under triaxial compression. The peak axial strain ϵ _1p is defined as the strain value at rupture in terms of the complete stress–strain curve. The modulus of elasticity (E _s) refers to the slope of the approximately linear part of the stress–strain curve, while the modulus of deformation (E ₅₀) is defined as the slope between the original point and the stress at half-peak strength.

The parameters $c1*$⁠, $c2*$ and $c3*$ of the viscoplastic-damage model are calculated from the measured angle of friction φ and a dilatancy angle of φ/4 using Equation 42. Hoek and Brown (1997) studied the relationship between the angle of dilation ψ and the angle of friction φ and recommended that for a good-quality sample, ψ can be taken as φ/4. The average value of the Young’s modulus $(E50+Es)/2$ for each sample is used in the numerical simulation (see Table 4). In the absence of the detailed testing procedure (described earlier) for obtaining the model parameters, the short-term loading test data for a confining pressure of 5 MPa are used to obtain the remaining model parameters. A Matlab software built-in standard Levenberg–Marquardt optimisation procedure has been utilised to estimate the set of material parameters. Table 4 shows the resulting set of material parameters used in the simulations.

Figure 5 shows the comparison between the short-term drained stress–strain response of the model and the experimental data with three different confining pressures. As can be seen from Figures 5(a)–5(c), the new coupled viscoplastic-damage model is remarkably successful in matching the general shape of the stress–strain response for all three confining pressures, at both peak stress and in the region of material softening.

Validation was also carried out with drained triaxial compression tests of saturated Meuse–Haute/Marne claystone. These tests were undertaken by Hu et al. (2014) on samples taken from a depth of about 521 m and were subjected to an in situ vertical stress of 12·7 MPa and horizontal stresses in the range 12·3–13·8 MPa. The axial strain rate in these tests was set at 10⁻⁷/s (corresponding to an average stress rate of 10⁻⁴ MPa/s). The measured angle of friction φ was found to be about 21·4°, and the cohesion c was found to be 3·94 MPa. The optimised parameters used in the analysis are shown in Table 3. Figure 6 shows the relationship between the radial strain and the axial strain at different confined pressures. This figure demonstrates that except for the case of high confining pressure (20 MPa), the new model can reasonably replicate the experimental data.

Table 5 shows the evolution of the relative residual norm for a confining pressure of 5 MPa for three typical load increments. These results provide a clear indication of the good rates of convergence. The number of iterations needed to meet the convergence tolerance of ϵ _r = 10⁻⁸ is between three and four per load step.

Yang and Jiang (2010) also conducted creep tests on sandstone in triaxial cells. In these tests, the deviator stress q = (σ ₁ − σ ₃) was kept constant. Figure 7 shows the comparison between the results of the viscoplastic-damage model and the experimental data for creep tests with a confined pressure of 5 MPa and a constant deviator stress of 160 MPa. The failure occurs after approximately 31·48 h. It can be seen that the general trend is well captured in terms of the time required to reach failure and the magnitude of axial strain. The evolution of the accumulative damage predicted by the model during the creep test is shown in Figure 8. It can be noticed that the accelerated damage rate follows a pattern similar to that of the accelerated creep strain. This is consistent with the experimental data of Heap et al. (2010), which show an accelerated rate of change in the pore volume of sandstone during tertiary creep.

Validation is also made against triaxial creep tests on Inada granite reported by Fujii et al. (1999). These tests were conducted under a confining pressure of 10 MPa. The testing programme investigated the effect of creep stress on the strain–time behaviour. The numerical simulations were carried out with the parameters for Inada granite shown in Table 4. Figure 9 demonstrates that the viscoplastic damage model is capable of predicting the time required to failure at different stress levels.

Inherent or structural anisotropy is observed in many types of sedimentary rocks. The strength of the rock material strongly depends on the loading orientation with respect to the microstructure (Pietruszczak et al., 2002). The model presented here deals with isotropic damage only. Therefore, it is applicable only to the case of monotonic loading. The present formulation assumes that the dilatancy angle is constant and that there is a unique relationship between the plastic shear strain and the volumetric strain increments. Alejano and Alonso (2005) showed the dependency of the dilatancy angle on the confining stress and the material plasticity. The model is formulated in terms of triaxial stress space, and validation is carried out against triaxial data. It could be generalised by incorporating an appropriate Lode’s angle dependency function such as that suggested by Eekelen (1980). The strength parameters c ₁, c ₂ and c ₃ can be expressed as a function of Lode’s angle. However, further validation is needed. Furthermore, temperature dependency could be easily incorporated within the hyperplasticity approach as explained by Houlsby and Puzrin (2000).

In this paper, the authors have combined the continuum damage approach within the framework of hyperplasticity for rate-dependent materials, thus encompassing viscoplasticity and damage within a single theory. A novel coupled viscoplastic-damage constitutive model has been derived from two scalar potentials: a free-energy potential that provides the elasticity law and a dissipation potential that provides the yield function, the direction of plastic flow and the evolution of a damage variable. No additional assumptions are required. The model is formulated with few parameters that can be derived from conventional laboratory testing. It has been shown that the new model is capable of capturing the secondary and tertiary creep observed in soft rocks. Comparison with triaxial data has shown that the model provides a good approximation to short-term stress–strain curves obtained for different confining pressures and to long-term creep data.

The authors would like to thank the UK Engineering and Physical Sciences Research Council (EPSRC) for the PhD studentship awarded to the second author (Grant reference: EP/D504376/1).