Modified Cam Clay bounding surface hyper-viscoplastic model

Dadras-Ajirlou, Davood; Grimstad, Gustav; Ghoreishian Amiri, Seyed Ali; Degago, Samson Abate; Houlsby, Guy T.

doi:10.1680/jgeot.23.00308

Clays exhibit complex mechanical behaviour with significant viscous, non-linear and hysteric characteristics, beyond the prediction capacity of the well-known modified Cam Clay (MCC) model. This paper extends the MCC model to address these important limitations. The proposed family of models is constructed entirely within the hyperplasticity framework deduced from thermodynamic extremal principles. More specifically, the previously developed MCC hyper-viscoplastic model based on the isotache concept is extended to incorporate multiple internal variables and to capture recent loading history, hysteresis and smooth response of the material. This is achieved by defining an inelastic free energy and an element that implements a bounding surface within hyperplasticity, resulting in pressure dependency in both reversible and irreversible processes with a unique critical state envelope, and only eight material parameters with a readily measurable viscous parameter. A kinematic hardening in the logistic differential form in stress space is derived that enables the proposed model to function effectively across a wide range of stresses. Based on this kinematic hardening rule, the current stress state acts as an asymptotic attractor for the back or shift stresses whose evolution rates are proportional to their current state.

INTRODUCTION

The modified Cam Clay (MCC) model (Burland, 1965; Roscoe & Burland, 1968), deduced from the unified and comprehensive behavioural framework of critical state soil mechanics (CSSM) (Roscoe et al., 1958; Schofield & Wroth, 1968), revolutionised the understanding of the mechanical behaviour of soil, particularly clay, by linking the compressional and shearing behaviour. Although it was originally energy based, with the normality rule for the inelastic flow imposed by Drucker's stability postulate (Drucker, 1959), the MCC model has also been re-appraised within the context of modern thermodynamics. Constitutive models lacking thermodynamic validity should not be confidently used for solving boundary value problems under various loading conditions, especially those involving unloading–reloading cycles. This is because they may spuriously generate or lose energy and subsequently violate the principle of conservation of energy. For instance, studies by Zytynski et al. (1978) and Borja et al. (1997) demonstrated that a system represented by such non-conservative models may not return to its initial state after a reversible process, due to a false generation or dissipation of energy.

Shortly after its introduction, the MCC model received a thermodynamic description (Houlsby, 1981), laying the foundation for its further development in a rigorous, systematic and energy-based manner. Notably, improvements have been made to the MCC model, such as lifting the unnecessary normality rule (Collins & Hilder, 2002; Collins & Kelly, 2002) while maintaining maximum dissipation. This has been done just by introducing a convex and versatile dissipation potential coupled with the current state of the material. In these extensions, the approach of thermodynamics with internal variables deduced from the extremum principles (Ziegler, 1977) has been employed. The systemised version of the approach, termed ‘Hyperplasticity’, which allows a systematic development of models from a basic version, has been presented by Houlsby & Puzrin (2006). Within the hyperplasticity framework, all components of a model are integrated into two characteristic potential functions: the force/dissipation potential defining the path-dependent/irreversible behaviour, and the free energy potential defining the path-independent/reversible behaviour.

In addition to the normality restriction, which could result in extreme dilatancy for some overconsolidated clays, the original MCC model suffers from other important limitations. Perhaps, as pointed out by its founders (Roscoe & Burland, 1968), the most profound deficiency, especially when prediction of soft clay behaviour is concerned, is the lack of time concept, preventing spontaneous change of the material state with a lapse of time. Adachi & Okano (1974), Kutter & Sathialingam (1992) and Yin & Graham (1999), among others, addressed the time-independency of the MCC model using the overstress theory of Perzyna (1963, 1966), which is also based on the restrictive stability postulate of Drucker (Drucker, 1959), resulting in the normality restriction. Using the hyperplasticity approach (Houlsby & Puzrin, 2006), a viscoplastic MCC model has recently been developed (Dadras-Ajirlou et al., 2023a, 2023b) that lifts the normality restriction of inelastic viscous flows while the dissipation potential is maximal and, simultaneously, the critical state (CS) envelope is unique under different loading rates. In this model, the viscous effects of clay behaviour are captured by way of the isotache concept (Šuklje, 1957; Leroueil, 2006), referring to a unique relationship between the inelastic strain rate (isotache) and the current state of the material. The versatility of the developed model has been enhanced by incorporating the spacing ratio (Wroth & Houlsby, 1985; Collins & Hilder, 2002), normally defined as the ratio of the normal consolidation pressure to the pressure at the CS (Yu, 1998). This addition allowed the regulation of shearing isotaches between the contractant (wet) and the dilatant (dry) zones, separated by a unique CS envelope. The inclusion of these features through the dissipative elastoplastic coupling (Collins, 2002) enabled the hyper-viscoplastic model to replicate effectively the general behaviour of clay under monotonic loading with different rates, including stress relaxation and creep. This capability has been exemplified in the simulation of the well-documented behaviour of Hong Kong marine deposits (HKMD) (Yin & Zhu, 1999; Zhu, 2000; Yin et al., 2002).

Other fundamental limitations of the MCC model and its offshoots are the abrupt stiffness degradation and the lack of memory of the recent loading history, failing to capture the hysteric behaviour during unloading–reloading cycles. Several efforts (e.g. Mróz et al., 1979; Stallebrass & Taylor, 1997; Houlsby, 2000; Kavvadas & Amorosi, 2000; Rouainia & Muir Wood, 2000; Einav & Puzrin, 2003) have been made to address these important features of clay behaviour – however, without considering the time and rate effects. More recently, there has been a focus on capturing the hysteresis behaviour of clay within the context of viscoplasticity, (e.g. Kimoto et al., 2015; Maranha et al., 2016; Jiang et al., 2017; Tafili & Triantafyllidis, 2020; Yuan & Whittle, 2021; Bathayian & Maleki, 2023; Shi et al., 2023). Nonetheless, despite their capabilities and sophistication, some of these advanced viscoplastic models have a large number of material parameters. This motivates further research on constitutive modelling of soil behaviour with fewer material parameters, while addressing the mentioned limitations. In this regard, hyperplasticity stands out as a promising framework due to its stringent mathematical structure, ensuring unconditional thermodynamic consistency and minimising ad-hoc assumptions.

The above-mentioned important limitations of the phenomenological models, including the MCC model and its offshoots, could stem from the crude assumption of homogeneous plastic flow within the representative elementary volume (REV) (Coussy, 1995). This is particularly relevant for clays with complex inhomogeneous microstructure, which can both cause and be influenced by the inhomogeneous plastic deformation. Mesoscopic mechanical models (Einav & Collins, 2008; Houlsby, 2020) suggest that the inhomogeneity of plastic deformation within REV leads to phenomena like the Bauschinger effect and hysteresis behaviour. Collins (2005) interpreted the inhomogeneous behaviour of geomaterials within REV on the mesoscale, including clay, as a signature of the formation of a network of strong and weak force chains (Radjai et al., 1998). In such a setting, it is hypothesised that the reversible energy is locked in the weak network trapped between strong networks. The trapped reversible energy is referred to as stored inelastic energy since its release can only occur through changes in the network of force chains at the mesoscopic level, induced by reversed plastic loading at the macroscopic or continuum level.

In the hyperplasticity framework, the reversible inelastic energy is integrated into the free energy potential (Collins & Houlsby, 1997; Houlsby & Puzrin, 2006). Following the orthogonality principle introduced by Ziegler (1977), a crucial element of the hyperplasticity framework, the inclusion of reversible inelastic energy leads to a conservative stress quantity. This quantity is termed as the back or shift stress whose evolution in stress space is known as the kinematic hardening rule. The kinematic hardening or the back-stress concept has been widely adopted as an effective phenomenological modelling tool for capturing the hysteresis behaviour and Bauschinger effect (e.g. Mróz, 1967; Chaboche, 1989; Puzrin & Houlsby, 2001; Krabbenhøft & Krabbenhøft, 2021). Nevertheless, integrating the kinematic hardening and the back-stress concept into viscoplastic models within soil mechanics comes with some important consequences.

Grimstad et al. (2020) and Dadras-Ajirlou et al. (2023b) have demonstrated that incorporating back-stress and related kinematic hardening under isotache viscosity for a system with a single internal variable (tantamount to single yield or dynamic surface in the context of viscoplasticity) leads to non-uniqueness of the CS envelope under different loading rates, even in triaxial stress space. A unique envelope for the CS condition is an essential paradigm in CSSM that provides a reference for the unified description of the general mechanical behaviour of soils with different densities and history. A large body of test data (e.g. Arulanandan et al., 1971; Vaid & Campanella, 1977; Adachi et al., 1995; Sheahan et al., 1996; Zhu, 2000; Hicher, 2016; Tafili et al., 2021) strongly supports this observation, suggesting that the CS envelope is largely independent of the rate of mechanical processes. These experimental observations align with the widely used first-order approximation of Coulomb's sliding friction at macroscale (Popova & Popov, 2015), which is generally considered to be almost independent of the rate of the shearing in an isothermal process with no chemical reactions. The uniqueness of the CS envelope under various loading rates should be distinguished from the more subtle matter of the uniqueness of the CS envelope under different shearing modes (Lode angle) (Dafalias & Taiebat, 2013; Li & Dafalias, 2015).

Moreover, based on the conventional practice of the plasticity theory, back-stress may be associated with a certain elastic domain that shifts by movement of the back-stress according to a kinematic hardening rule. A general obstacle in this context is the establishment of the initial position of the back-stress in stress space for each material point within the boundary value problem's domain. Any viscoplastic processes, including pure creep or stress relaxation, necessitates the stress state to be situated outside the considered elastic domain (Chaboche, 2008).

Introducing multiple kinematic surfaces may entail additional viscous parameters (e.g. Shi et al., 2023). Each of these additional viscosities is associated with a specific surface. However, establishing and evaluating these extra viscous parameters through a geotechnical laboratory or in situ tests can be challenging, and may not be feasible. For a system with a single internal variable, expressing the scalar-valued dissipation function as sums of homogeneous functions with varying degrees of homogeneity results in a non-homogeneity of the dissipation function which requires a special mathematical treatment; see Puzrin & Houlsby (2003). However, a single viscosity can be sufficient to represent the essential characteristics of the material's behaviour. In the context of the creep behaviour of clay and peat, Ziegler (1972, 1981) presents an argument that aligns with the observations and statements made by Barden (1968), who advocated for simplifying a series of viscous elements with a single element having minimal viscous parameters.

The work presented in this paper builds upon the prior research (Dadras-Ajirlou et al., 2023a, 2023b) and specifically addresses the non-linear, smooth (gradual stiffness degradation), and viscous behaviour of clay with recent loading history effect. To achieve this, the novel concept of the bounding surface in hyperplasticity (Houlsby & Richards, 2023), which differs fundamentally from the bounding surface in plasticity (e.g. Dafalias & Popov, 1975; Dafalias, 1986; Hashiguchi, 1989), is utilised for the first time in viscoplastic modelling using the isotache viscosity. This treatment offers a minimum number of additional material parameters with only a single viscous parameter, in line with the argument by Ziegler (1972, 1981) and Barden (1968). Specifically, the classical force, flow and free energy potentials of the basic MCC hyper-viscoplastic model (Dadras-Ajirlou et al., 2023a) are enriched with additional internal variables and mechanisms using the bounding surface concept. Subsequently, the incremental formulation, including inelastic flow and kinematic hardening, are derived. To incorporate the spacing ratio and exert control over the distribution of deviatoric isotaches in the contractant and the dilatant zones, a more generalised form for force/flow potential is developed. This generalisation is designed to ensure uniqueness of the CS envelope under various loading rates. The paper concludes with an evaluation of the efficacy of the proposed hyper-viscoplastic model in capturing the monotonic and hysteresis time-dependent behaviour of HKMD (Yin & Zhu, 1999; Zhu, 2000; Yin et al., 2002) and the recent loading history effect on the time-dependent behaviour of a saturated clay taken from an earth dam core (Hicher, 2016).

For clarity, the model is constructed under the same notation as the previous work for an axisymmetric system (conventional triaxial specimen) confined to isothermal and infinitesimal strain conditions. Therefore, the total strain can be expressed as sums of volumetric (ε_v = ε_a + 2ε_r) and deviatoric (ε_s = 2(ε_a − ε_r)/3) strains. The corresponding work-conjugate stresses are then the spherical (p = (σ_a + 2σ_r)/3) and the deviatoric (q = σ_a − σ_r) stresses. The subscripts ‘a’ and ‘r’ refer to the axial and radial components of stress and strain. Positive values are assigned to compressive strains and stresses. In the following, the notation (:) denotes the time derivative, all stresses are effective stresses, and inelastic volumetric (ε^p_v) and deviatoric (ε^p_s) strains are adopted as the internal variable of the system. The available space does not permit an outline of the hyperplasticity framework. For a detailed description of the framework, its thermodynamic rudiments, definition of the technical terminologies and the systematic procedure for derivation of components of a constitutive model from different potentials, such as force, flow and free energy potentials, reference is made to Houlsby & Puzrin (2006).

FORCE, FLOW AND FREE ENERGY POTENTIALS

The classical force potential (z) or the dissipation function (d) has already been constructed based on the isotache scaling (Dadras-Ajirlou et al., 2023a):

z = \frac{d}{n} = \frac{r p_{0}}{n} {(\frac{{\dot{ε}}_{v}^{p} + \sqrt{{({\dot{ε}}_{v}^{p})}^{2} + {(M {\dot{ε}}_{s}^{p})}^{2}}}{2 r})}^{n}

1

where M is the frictional coefficient at the CS $({\dot{ε}}_{v}^{p} = 0)$ ⁠, and p₀ is the isotropic pre-consolidation pressure with the definition:

p_{0} = p_{ref} \exp (\frac{ε_{v}^{p}}{λ - κ})

2

where p_ref is the value of p₀ at ε^p_v = 0, and λ and κ are the typical CSSM compressional and swelling indexes on the bi-logarithmic plane of the specific volume against the spherical effective stress.

The isotache scaling in equation (1) is considered by scaling a reference creep power (rp₀) with the ratio of current isotache per the reference isotache (reference strain rate r) to the power of rate sensitivity parameter n. For the current isotache, a mechanism without inelastic swelling has been utilised to impose the unidirectional nature of creep in soil, which as a rheological phenomenon is a compressive and fully dissipative process under a progressive and one-way motion. Note that due to the pure scaling operation in equation (1), energy always dissipates, regardless of the current state of material, or absence of the external power, leading to unending creep or relaxation. However, due to the exponential form of equation (2) and appropriate values for n, this is not a practical issue (Grimstad et al., 2021; Dadras-Ajirlou et al., 2023a). Nevertheless, one can simply add an extra term similar to equation (1), but with n = 1 similar to the example provided by Grimstad et al. (2021) to limit the dissipation by introducing a yielding criterion. For simplicity and practicality, this limit is not implemented in the current study.

The force potential (or dissipation function) described in equation (1), featuring strict convexity and a positive homogeneity degree of n ≥ 1, characterises path-dependent/irreversible behaviour due to its positive semi-definite nature. To define material behaviour fully, path-independent/reversible behaviour must also be established. To do so, the Helmholtz free energy potential (f) proposed by Houlsby et al. (2005) has been utilised:

f = κ p_{a} \exp (\frac{ε_{v} - ε_{v}^{p} + (3 / 2) g {(ε_{s} - ε_{s}^{p})}^{2}}{κ})

3

where g denotes dimensionless material parameters defining the elastic shear modulus, and p_a is an arbitrary reference pressure (preferably atmospheric pressure). This strictly convex Helmholtz free energy potential is particularly suitable for modelling the time-dependent behaviour of clay, such as creep or stress relaxation. Its exponential form, which couples volumetric and shear strains, provides a reference-independent and pressure-dependent state for the material (Dadras-Ajirlou et al., 2023a). These features enable the model, with the laboratory-estimated viscous and reversible material parameters, to reasonably represent the material behaviour under the in situ condition, with a significantly different timescale and pressure. When extending the model with the bounding surface concept, these attributes should be incorporated in the plastic part of free energy, which also influences the path-independent/reversible behaviour.

The above-mentioned features also apply to the Gibbs free energy potential, the Legendre transform of the Helmholtz free energy defined in equation (3). For a detailed description of the role of the Legendre transformation, refer to the appendix in Collins & Houlsby (1997) and Houlsby & Puzrin (2006). The choice between the Helmholtz and Gibbs forms of the free energy potential is influenced principally by modelling preferences. In this context, the main justification for choosing the Helmholtz form is not only to provide an aesthetically pleasing and consistent strain-based description of the system under isothermal processes but also to simplify the construction of the plastic part of free energy as a scalar-valued function of plastic strains. As presented below, the simplest form for the plastic free energy with the same attributes as those of the elastic free energy presented in equation (3), is an exponential function of plastic volumetric and deviatoric strains.

Now, by considering the earlier construction as the bounding condition (equations (1)–(3)), N − 1 additional mechanisms (including inelastic swelling) can be incorporated by increasing the number of internal variables using the bounding surface concept (Houlsby & Richards, 2023):

z = \frac{d}{n} = \frac{r p_{0}}{n} {\{\frac{{\dot{ε}}_{v N}^{p} + \sqrt{\sum_{i = 1}^{N} [{(K_{i} {\dot{ε}}_{v i}^{p})}^{2} + {(M K_{i} {\dot{ε}}_{s i}^{p})}^{2}]}}{2 r}\}}^{n}

4

p_{0} = p_{ref} \exp (\frac{\sum_{i = 1}^{N} ε_{v i}^{p}}{λ - κ})

5

\begin{aligned} f = & κ p_{a} \exp (\frac{ε_{v} - \sum_{i = 1}^{N} ε_{v i}^{p} + (3 / 2) g {(ε_{s} - \sum_{i = 1}^{N} ε_{s i}^{p})}^{2}}{κ}) \\ + \underset{f_{p}}{\underset{⏟}{\sum_{i = 1}^{N} (p_{a} / k_{p}) \exp \{k_{p} [H_{i} ε_{v i}^{p} + (3 / 2) g_{p} {(H_{i} ε_{s i}^{p})}^{2}]\}}} \end{aligned}

6

All the parameters are the same as presented previously for equations (1) and (3), except the dimensionless parameters k_p, g_p, and the weight functions K_i and H_i, chosen to be:

K_{i} = \frac{i}{N}

7

H_{i} = 1 - \frac{i}{N}

8

The parameters k_p, g_p and weight function H_i control the storage and release of reversible inelastic energy $({\dot{f}}_{p})$ ⁠. Notice how the bounding condition is considered at i = N. The force potential (equation (4)) is still positive semi-definite with a positive homogeneity degree of n ≥ 1 in all variables. Moreover, the plastic part of the Helmholtz free energy (f_p) is structured similarly to the elastic part to maintain the pressure dependence and the reference independence attributes, as discussed earlier. The above arrangement preserves the isotache scaling properties. Further elaborations are provided in the following and the proceeding section.

To understand the implications of the recent extension, an alternative expression in stress space is presented. To do so, by following the systematic procedures, the flow potential (w) is first computed in terms of the spherical and deviatoric dissipative stresses. Subsequently, w is transformed to the p–q stress space using Ziegler's theory of thermodynamic orthogonality (hereinafter called ZO condition). Based on the Legendre transformation, w can be expressed in terms of homogeneous force potential z with homogeneity order of n as (Houlsby & Puzrin, 2006):

w = d - z = (n - 1) z

9

In passing, it should also be realised that w = 0 for n = 1 (homogeneity degree of one for d or z, z = d), which signifies the rate-independent bounding surface form for the hyperplastic MCC model (Houlsby, 1981; Collins & Houlsby, 1997). Table 1 summarises the relationships between the systematic extensions of the MCC model to the viscoplastic and bounding surface forms.

Table 1.

Summary of potential functions defining different versions of the MCC model

Model version	Free energy potential	Force potential
Hyper-plastic MCC (equivalent to viscoplastic with n = 1)	$f = κ p_{a} \exp [\frac{ε_{v} - ε_{v}^{p} + (3 / 2) g {(ε_{s} - ε_{s}^{p})}^{2}}{κ}]$	$z = p_{0} [\frac{{\dot{ε}}_{v}^{p} + \sqrt{{({\dot{ε}}_{v}^{p})}^{2} + {(M {\dot{ε}}_{s}^{p})}^{2}}}{2}]$
Hyper-viscoplastic MCC		$z = \frac{r p_{0}}{n} {[\frac{{\dot{ε}}_{v}^{p} + \sqrt{{({\dot{ε}}_{v}^{p})}^{2} + {(M {\dot{ε}}_{s}^{p})}^{2}}}{2 r}]}^{n}$
Hyper-plastic bounding surface MCC (equivalent to viscoplastic with n = 1)	$\begin{aligned} f = & κ p_{a} \exp [\frac{ε_{v} - \sum_{i = 1}^{N} ε_{v i}^{p} + (3 / 2) g {(ε_{s} - \sum_{i = 1}^{N} ε_{s i}^{p})}^{2}}{κ}] \\ + \sum_{i = 1}^{N} (\frac{p_{a}}{k_{p}}) \exp \{k_{p} [H_{i} ε_{v i}^{p} + \frac{3}{2} g_{p} {(H_{i} ε_{s i}^{p})}^{2}]\} \end{aligned},$ $H_{i} = 1 - \frac{i}{N}$	$z = p_{0} \{\frac{{\dot{ε}}_{v N}^{p} + \sqrt{\sum_{i = 1}^{N} [{(K_{i} {\dot{ε}}_{v i}^{p})}^{2} + {(M K_{i} {\dot{ε}}_{s i}^{p})}^{2}]}}{2}\}$ ⁠, $K_{i} = \frac{i}{N}$
Hyper-viscoplastic bounding surface MCC		$z = \frac{r p_{0}}{n} {\{\frac{{\dot{ε}}_{v N}^{p} + \sqrt{\sum_{i = 1}^{N} [{(K_{i} {\dot{ε}}_{v i}^{p})}^{2} + {(M K_{i} {\dot{ε}}_{s i}^{p})}^{2}]}}{2 r}\}}^{n},$ $K_{i} = \frac{i}{N}$

Model version	Free energy potential	Force potential
Hyper-plastic MCC (equivalent to viscoplastic with n = 1)	$f = κ p_{a} \exp [\frac{ε_{v} - ε_{v}^{p} + (3 / 2) g {(ε_{s} - ε_{s}^{p})}^{2}}{κ}]$	$z = p_{0} [\frac{{\dot{ε}}_{v}^{p} + \sqrt{{({\dot{ε}}_{v}^{p})}^{2} + {(M {\dot{ε}}_{s}^{p})}^{2}}}{2}]$
Hyper-viscoplastic MCC		$z = \frac{r p_{0}}{n} {[\frac{{\dot{ε}}_{v}^{p} + \sqrt{{({\dot{ε}}_{v}^{p})}^{2} + {(M {\dot{ε}}_{s}^{p})}^{2}}}{2 r}]}^{n}$
Hyper-plastic bounding surface MCC (equivalent to viscoplastic with n = 1)	$\begin{aligned} f = & κ p_{a} \exp [\frac{ε_{v} - \sum_{i = 1}^{N} ε_{v i}^{p} + (3 / 2) g {(ε_{s} - \sum_{i = 1}^{N} ε_{s i}^{p})}^{2}}{κ}] \\ + \sum_{i = 1}^{N} (\frac{p_{a}}{k_{p}}) \exp \{k_{p} [H_{i} ε_{v i}^{p} + \frac{3}{2} g_{p} {(H_{i} ε_{s i}^{p})}^{2}]\} \end{aligned},$ $H_{i} = 1 - \frac{i}{N}$	$z = p_{0} \{\frac{{\dot{ε}}_{v N}^{p} + \sqrt{\sum_{i = 1}^{N} [{(K_{i} {\dot{ε}}_{v i}^{p})}^{2} + {(M K_{i} {\dot{ε}}_{s i}^{p})}^{2}]}}{2}\}$ , $K_{i} = \frac{i}{N}$
Hyper-viscoplastic bounding surface MCC		$z = \frac{r p_{0}}{n} {\{\frac{{\dot{ε}}_{v N}^{p} + \sqrt{\sum_{i = 1}^{N} [{(K_{i} {\dot{ε}}_{v i}^{p})}^{2} + {(M K_{i} {\dot{ε}}_{s i}^{p})}^{2}]}}{2 r}\}}^{n},$ $K_{i} = \frac{i}{N}$

The generalised dissipative stresses (χ_i) are obtained from the derivation of the force potential with respect to the corresponding work-conjugate variables (inelastic strain rates) as

\begin{aligned} χ_{p i} = \frac{\partial z}{\partial {\dot{ε}}_{v i}^{p}} = (\frac{p_{0}}{2}) {(\frac{n z}{r p_{0}})}^{(n - 1) / n} \\ \times \{⟨i - (N - 1)⟩ + \frac{K_{i}^{2} {\dot{ε}}_{v i}^{p}}{\sqrt{\sum_{i = 1}^{N} [{(K_{i} {\dot{ε}}_{v i}^{p})}^{2} + {(M K_{i} {\dot{ε}}_{s i}^{p})}^{2}]}}\} \end{aligned}

10

\begin{aligned} χ_{q i} = \frac{\partial z}{\partial {\dot{ε}}_{s i}^{p}} = (\frac{p_{0}}{2}) {(\frac{n z}{r p_{0}})}^{(n - 1) / n} \\ \times \{\frac{{(M K_{i})}^{2} {\dot{ε}}_{s i}^{p}}{\sqrt{\sum_{i = 1}^{N} [{(K_{i} {\dot{ε}}_{v i}^{p})}^{2} + {(M K_{i} {\dot{ε}}_{s i}^{p})}^{2}]}}\} \end{aligned}

11

where 〈·〉 in equation (10) is the Macaulay bracket. Now by combining equations (10) and (11), and applying the relation between the homogeneous force and flow potentials (equation (9)), w can be found after moderate algebra as

\begin{aligned} w = & r p_{0} (\frac{n - 1}{n}) {(\frac{χ_{p N}}{p_{0}} \{\frac{\sum_{i = 1}^{N} [{(M χ_{p i} / K_{i})}^{2} + {(χ_{q i} / K_{i})}^{2}]}{{(M χ_{p N})}^{2}}\})}^{n / (n - 1)} \\ = & r p_{0} (\frac{n - 1}{n}) {(\frac{p_{eq}}{p_{0}})}^{n / (n - 1)} \end{aligned}

12

where p_eq is called equivalent pressure. Note that the flow potential w exhibits a positive homogeneity degree of n/(n − 1) with respect to the dissipative stresses. This is expected based on the scaling property of Legendre-conjugate homogeneous functions (1/n + (n − 1)/n = 1).

To apply the ZO condition and complete the transformation, p, q and the generalised conservative stresses $(\bar{χ_{i}})$ are computed by taking the derivatives of the Helmholtz free energy (equation (6)) with respect to the corresponding conjugate variables:

\begin{aligned} p = \frac{\partial f}{\partial ε_{v}} = p_{a} \exp [\frac{ε_{v} - \sum_{i = 1}^{N} ε_{v i}^{p} + \frac{3}{2} g {(ε_{s} - \sum_{i = 1}^{N} ε_{s i}^{p})}^{2}}{κ}] \end{aligned}

13

q = \frac{\partial f}{\partial ε_{s}} = 3 g p (ε_{s} - \sum_{i = 1}^{N} ε_{s i}^{p})

14

\begin{aligned} \bar{χ_{p i}} = & - \frac{\partial f}{\partial ε_{v i}^{p}} = p - p_{a} H_{i} \exp \{k_{p} [H_{i} ε_{v i}^{p} + \frac{3}{2} g_{p} {(H_{i} ε_{s i}^{p})}^{2}]\} \\ = & p - p_{b i} \end{aligned}

15

\bar{χ_{q i}} = - \frac{\partial f}{\partial ε_{s i}^{p}} = q - 3 g p_{b i} H_{i} ε_{s i}^{p} = q - q_{b i}

16

where p_bi and q_bi are the spheric and deviatoric back-stresses. Notice that p_bN = q_bN = 0 since H_N = 0. Now by imposing ZO condition:

\{\begin{matrix} \bar{χ_{p i}} = χ_{p i} \\ \bar{χ_{q i}} = χ_{q i} \end{matrix}

17

the equivalent pressure p_eq in equation (12) can be expressed in the p–q space:

\begin{aligned} p_{eq} & = p (\frac{\sum_{i = 1}^{N} \{{[M (p - p_{b i}) / K_{i}]}^{2} + {[(q - q_{b i}) / K_{i}]}^{2}\}}{M^{2} p^{2}}) \\ = p (1 + \frac{η^{2}}{M^{2}}) \end{aligned}

18

where, by considering that p_bN = q_bN = 0 (since H_N = 0), η, the stress ratio, can be defined as:

η = \frac{\sqrt{q^{2} + \sum_{i = 1}^{N - 1} \{{[M (p - p_{b i}) / K_{i}]}^{2} + {[(q - q_{b i}) / K_{i}]}^{2}\}}}{p}

19

It can now be reasoned that under the condition of ZO, the pressure dependency emerges within all formulations. The pressure dependency feature is due to the exponential coupling of the volumetric and deviatoric strains in both elastic and inelastic parts of the Helmholtz free energy (equation (6)), as well as the imposition of the unidirectional nature of time-dependent behaviour of clay (in a phenomenological sense based on the isotache concept) within the force potential (equation (4)). The latest criterion is established through two choices made at the bounding condition (i = N), one in the force potential (z) and the other in the inelastic part of the Helmholtz free energy (f_p). The first choice enforces zero dissipation for the pure plastic swelling at the bounding condition $({\dot{ε}}_{v N}^{p} + |{\dot{ε}}_{v N}^{p}| = 0)$ ⁠, while the second ensures zero back-stress at the bounding condition (H_N = 0). Based on these choices and the exponential form of the free energy, pressure remains strictly positive, which is required by the isotache scaling (Dadras-Ajirlou et al., 2023a).

INCREMENTAL FORMULATION

The incremental formulation can be derived directly from two potential functions, the free energy potential representing reversible response, and the force/dissipation potential representing irreversible response. However, to facilitate the implementation of the model into numerical codes by employing the common numerical integration schemes, the flow potential can be used instead of the force potential to derive the irreversible incremental response completely in terms of stress.

Based on equations (13) and (14) the following time differentials can be derived to describe the evolution of p and q stresses over time:

\begin{aligned} \dot{p} = \frac{\partial^{2} f}{\partial ε_{v}^{2}} {\dot{ε}}_{v} + \frac{\partial^{2} f}{\partial ε_{v} \partial ε_{s}} {\dot{ε}}_{s} + \sum_{i = 1}^{N} (\frac{\partial^{2} f}{\partial ε_{v} \partial ε_{v i}^{p}} {\dot{ε}}_{v i}^{p} + \frac{\partial^{2} f}{\partial ε_{v} \partial ε_{s i}^{p}} {\dot{ε}}_{s i}^{p}) \end{aligned}

20

\begin{aligned} \dot{q} = \frac{\partial^{2} f}{\partial ε_{s} \partial ε_{v}} {\dot{ε}}_{v} + \frac{\partial^{2} f}{\partial ε_{s}^{2}} {\dot{ε}}_{s} + \sum_{i = 1}^{N} (\frac{\partial^{2} f}{\partial ε_{s} \partial ε_{v i}^{p}} {\dot{ε}}_{v i}^{p} + \frac{\partial^{2} f}{\partial ε_{s} \partial ε_{s i}^{p}} {\dot{ε}}_{s i}^{p}) \end{aligned}

21

Now since:

\frac{\partial^{2} f}{\partial ε_{v}^{2}} = - \frac{\partial^{2} f}{\partial ε_{v} \partial ε_{v i}^{p}}

22

\frac{\partial^{2} f}{\partial ε_{s}^{2}} = - \frac{\partial^{2} f}{\partial ε_{s} \partial ε_{s i}^{p}}

23

\frac{\partial^{2} f}{\partial ε_{v} \partial ε_{s}} = \frac{\partial^{2} f}{\partial ε_{s} \partial ε_{v}} = - \frac{\partial^{2} f}{\partial ε_{s} \partial ε_{v i}^{p}} = - \frac{\partial^{2} f}{\partial ε_{v} \partial ε_{s i}^{p}}

24

the above system of differentials can be simplified as

\begin{aligned} [\begin{matrix} \dot{p} \\ \dot{q} \end{matrix}] = & \underset{D^{e}}{\underset{⏟}{[\begin{matrix} \frac{\partial^{2} f}{\partial ε_{v}^{2}} & \frac{\partial^{2} f}{\partial ε_{v} \partial ε_{s}} \\ \frac{\partial^{2} f}{\partial ε_{v} \partial ε_{s}} & \frac{\partial^{2} f}{\partial ε_{s}^{2}} \end{matrix}]}} [\begin{matrix} {\dot{ε}}_{v} - \sum_{i = 1}^{N} {\dot{ε}}_{v i}^{p} \\ {\dot{ε}}_{s} - \sum_{i = 1}^{N} {\dot{ε}}_{s i}^{p} \end{matrix}] \\ = & [\begin{matrix} \frac{p}{κ} & \frac{q}{κ} \\ \frac{q}{κ} & \frac{p}{κ} (3 g κ + \frac{q^{2}}{p^{2}}) \end{matrix}] [\begin{matrix} {\dot{ε}}_{v} - \sum_{i = 1}^{N} {\dot{ε}}_{v i}^{p} \\ {\dot{ε}}_{s} - \sum_{i = 1}^{N} {\dot{ε}}_{s i}^{p} \end{matrix}] \end{aligned}

25

where $D^{e}$ is the elasticity tensor whose components are expressed in terms of p and q based on the definitions provided in equations (13) and (14).

The viscoplastic strain rates, based on the definition (Houlsby & Puzrin, 2006), can be directly derived from the flow potential w (equation (12)). In this regard for the bounding viscoplastic flow, there are:

{\dot{ε}}_{v N}^{p} = \frac{\partial w}{\partial χ_{p N}} = r {(\frac{p_{eq}}{p_{0}})}^{1 / (n - 1)} [1 - {(\frac{η}{M})}^{2}]

26

{\dot{ε}}_{s N}^{p} = \frac{\partial w}{\partial χ_{q N}} = r {(\frac{p_{eq}}{p_{0}})}^{1 / (n - 1)} (\frac{2 q}{M^{2} p})

27

and for the other viscoplastic flows:

{\dot{ε}}_{v i}^{p} = \frac{\partial w}{\partial χ_{p i}} = r {(\frac{p_{eq}}{p_{0}})}^{1 / (n - 1)} [\frac{2 (p - p_{b i})}{K_{i}^{2} p}]

28

{\dot{ε}}_{s i}^{p} = \frac{\partial w}{\partial χ_{q i}} = r {(\frac{p_{eq}}{p_{0}})}^{1 / (n - 1)} [\frac{2 (q - q_{b i})}{{(M K_{i})}^{2} p}]

29

In the above equations, the ZO condition (equation (17)) has been applied after the derivations of w with respect to the corresponding conjugate dissipative stresses. η in equation (26) refers to the definition provided in equation (19) for the stress ratio in the p–q stress space.

There are also two hardening rules controlling the incremental response of the material. The isotropic hardening rule refers to the evolution of p₀ during an inelastic process, which according to equation (5) can be expressed as

{\dot{p}}_{0} = p_{0} (\frac{\sum_{i = 1}^{N} {\dot{ε}}_{v i}^{p}}{λ - κ})

30

The kinematic hardening rule can be expressed as a time differential of the back stresses p_bi and q_bi defined in equations (15) and (16), respectively:

\begin{aligned} [\begin{matrix} {\dot{p}}_{b i} \\ {\dot{q}}_{b i} \end{matrix}] & = [\begin{matrix} - \frac{\partial^{2} f^{p}}{\partial ε_{v i}^{2}} & - \frac{\partial^{2} f^{p}}{\partial ε_{v i} \partial ε_{s i}} \\ - \frac{\partial^{2} f^{p}}{\partial ε_{v i} \partial ε_{s i}} & - \frac{\partial^{2} f^{p}}{\partial ε_{s i}^{2}} \end{matrix}] [\begin{matrix} {\dot{ε}}_{v i}^{p} \\ {\dot{ε}}_{s i}^{p} \end{matrix}] \\ = H_{i} k_{p} [\begin{matrix} p_{b i} & q_{b i} \\ q_{b i} & p_{b i} (\frac{3 g_{p}}{k_{p}} + \frac{q_{b i}^{2}}{p_{b i}^{2}}) \end{matrix}] [\begin{matrix} {\dot{ε}}_{v i}^{p} \\ {\dot{ε}}_{s i}^{p} \end{matrix}] \end{aligned}

31

in which the same properties as those expressed in equations (22)–(24) have been applied for f^p. Upon replacing the viscoplastic strain rate with their definitions in equations (28) and (29) in which the ZO condition is imposed, the kinematic hardening rule transforms to:

\begin{aligned} [\begin{matrix} {\dot{p}}_{b i} \\ {\dot{q}}_{b i} \end{matrix}] = & \frac{2 H_{i} k_{p} r}{K_{i} p} {(\frac{p_{eq}}{p_{0}})}^{1 / (n - 1)} \\ [\begin{matrix} p_{b i} & q_{b i} \\ q_{b i} & p_{b i} (\frac{3 g_{p}}{k_{p}} + \frac{q_{b i}^{2}}{p_{b i}^{2}}) \end{matrix}] [\begin{matrix} p - p_{b i} \\ \frac{q - q_{b i}}{M^{2}} \end{matrix}] \end{aligned}

32

Equation (32) is a system of the logistic differential equations in terms of back-stresses, meaning that the p–q stress state acts as an attractor for the back-stresses whose evolution over time is exponential. As a result, after sufficient deformation or elapsed time during, for example, the creep/relaxation process, the short-term memory of material in terms of the relative location of back-stresses with the current stress state in stress space diminishes. At such bounding states (states with faded short-term memory), like the CS or the creep/relaxation after a long time, the response of the bounding surface model is quite similar to the previous version of the model (Dadras-Ajirlou et al., 2023a). This is because, based on equation (32), the back-stresses asymptotically approach the current stress state (the stable steady state of the logistic differential equation (32)) during sustained monotonic loading, resulting in both infinitesimal inelastic strain rates (equations (28) and (29)) and dissipation shares for the mechanisms other than the bounding mechanism. Consequently, the properties related to these bounding states are the same as in the previous version of the model. In this regard, under a pure creep process, isotropic hardening, equation (30), can be understood as the evolution of the material's resistance to further compression over a relatively long time – that is, the long-term memory of the material.

It is now evident that by implementing the choices made in the previous section under the condition of ZO and arriving at the viscoplastic flows (equations (26)–(29)) and the kinematic hardening rule (equation (32)), a unique CS envelope as a bounding condition can emerge. At the CS, where the material continuously deforms under constant stress state (q = Mp), the back-stresses are asymptotically attracted towards the constant stress state. In this case, according to equation (32), p_bi ≃ p and q_bi ≃ q, and subsequently based on equations (28), (29) and (19), ${\dot{ε}}_{v i}^{p} ≃ 0$ ⁠, ${\dot{ε}}_{s i}^{p} ≃ 0 (i \neq N)$ ⁠, and η ≃ M, respectively. As a result, the determining internal variables at the CS are the internal variables associated with the bounding condition (i = N) whose volumetric rate is asymptotically zero, with the deviatoric rate continuously evolving, according to equations (26) and (27), respectively. Indeed, at the CS, the response of the bounding surface model is asymptotically equal to the response of the basic model (Dadras-Ajirlou et al., 2023a).

Similarly, the response of the bounding surface model is asymptotically equal to the response of the basic model under the creep process, meaning that the material properties for the creep process remain consistent with those established previously based on the time resistance concept (Janbu, 1969, 1985; Vermeer & Neher, 1999; Grimstad et al., 2010; Jostad & Yannie, 2017). Further details are provided in the following.

Figure 1 displays the performance of the kinematic hardening rule (equation (32)) during a creep process under a certain constant p–q stress state. The largest surface refers to the bounding condition, and the two small, dashed surfaces denote the two additional mechanisms incorporated into the dissipation or force potential. The corresponding viscoplastic flows are also shown in Fig. 1, where the solid grey vector is associated with the flow at the bounding condition (equations (26) and (27)) and the dashed vectors correspond to the two additional mechanisms (equations (28) and (29)). The solid black surface defined by equation (18) is always convex since it is a conic combination of the bounding and the two additional mechanisms in the forms of the MCC-type ellipsoids. The solid black vectors represent the overall viscoplastic flow as the resultant of all viscoplastic flows shown in grey. As observed, due to creep, the two back-stresses are approaching the current p–q stress state, causing two kinematic surfaces and their corresponding viscoplastic flows to converge, and eventually the corresponding viscoplastic flows diminish. Consequently, the short-term memory of the material fades from Fig. 1(a) to Fig. 1(d).

Fig. 1.

View large Download slide

Attraction of back-stresses by current stress state and erasing of short-term memory during a pure creeping process (a)–(b)–(c)–(d). Back-stresses are shown by plus and asterisk signs and M = 1·2

It is important to note that the kinematic and bounding surfaces depicted in Fig. 1. are merely theoretical representations and are not actual components of the proposed model. The sole surface with a physical significance is the dynamic surface, shown by the solid black curve defined in equation (18), encompassing all the back-stresses. Based on equation (18), the kinematic surfaces can be expressed as

p_{eq, i} = p [{(\frac{p - p_{b i}}{K_{i} p})}^{2} + {(\frac{q - q_{b i}}{M K_{i} p})}^{2}]

33

where p_eq,i is the equivalent pressure or size of the ith dynamic kinematic surface that is a function of distance between the current stress state and the ith back-stress.

In the context of viscoplasticity, the kinematic surfaces shown in Fig. 1 are not necessarily always bounded by the bounding surface. In fact, due to the ZO condition (equation (17)), the bounding surface only bounds the back-stresses after a sufficiently long elapsed time and reduction in rate of change of material state. For instance, as shown in Fig. 2, during reversed loading or stress relaxation immediately after fading of the short-term memory during the recent sustained loading, a large domain of the kinematic surfaces can be located outside the bounding surface (generation of a new short-term memory). This is because the distance between the current stress state and the back-stress state suddenly increases. However, after sufficient lapse of time and decrease in the rate of change of stress, the back-stresses are asymptotically attracted to the current stress state and subsequently to the bounding surface, which always passes through the origin and current stress state (due to the unidirectionality condition imposed for the bounding internal variable in the force potential, as explained at the end of the preceding section).

Fig. 2.

View large Download slide

Generation of new short-term memory due to the reversed loading immediately after the state shown in Fig. 1(d). Note that M = 1·2 and back-stresses shown by plus and asterisk signs are overlapping due to their attraction during the previous creeping process

GENERALISED FORCE AND FLOW POTENTIALS

In this section, the force potential defined in equation (4), and subsequently, the flow potential computed in equation (12) as the conjugate of the force potential, are generalised to incorporate the spacing ratio. This extension offers flexibility to adjust the distribution of deviatoric isotaches and enhances the model's predictions of shearing behaviour under different loading rates for different types of clays.

Continuing in line with previous works (Dadras-Ajirlou et al., 2023a, 2023b), the force potential is generalised by incorporating a transition function (T) and the spacing ratio (R) as

\begin{aligned} z = \frac{d}{n} = \frac{r p_{0}}{n} \times {(\frac{{\dot{ε}}_{v N}^{p} + \sqrt{{(T {\dot{ε}}_{v N}^{p})}^{2} + {(M {\dot{ε}}_{s N}^{p})}^{2} + \sum_{i = 1}^{N - 1} ({(K_{i} {\dot{ε}}_{v i}^{p})}^{2} + {(M K_{i} {\dot{ε}}_{s i}^{p})}^{2})}}{R r})}^{n} \end{aligned}

34

where T is a variant of the logistic function with the state variable S as a varying input:

T = \frac{R}{2} + (\frac{R}{2} - 1) \tanh (S)

35

in which R is greater than one and the state variable S is defined as:

S = {(\frac{M}{η})}^{2} - {(\frac{η}{M})}^{2}

36

where η is defined in equation (19).

S is a state variable since based on equations (13)–(16), S, under ZO condition (equation (17)), is linked to the Helmholtz free energy representing the current state of the material. The state variable S is indeed the distance of the current state of the material from the CS emerging as a bounding condition after sustained continuous deformation (details in the preceding section). The function T defined in equation (34) is called the transition function (Dadras-Ajirlou et al., 2023a) due to the property of hyperbolic tangent function (tanh) that smoothly connects the contractant (M > η) and the dilatant (η > M) zones, separated by the critical state line passing through the desired spacing ratio $(q = M \bar{p_{eq}} / R)$ shown in Fig. 3(b) through the evolution of η.

Fig. 3.

View large Download slide

Effect of the spacing ratio R = 3 on the bounding and kinematic surfaces during a creep process (a)–(b) with (a) short-term memory and (b) faded short-term memory. Note that M = 1·2, the back-stresses are shown by plus and asterisk signs, and ${\bar{p}}_{eq}$ is the equivalent pressure of the bounding surface defined exclusively based on the current stress state (as in Dadras-Ajirlou et al. (2023a)) without any back-stresses

It should also be noted that the force potential z or the dissipation function d in equation (34), which quantifies the amount of dissipated energy, are positive semi-definite for any rates of internal variables, like those defined previously in equations (1) and (4). In this regard, it should be also emphasised that the extreme pure delayed plastic swelling (where p approaches zero) is discarded in all force potentials as a possible dissipative process. With these considerations in mind, by following procedures similar to those in the previous section, the flow potential can be obtained after a moderately tedious algebra as

w = r p_{0} (\frac{n - 1}{n}) {(\frac{p_{eq}}{p_{0}})}^{n / (n - 1)}

37

with the following definition for the equivalent pressure in the dissipative stress space:

p_{eq} = \frac{R χ_{p} \sqrt{(M^{2} - η_{χ}^{2}) + η_{χ}^{2} {(M T + \sqrt{M^{2} + (T^{2} - 1) η_{χ}^{2}})}^{2}}}{T (M^{2} - η_{χ}^{2}) + \sqrt{(M^{2} - η_{χ}^{2}) + η_{χ}^{2} {(M T + \sqrt{M^{2} + (T^{2} - 1) η_{χ}^{2}})}^{2}}}

38

where η_χ is the stress ratio in the dissipative stress space:

η_{χ} = \frac{\sqrt{χ_{q N}^{2} + \sum_{i = 1}^{N - 1} [{(M χ_{p i} / K_{i})}^{2} + {(χ_{q i} / K_{i})}^{2}]}}{χ_{p N}}

39

Notice that T defined in equation (35) is the function of current stress state (through η in equation (36)) causing the inelastic flow to be practically non-associated, as in the basic model (Dadras-Ajirlou et al., 2023a). Fig. 3 demonstrates the effect of spacing ratio R on the relative location of critical state. As observed, the bounding surface (solid grey curve), which controls the CS as a bounding condition, has passed the CS at the desired relative location.

It is essential to highlight that incorporating the spacing ratio in the isotache viscoplasticty using the commonly employed approach of Collins & Hilder (2002) and Collins (2003) in plasticity, where the volumetric and shearing dissipative mechanism is scaled by a function of the pre-consolidation pressure (p₀), results in the loss of uniqueness of the CS envelope under different loading rates. This has been demonstrated by Grimstad et al. (2020), Grimstad et al. (2021) and Dadras-Ajirlou et al. (2023b).

MODEL PARAMETERS

The proposed hyper-viscoplastic model has eight dimensionless material parameters, as detailed in Table 2. In comparison to the previous version of the model (Dadras-Ajirlou et al., 2023a), only two additional parameters are introduced to control the kinematic hardening and capture hysteresis behaviour. These parameters can readily be determined through conventional triaxial, oedometer, or isotropic compression tests.

Table 2.

Dimensionless parameters of the model and their values for HKMD and a compacted clay (CC)

Model parameters	Description	Required test	HKMD*	CC†
Calibrated parameters for N = 10 (ten internal variables)
κ	Swelling index	IC or oedometer tests with URC	0·0102	0·01
g	Shear modulus coefficient	T(U) or oedometer tests with URC	125	100
k_p	Kinematic hardening coefficient	T(U) or oedometer tests with URC	4000	700
g_p	Kinematic hardening coefficient	T(U) or oedometer tests with URC	900	350
R	Spacing ratio	T(U)	2·175	2·1
Laboratory estimated parameters
μ	Creep index	IC or oedometer creep tests	0·0025	0·003
λ	Compression index	IC or oedometer tests	0·0792	0·092
M	Slope of critical state line in p–q stress space	T(U)	1·265	1·07

Model parameters	Description	Required test	HKMD*	CC†
Calibrated parameters for N = 10 (ten internal variables)
κ	Swelling index	IC or oedometer tests with URC	0·0102	0·01
g	Shear modulus coefficient	T(U) or oedometer tests with URC	125	100
k_p	Kinematic hardening coefficient	T(U) or oedometer tests with URC	4000	700
g_p	Kinematic hardening coefficient	T(U) or oedometer tests with URC	900	350
R	Spacing ratio	T(U)	2·175	2·1
Laboratory estimated parameters
μ	Creep index	IC or oedometer creep tests	0·0025	0·003
λ	Compression index	IC or oedometer tests	0·0792	0·092
M	Slope of critical state line in p–q stress space	T(U)	1·265	1·07

*The measured parameters are taken from original reports (Yin & Zhu, 1999; Zhu, 2000; Yin et al., 2002).

†The measured parameters are taken from Shi et al. (2023).

Note: IC, isotropic consolidation; T(U), triaxial (undrained preferably); URC, unloading–reloading cycle.

The rate sensitivity parameter (n) and a reference strain rate (reference isotache r), as mentioned previously, control the isotache scaling of a reference power (rp₀). Following the approach of Grimstad et al. (2010) and Jostad & Yannie (2017), among others, n is defined as

n = 1 + \frac{μ}{λ - κ}

40

where μ is the creep index. Based on the previous discussion, since the pure creep process occurs exclusively at the bounding condition, μ should be evaluated from a creep test. In this regard, Jostad & Yannie (2017), among others, presented a clear procedure for estimation of μ based on the time resistance concept in which 1/μ is the slope of creep response in terms of strain rate against the elapsed time.

r is the average volumetric strain rate under the creep process at a certain stress ratio. It can be estimated from conventional 24 h incremental loading oedometer or isotropic compression tests. For instance, the creep strain rate under the sustained K₀ stress state in the oedometer test, where the short-term memory has been faded, can be computed as

{{\dot{ε}}_{v}^{p}|}_{oed} = {{\dot{ε}}_{v N}^{p}|}_{K_{0}} = {\frac{\partial w}{\partial χ_{p N}}|}_{K_{0}} = r {(\frac{p_{eq}}{p_{0}})}^{1 / (n - 1)} {\frac{\partial p_{eq}}{\partial χ_{p N}}|}_{K_{0}}

41

The creep strain rate under the K₀ stress state, based on Bjerrum's equivalent time concept (Bjerrum, 1967), can be defined as (Jostad & Yannie, 2017):

{{\dot{ε}}_{v}^{p}|}_{oed} = \frac{μ}{τ}

42

where τ is the equivalent time normally equal to 24 h. Now by assuming p_eq = p₀ after 24 h K₀ compression, the reference isotache (r) can be computed by combining equations (41) and (42) as

r = \frac{μ}{τ {(\partial p_{eq} / \partial χ_{p N})|}_{K_{0}}}

43

Note that in equation (43), the derivative is first taken with respect to the spherical dissipative stress associated with the bounding condition (χ_pN). Then the K₀ stress state is applied by imposing the ZO condition:

η_{K_{0}} = \frac{χ_{q N}}{χ_{p N}} = \frac{q}{p} = \frac{3 (1 - K_{0})}{1 + 2 K_{0}}

44

where K₀ is the radial to the axial stress ratio under sustained oedometric loading estimated by the model, not the oedometer test. To incorporate a material representative value for K₀, the proposed model can be enriched with an extra degree of freedom as practised by Collins & Hilder (2002) and Rollo & Amorosi (2020), among others.

k_p and g_p, as mentioned earlier, control the performance of the kinematic hardening rule. Since kinematic hardening rule plays a crucial role in determining the reversible response through the release and storage of inelastic free energy (hysteresis), k_p and g_p should be calibrated alongside their counterparts, κ and g, which govern the release and storage of the elastic free energy. In the best scenario, k_p and κ can be first calibrated based on an isotropic hysteresis response, and then g_p and g are calibrated for a shear hysteresis response. This calibration process has been employed for the HKMD since both isotropic and shear hysteresis responses were available. It 's also important to highlight that, in the bounding surface form, unlike the multi-surface form (e.g. Puzrin & Houlsby, 2001, 2003; Einav & Puzrin, 2004; Apriadi et al., 2013), all dissipative mechanisms are coupled. They are neither mutually exclusive nor sequential. As a result, the considered number of internal variables significantly influences the values of k_p and g_p. Based on the authors’ experience, between five and ten numbers of internal variables are sufficient for effectively capturing non-linear and hysteresis behaviour, even for small loading cycles. Potentially, even smaller numbers of internal variables may be sufficient if kinematic hardening moduli associated with a certain internal variable (or back-stress) are discretely calibrated (no use of weight functions), as shown by Houlsby & Richards (2023). This will remain to be explored further in future endeavours.

EVALUATION OF MODEL

In this section, the efficacy of the proposed model in simulating the response of the reconstituted HKMD (Yin & Zhu, 1999; Zhu, 2000; Yin et al., 2002) and a saturated compacted clay (Hicher, 2016) under different loading conditions is explored. The model parameters for these clayey soils are presented in Table 2.

The first test to consider is the 24 h isotropic consolidation test on HKMD, which includes an unloading–reloading cycle. Based on this test, which is completely uncoupled from the shear behaviour, the parameters k_p and κ are calibrated as those presented in Table 2. As shown in Fig. 4, the model response can reasonably capture the discrete experimental response, including the unloading–reloading cycle.

Fig. 4.

View large Download slide

Comparison between the experimental (Yin et al., 2002) and the simulated 24 h isotropic loading and unloading test on reconstituted HKMD

Next, to calibrate the parameters g_p, g and R, a shearing test with loading and unloading cycles is simulated. The available test for HKMD is an undrained triaxial test with complex loading stages, as detailed in Table 3. Fig. 5 shows the comparison with the test data and their simulation using the calibrated parameters. Despite minor differences from the lab data, the overall performance of the model is remarkable, especially for the unloading–reloading responses, which are challenging to capture.

Fig. 5.

View large Download slide

Comparison between the experimental (Yin et al., 2002) and the simulated results of multi-stage undrained triaxial compression test on reconstituted and normally consolidated HKMD in terms of: (a) deviatoric stress plotted against axial strain; (b) excess pore water pressure plotted against axial strain; (c) deviatoric stress plotted against mean effective stress

Table 3.

Loading history of the multi-stage triaxial compression test on HKMD (Yin et al., 2002)

Schedule	Loading	Unloading	Reloading	Relaxation	Loading	Relaxation	Loading	Relaxation
Axial strain rate: 1/min	0·1%	−0·1%	0·1%	0	0·01%	0	0·001%	0
Duration: min	29	7	20	2540	232	1320	830	705

Now, to evaluate the performance of the proposed model across different stress levels, the responses of a set of undrained triaxial tests with different initial conditions are predicted. These tests were conducted on samples with different overconsolidation ratios (OCR = p_max/p). Following the original report (Zhu, 2000), to construct the history of the material and achieve the desired OCRs before the undrained triaxial shearing, the isotropic loading and unloading stages, both under 36 h, are performed first. Fig. 6 displays the comparison between the measured and predicted responses. As can be observed, the model can function under different initial stress states before shearing, thanks to the logistic differential form of the kinematic hardening rule (equation (32)). Considering the challenge of capturing the behaviour of HKMD, as shown in earlier attempts (Yin et al., 2002; Bodas Freitas et al., 2011; Yang et al., 2016; Islam & Gnanendran, 2017; Shahbodagh et al., 2020; Mánica et al., 2021; Yao et al., 2015; Qiao et al., 2016; Bathayian & Maleki, 2023), the overall performance of the model in both rate-dependent monotonic and hysteresis behaviour of HKMD is reasonably good.

Fig. 6.

View large Download slide

Comparison between the experimental (Yin et al., 2002) and the simulated undrained triaxial compression tests on reconstituted HKMD with different overconsolidation ratios (OCRs) under constant axial strain rate of 1·5%/h in terms of (a) the stress–strain response and (b) the normalised effective stress path. Note that p_max is the maximum isotropic consolidation pressure, the initial stress for OCR = 1 is p = 400 kPa, and for the other OCRs is p = 100 kPa

To assess the proposed model's effectiveness in capturing the effects of recent loading, simulations are extended using the tests conducted by Hicher (2016) on a saturated compacted clay from an earth dam core at a depth between 5 and 8 m. This experimental study is unusual and provides insights into the influence of recent loading history on the creep and stress relaxation behaviour of clay. The measurable parameters of the model in Table 2 are selected based on the values reported by Li et al. (2022) and Shi et al. (2023).

The first sets of tests to consider are the undrained triaxial tests under various loading rates. The comparison between the test and simulation results is illustrated in Fig. 7. Fig. 8 shows the comparison for creep test under an axisymmetric undrained condition at different deviatoric stress levels obtained under a strain rate of 1%/min. These results indicate the capability of the model to capture the conventional rate and creep behaviour of the saturated compacted clay based on selected values for parameters, particularly the single viscous parameter (μ).

Fig. 7.

View large Download slide

Comparison between the experimental (Hicher, 2016) and the simulated undrained triaxial compression tests on saturated compacted clay under different strain rates in terms of (a) deviatoric stress plotted against axial strain response and (b) the effective stress path in p–q stress space

Fig. 8.

View large Download slide

Comparison between the experimental (Hicher, 2016) and the simulated undrained triaxial creep tests at different deviatoric stresses obtained under strain rate of 1%/min in terms of (a) axial strain and (b) excess pore water pressure plotted against time

The results of the undrained triaxial test, showing the effects of recent loading history on the time-dependent behaviour, are compared with the simulation results in Fig. 9. This test, conducted at a strain rate of 1%/min, includes two unloading–reloading cycles with nine stress relaxation phases (represented by markers SR 1 to 9) conducted at various stages of the loading, unloading and reloading processes. Notably, after each loading–unloading cycle, the difference between the test results and simulation in the loading processes increases. Interestingly, the test exhibits a stiffer response than the simulation, even for strains below 2%, contrary to the triaxial test under the same conditions in Fig. 7, where the simulation shows a stiffer response, particularly for strains below 2%.

Fig. 9.

View large Download slide

Comparison between the experimental (Hicher, 2016) and the simulated undrained triaxial test under strain rate of 1%/min with two unloading–reloading cycles and nine relaxation phases (SR 1–SR 9) conducted at different stages of the test. A full-colour version of this figure can be found on the ICE Virtual Library (www.icevirtuallibrary.com)

The detailed effect of recent loading history on relaxation behaviour, in terms of deviatoric stress increment, is depicted for each phase in Fig. 10. The relaxation immediately after loading at phases SR 1, SR 3 and SR 7 resulted in a decrease in deviatoric stress in both the test and simulation results. In contrast, immediately after unloading at phases SR 2 and SR 9, deviatoric stress increases in both the test and simulation. Interestingly, at intermediate stages of the unloading and reloading processes, specifically at phases SR 5 and SR 6, the relaxation in both the test and simulation exhibits an opposite response compared to those occurring at SR 3 and SR 2, respectively, which occurred at the initial stages of the unloading and reloading processes. Particularly noteworthy are the relaxation phases SR 4 and SR 8, where the model was less successful in capturing the response. In these phases, which occurred after stress decrease due to relaxation in phases SR 3 and SR 7 followed by slight unloading (0·1% decrease of axial strain), deviatoric stress initially increases and then decreases. The model can capture this particular behaviour to some extent for phase SR 8, but for phase SR 4, it is unsuccessful. However, overall, the model performs reasonably well in capturing hysteresis, rate and time-dependent behaviour with the recent loading history effect, especially when compared to the limited attempts in the literature (Li et al., 2022; Shi et al., 2023) to predict the challenging behaviour of this clay (Hicher, 2016) under complex loading conditions.

Fig. 10.

View large Download slide

Comparison between the experimental (Hicher, 2016) and the simulated nine undrained triaxial relaxation phases (SR 1–SR 9) along different stages of the undrained triaxial test with complex loading schemes presented in Fig. 9

CONCLUSION

This study has extended the previously developed hyper-viscoplastic model (Dadras-Ajirlou et al., 2023a) by incorporating inelastic free energy and the bounding surface concept (Houlsby & Richards, 2023) to capture the non-linear and hysteresis behaviour of clay. The models are mathematically rigorous and are guaranteed to obey the laws of thermodynamics, since they have been built entirely within the hyperplasticity framework (Houlsby & Puzrin, 2006) by specifying two potential functions, one for reversible and the other for irreversible behaviour. In addition, the proposed model with inelastic free energy complies with critical state soil mechanics (CSSM) and the isotache concept. As a result, the main limitation of the basic hyper-viscoplastic model (Dadras-Ajirlou et al., 2023a, 2023b), namely, discarding inelastic free energy to secure a unique critical state (CS) envelope, has been lifted. Another salient feature of the model is the pressure dependency (a common feature of geomaterials behaviour) in both reversible and irreversible behaviour. There are overall eight dimensionless material parameters in the proposed model. Three of these parameters pertain to the bounding conditions of the pure creep and the CS and can be measured in a straightforward way from conventional laboratory tests. While some of the other five parameters can be estimated from laboratory tests, it is recommended that they should be calibrated collectively against the conventional undrained triaxial, oedometer or isotropic compression tests with at least a cycle of unloading and reloading. The efficacy of the model in capturing the hysteresis and the effect of recent loading history on the time-dependent behaviour of clays has been demonstrated by simulating a set of laboratory tests on the reconstituted HKMD (Yin & Zhu, 1999; Zhu, 2000; Yin et al., 2002) and a saturated compacted clay (Hicher, 2016).

Clearly, the ultimate goal of any constitutive model is its application in solving engineering boundary value problems. To this end, the numerical integration of the incremental formulation over strain and time increments is crucial. Based on the simple incremental formulation of the proposed model, the most commonly used algorithms, such as explicit modified-forward Euler with error control (Sloan, 1987) or implicit backward Euler (e.g. Simo & Hughes, 1998; Heeres et al., 2002), can be employed for the numerical integration. In fact, their application in the integration of the proposed model is straightforward because of deliberate omission of the yield criterion for pragmatic reasons. In addition to these conventional algorithms, integration of the model can be undertaken directly from the hyperplasticity relationships, rather than the incremental formulation, as shown by Ghoreishian Amiri et al. (2023).

Acknowledgements

This paper was motivated by general studies of the mechanical behaviour of soft natural soils in a research project named SAUNA (‘safety of urbanised natural slopes’). SAUNA is an industry–academic collaborative project on the safety assessment of natural slopes in Norway and is sponsored by the Norwegian Public Roads Administration (SVV), the Norwegian Water Resources and Energy Directorate (NVE), Bane NOR and the Norwegian University of Science and Technology (NTNU). This work was also supported by the Research Council of Norway through its Centres of Excellence funding scheme, grant number 262644- Porelab. All support is sincerely acknowledged. Any opinions, findings, conclusions, or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the SAUNA project collaborators or the Research Council of Norway. Special gratitude is extended to Arnstein Watn, Professor Steinar Nordal and Professor Gudmund Reidar Eiksund, for their contributions in enabling this research.

NOTATION

d: scalar valued dissipation function
f: Helmholtz free energy potential
g: dimensionless coefficient for elastic shear modulus
g_p: dimensionless coefficient for kinematic hardening modulus
H_i: weight function
K_i: weight function
k_p: dimensionless kinematic coefficient for hardening modulus
M: slope of critical state line in p–q stress space
n: rate sensitivity parameter (dimensionless)
p: mean effective pressure
p₀: isotropic pre-consolidation pressure
p_a: reference pressure (atmospheric pressure)
p_b: spherical back-stress
p_eq: equivalent pressure
q: deviatoric stress invariant
q_b: deviatoric back-stress
R: spacing ratio
r: arbitrary reference strain rate
S: state variable
T: transition function
w: flow potential
z: force potential
ε_s, ε_s^p: total and plastic deviatoric strains
ε_v, ε_v^p: total and plastic volumetric strains
η: stress ratio in p–q stress space
η_χ: stress ratio in dissipative stress space
κ: swelling index
λ: compression index
μ: creep index
τ: arbitrary reference equivalent time (normally 24 h)
χ_p, χ_q: spherical and deviatoric dissipative stresses

REFERENCES

Adachi

,

T.

&

Okano

,

M.

(

1974

).

A constitutive equation for normally consolidated clay

.

Soils Found.

14

, No.

4

,

55

–

73

.

Google Scholar

Crossref

Adachi

,

T.

,

Oka

,

F.

,

Hirata

,

T.

,

Hashimoto

,

T.

,

Nagaya

,

J.

,

Mimura

,

M.

&

Pradhan

,

T. B. S.

(

1995

).

Stress–strain behavior and yielding characteristics of eastern Osaka clay

.

Soils Found.

35

, No.

3

,

1

–

13

.

Google Scholar

Crossref

Apriadi

,

D.

,

Likitlersuang

,

S.

&

Pipatpongsa

,

T.

(

2013

).

Loading path dependence and non-linear stiffness at small strain using rate-dependent multisurface hyperplasticity model

.

Comput. Geotech.

49

,

100

–

110

.

Google Scholar

Crossref

Arulanandan

,

K.

,

Shen

,

C. K.

&

Young

,

R. B.

(

1971

).

Undrained creep behaviour of a coastal organic silty clay

.

Géotechnique

21

, No.

4

,

359

–

375

,

https://doi.org/10.1680/geot.1971.21.4.359

.

Google Scholar

Crossref

Barden

,

L.

(

1968

).

Primary and secondary consolidation of clay and peat

.

Géotechnique

18

, No.

1

,

1

–

24

,

https://doi.org/10.1680/geot.1968.18.1.1

.

Google Scholar

Crossref

Bathayian

,

S. M. H.

&

Maleki

,

M.

(

2023

).

Kinematic hardening based coupled elastoplastic–viscoplastic model for describing time-dependent behavior of soils subjected to non-monotonic loadings

.

Comput. Geotech.

161

,

105602

.

Google Scholar

Crossref

Bjerrum

,

L.

(

1967

).

Engineering geology of Norwegian normally-consolidated marine clays as related to settlements of buildings

.

Géotechnique

17

, No.

2

,

83

–

118

,

https://doi.org/10.1680/geot.1967.17.2.83

.

Google Scholar

Crossref

Bodas Freitas

,

T. M.

,

Potts

,

D. M.

&

Zdravkovic

,

L.

(

2011

).

A time dependent constitutive model for soils with isotach viscosity

.

Comput. Geotech.

38

, No.

6

,

809

–

820

.

Google Scholar

Crossref

Borja

,

R. I.

,

Tamagnini

,

C.

&

Amorosi

,

A.

(

1997

).

Coupling plasticity and energy-conserving elasticity models for clays

.

J. Geotech. Geoenviron. Engng

123

, No.

10

,

948

–

957

.

Google Scholar

Crossref

Burland

,

J. B.

(

1965

).

The yielding and dilation of clay

.

Géotechnique

15

, No.

2

,

211

–

214

,

https://doi.org/10.1680/geot.1965.15.2.211

.

Google Scholar

Crossref

Chaboche

,

J. L.

(

1989

).

Constitutive equations for cyclic plasticity and cyclic viscoplasticity

.

Int. J. Plast.

5

, No.

3

,

247

–

302

.

Google Scholar

Crossref

Chaboche

,

J. L.

(

2008

).

A review of some plasticity and viscoplasticity constitutive theories

.

Int. J. Plast.

24

, No.

10

,

1642

–

1693

.

Google Scholar

Crossref

Collins

,

I. F.

(

2002

).

Associated and non-associated aspects of the constitutive laws for coupled elastic/plastic materials

.

Int. J. Geomech.

2

, No.

2

,

259

–

267

.

Google Scholar

Crossref

Collins

,

I. F.

(

2003

).

A systematic procedure for constructing critical state models in three dimensions

.

Int. J. Solids Structs

40

, No.

17

,

4379

–

4397

.

Google Scholar

Crossref

Collins

,

I. F.

(

2005

).

The concept of stored plastic work or frozen elastic energy in soil mechanics

.

Géotechnique

55

, No.

5

,

373

–

382

,

https://doi.org/10.1680/geot.2005.55.5.373

.

Google Scholar

Crossref

Collins

,

I. F.

&

Hilder

,

T.

(

2002

).

A theoretical framework for constructing elastic/plastic constitutive models of triaxial tests

.

Int. J. Numer. Analyt. Methods Geomech.

26

, No.

13

,

1313

–

1347

.

Google Scholar

Crossref

Collins

,

I. F.

&

Houlsby

,

G. T.

(

1997

).

Application of thermomechanical principles to the modelling of geotechnical materials

.

Proc. R. Soc. A: Math. Phys. Engng Sci.

453

, No.

1964

,

1975

–

2001

.

Google Scholar

Crossref

Collins

,

I. F.

&

Kelly

,

P. A.

(

2002

).

A thermomechanical analysis of a family of soil models

.

Géotechnique

52

, No.

7

,

507

–

518

,

https://doi.org/10.1680/geot.2002.52.7.507

.

Google Scholar

Crossref

Coussy

,

O.

(

1995

).

Mechanics of porous continua

.

Chichester, UK

:

Wiley

.

Google Scholar

Dadras-Ajirlou

,

D.

,

Grimstad

,

G.

&

Ghoreishian Amiri

,

S. A.

(

2023

a).

On the isotache viscous modelling of clay behaviour using the hyperplasticity approach

.

Géotechnique

73

, No.

12

,

1112

–

1124

,

https://doi.org/10.1680/jgeot.21.00245

.

Google Scholar

Crossref

Dadras-Ajirlou

,

D.

,

Grimstad

,

G.

,

Ghoreishian Amiri

,

S. A.

&

Nordal

,

S.

(

2023

b).

A set of hyper-viscoplastic critical state models with different friction mobilisation criteria

.

Int. J. Solids Structs

273

,

112267

,

https://doi.org/10.1016/j.ijsolstr.2023.112267

.

Google Scholar

Crossref

Dafalias

,

Y. F.

(

1986

).

Bounding surface plasticity. I: Mathematical foundation and hypoplasticity

.

J. Engng Mech.

112

, No.

9

,

966

–

987

.

Google Scholar

Dafalias

,

Y. F.

&

Popov

,

E. P.

(

1975

).

A model of nonlinearly hardening materials for complex loading

.

Acta Mech.

21

, No.

3

,

173

–

192

.

Google Scholar

Crossref

Dafalias

,

Y. F.

&

Taiebat

,

M.

(

2013

).

Anatomy of rotational hardening in clay plasticity

.

Géotechnique

63

, No.

16

,

1406

–

1418

,

https://doi.org/10.1680/geot.12.P.197

.

Google Scholar

Crossref

Drucker

,

D. C.

(

1959

).

A definition of stable inelastic material

.

J. Appl. Mech.

26

, No.

1

,

101

–

106

.

Google Scholar

Crossref

Einav

,

I.

&

Collins

,

I. F.

(

2008

).

A thermomechanical framework of plasticity based on probabilistic micromechanics

.

J. Mech. Mater. Struct.

3

, No.

5

,

867

–

892

.

Google Scholar

Crossref

Einav

,

I.

&

Puzrin

,

A. M.

(

2003

).

Evaluation of continuous hyperplastic critical state (CHCS) model

.

Géotechnique

53

, No.

10

,

901

–

913

,

https://doi.org/10.1680/geot.2003.53.10.901

.

Google Scholar

Crossref

Einav

,

I.

&

Puzrin

,

A. M.

(

2004

).

Continuous hyperplastic critical state (CHCS) model: derivation

.

Int. J. Solids Structs

41

, No.

1

,

199

–

226

.

Google Scholar

Crossref

Ghoreishian Amiri

,

S. A.

,

Dadras-Ajirlou

,

D.

&

Grimstad

,

G.

(

2023

). Direct integration method for hyperplastic models. In

10th European conference on numerical methods in geotechnical engineering

(eds

L.

Zdravkovic

,

S.

Kontoe

,

D.

Taborda

and

A.

Tsiampousi

),

https://doi.org/10.53243/NUMGE2023-32

.

London, UK

:

International Society for Soil Mechanics and Geotechnical Engineering

.

Google Scholar

Grimstad

,

G.

,

Degago

,

S. A.

,

Nordal

,

S.

&

Karstunen

,

M.

(

2010

).

Modeling creep and rate effects in structured anisotropic soft clays

.

Acta Geotech.

5

, No.

1

,

69

–

81

.

Google Scholar

Crossref

Grimstad

,

G.

,

Dadras-Ajirlou

,

D.

&

Amiri

,

S. A. G.

(

2020

).

Modelling creep in clay using the framework of hyper-viscoplasticity

.

Géotechnique Lett.

10

, No.

3

,

404

–

408

,

https://doi.org/10.1680/jgele.20.00004

.

Google Scholar

Crossref

Grimstad

,

G.

,

Long

,

M.

,

Dadras-Ajirlou

,

D.

&

Amiri

,

S. A. G.

(

2021

).

Investigation of development of the earth pressure coefficient at rest in clay during creep in the framework of hyper-viscoplasticity

.

Int. J. Geomech.

21

, No.

1

,

04020235

.

Google Scholar

Crossref

Hashiguchi

,

K.

(

1989

).

Subloading surface model in unconventional plasticity

.

Int. J. Solids Structs

25

, No.

8

,

917

–

945

.

Google Scholar

Crossref

Heeres

,

O. M.

,

Suiker

,

A. S. J.

&

de Borst

,

R.

(

2002

).

A comparison between the Perzyna viscoplastic model and the consistency viscoplastic model

.

Eur. J. Mech. – A/Solids

21

, No.

1

,

1

–

12

.

Google Scholar

Crossref

Hicher

,

P. Y.

(

2016

).

Experimental study of viscoplastic mechanisms in clay under complex loading

.

Géotechnique

66

, No.

8

,

661

–

669

,

https://doi.org/10.1680/jgeot.15.P.203

.

Google Scholar

Crossref

Houlsby

,

G. T.

(

1981

).

Study of plasticity theories and their applicability to soils

.

PhD thesis

,

University of Cambridge

,

Cambridge, UK

.

Houlsby

,

G. T.

(

2000

). Critical state models and small-strain stiffness. In

Developments in theoretical geomechanics: Proceedings of the Booker memorial symposium,

Sydney, NSW, Australia

(eds

D. W.

Smith

and

J. P.

Carter

), pp.

295

–

312

.

Rotterdam, The Netherlands

:

A.A. Balkema

.

Google Scholar

Houlsby

,

G. T.

(

2020

). Hyperplasticity: from micro to macro. In

Views on microstructures in granular materials

(eds

P.

Giovine

,

P. M.

Mariano

and

G.

Mortara

), pp.

87

–

106

.

Cham, Switzerland:

Springer International Publishing

.

Google Scholar

Crossref

Houlsby

,

G. T.

&

Puzrin

,

A. M.

(

2006

).

Principles of hyperplasticity: an approach to plasticity theory based on thermodynamic principles

.

London, UK

:

Springer

.

Google Scholar

Houlsby

,

G. T.

&

Richards

,

I. A.

(

2023

).

Multi-surface and bounding surface models in hyperplasticity

.

Comput. Geotech.

156

,

105143

,

https://doi.org/10.1016/j.compgeo.2022.105143

.

Google Scholar

Crossref

Houlsby

,

G. T.

,

Amorosi

,

A.

&

Rojas

,

E.

(

2005

).

Elastic moduli of soils dependent on pressure: a hyperelastic formulation

.

Géotechnique

55

, No.

5

,

383

–

392

,

https://doi.org/10.1680/geot.2005.55.5.383

.

Google Scholar

Crossref

Islam

,

M. N.

&

Gnanendran

,

C. T.

(

2017

).

Elastic-viscoplastic model for clays: development, validation, and application

.

J. Engng Mech.

143

, No.

10

,

04017121

.

Google Scholar

Janbu

,

N.

(

1969

).

The resistance concept applied to deformations of soils

.

Proceedings of the 7th international conference on soil mechanics and foundation engineering

, pp.

191

–

196

.

Mexico City, Mexico

:

Sociedad Mexicana de Mecánica de Suelos

.

Google Scholar

Janbu

,

N.

(

1985

).

Soil models in offshore engineering

.

Géotechnique

35

, No.

3

,

241

–

281

,

https://doi.org/10.1680/geot.1985.35.3.241

.

Google Scholar

Crossref

Jiang

,

J.

,

Ling

,

H. I.

,

Kaliakin

,

V. N.

,

Zeng

,

X.

&

Hung

,

C.

(

2017

).

Evaluation of an anisotropic elastoplastic–viscoplastic bounding surface model for clays

.

Acta Geotech.

12

, No.

2

,

335

–

348

.

Google Scholar

Crossref

Jostad

,

H. P.

&

Yannie

,

J.

(

2017

).

A procedure for determining long-term creep rates of soft clays by triaxial testing

.

Eur. J. Environ. Civ. Engng

26

, No.

7

,

2600

–

2615

.

Google Scholar

Crossref

Kavvadas

,

M.

&

Amorosi

,

A.

(

2000

).

A constitutive model for structured soils

.

Géotechnique

50

, No.

3

,

263

–

273

,

https://doi.org/10.1680/geot.2000.50.3.263

.

Google Scholar

Crossref

Kimoto

,

S.

,

Khan

,

B. S.

,

Mirjalili

,

M.

&

Oka

,

F.

(

2015

).

Cyclic elastoviscoplastic constitutive model for clay considering nonlinear kinematic hardening rules and structural degradation

.

Int. J. Geomech.

15

, No.

5

,

A4014005

.

Google Scholar

Crossref

Krabbenhøft

,

K.

&

Krabbenhøft

,

J.

(

2021

).

Simplified kinematic hardening plasticity framework for constitutive modelling of soils

.

Comput. Geotech.

138

,

104146

.

Google Scholar

Crossref

Kutter

,

B. L.

&

Sathialingam

,

N.

(

1992

).

Elastic-viscoplastic modelling of the rate-dependent behaviour of clays

.

Géotechnique

42

, No.

3

,

427

–

441

,

https://doi.org/10.1680/geot.1992.42.3.427

.

Google Scholar

Crossref

Leroueil

,

S.

(

2006

).

The isotache approach. Where are we 50 years after its development by Professor Šuklje? (2006 Prof. Šuklje's memorial lecture)

. In

Proceedings of the 13th Danube European conference on geotechnical engineering

(ed. J. Logar), pp.

55

–

88

.

Ljubljana, Slovenia

:

Slovenian Geotechnical Society

.

Google Scholar

Li

,

X. S.

&

Dafalias

,

Y. F.

(

2015

).

Dissipation consistent fabric tensor definition from DEM to continuum for granular media

.

J. Mech. Phys. Solids

78

,

141

–

153

.

Google Scholar

Crossref

Li

,

J.

,

Zhang

,

D. M.

,

Hicher

,

P. Y.

&

Wei

,

C. F.

(

2022

).

Three-dimensional modelling of stress relaxation of soft clay under complex loading conditions

.

Eur. J. Environ. Civ. Engng

26

, No.

7

,

2496

–

2508

.

Google Scholar

Crossref

Mánica

,

M. A.

,

Gens

,

A.

,

Ovando-Shelley

,

E.

,

Botero

,

E.

&

Vaunat

,

J.

(

2021

).

An effective combined framework for modelling the time-dependent behaviour of soft structured clays

.

Acta Geotech.

16

, No.

2

,

535

–

550

.

Google Scholar

Crossref

Maranha

,

J. R.

,

Pereira

,

C.

&

Vieira

,

A.

(

2016

).

A viscoplastic subloading soil model for rate-dependent cyclic anisotropic structured behaviour

.

Int. J. Numer. Analyt. Methods Geomech.

40

, No.

11

,

1531

–

1555

.

Google Scholar

Crossref

Mróz

,

Z.

(

1967

).

On the description of anisotropic workhardening

.

J. Mech. Phys. Solids

15

, No.

3

,

163

–

175

.

Google Scholar

Crossref

Mróz

,

Z.

,

Norris

,

V. A.

&

Zienkiewicz

,

O. C.

(

1979

).

Application of an anisotropic hardening model in the analysis of elasto–plastic deformation of soils

.

Géotechnique

29

, No.

1

,

1

–

34

,

https://doi.org/10.1680/geot.1979.29.1.1

.

Google Scholar

Crossref

Perzyna

,

P.

(

1963

).

The constitutive equations for rate sensitive plastic materials

.

Q. Appl. Math.

20

, No.

4

,

321

–

332

.

Google Scholar

Crossref

Perzyna

,

P.

(

1966

). Fundamental problems in viscoplasticity. In

Advances in applied mechanics

(eds

G. G.

Chernyi

,

H. L.

Dryden

,

P.

Germain

,

L.

Howarth

,

W.

Olszak

,

W.

Prager

,

R. F.

Probstein

and

H.

Ziegler

), vol.

9

, pp.

243

–

377

.

New York, NY, USA

:

Academic Press

.

Google Scholar

Popova

,

E.

&

Popov

,

V. L.

(

2015

).

The research works of Coulomb and Amontons and generalized laws of friction

.

Friction

3

, No.

2

,

183

–

190

.

Google Scholar

Crossref

Puzrin

,

A. M.

&

Houlsby

,

G. T.

(

2001

).

Fundamentals of kinematic hardening hyperplasticity

.

Int. J. Solids Structs

38

, No.

21

,

3771

–

3794

.

Google Scholar

Crossref

Puzrin

,

A. M.

&

Houlsby

,

G. T.

(

2003

).

Rate-dependent hyperplasticity with internal functions

.

J. Engng Mech.

129

, No.

3

,

252

–

263

.

Google Scholar

Qiao

,

Y.

,

Ferrari

,

A.

,

Laloui

,

L.

&

Ding

,

W.

(

2016

).

Nonstationary flow surface theory for modeling the viscoplastic behaviors of soils

.

Comput. Geotech.

76

,

105

–

119

.

Google Scholar

Crossref

Radjai

,

F.

,

Wolf

,

D. E.

,

Jean

,

M.

&

Moreau

,

J. J.

(

1998

).

Bimodal character of stress transmission in granular packings

.

Phys. Rev. Lett.

80

, No.

1

,

61

–

64

.

Google Scholar

Crossref

Rollo

,

F.

&

Amorosi

,

A.

(

2020

).

SANICLAY-T: simple thermodynamic-based anisotropic plasticity model for clays

.

Comput. Geotech.

127

,

103770

.

Google Scholar

Crossref

Roscoe

,

K. H.

&

Burland

,

J. B.

(

1968

).

On the generalized stress–strain behaviour of wet clay

. In

Engineering plasticity

(eds

J.

Heyman

and

F. A.

Leckie

), pp.

535

–

609

.

Cambridge, UK

:

Cambridge University Press

.

Google Scholar

Roscoe

,

K. H.

,

Schofield

,

A. N.

&

Wroth

,

C. P.

(

1958

).

On the yielding of soils

.

Géotechnique

8

, No.

1

,

22

–

53

,

https://doi.org/10.1680/geot.1958.8.1.22

.

Google Scholar

Crossref

Rouainia

,

M.

&

Muir Wood

,

D.

(

2000

).

A kinematic hardening constitutive model for natural clays with loss of structure

.

Géotechnique

50

, No.

2

,

153

–

164

,

https://doi.org/10.1680/geot.2000.50.2.153

.

Google Scholar

Crossref

Schofield

,

A. N.

&

Wroth

,

C. P.

(

1968

).

Critical state soil mechanics

.

New York, NY, USA

:

McGraw-Hill

.

Google Scholar

Shahbodagh

,

B.

,

Mac

,

T. N.

,

Esgandani

,

G. A.

&

Khalili

,

N.

(

2020

).

A bounding surface viscoplasticity model for time-dependent behavior of soils including primary and tertiary creep

.

Int. J. Geomech.

20

, No.

9

,

04020143

.

Google Scholar

Crossref

Sheahan

,

T. C.

,

Ladd

,

C. C.

&

Germaine

,

J. T.

(

1996

).

Rate-dependent undrained shear behavior of saturated clay

.

J. Geotech. Engng

122

, No.

2

,

99

–

108

.

Google Scholar

Crossref

Shi

,

Z.

,

Muir Wood

,

D.

,

Huang

,

M.

&

Hambleton

,

J. P.

(

2023

).

Tay creep: a multi-mechanism model for rate-dependent deformation of soils

.

Géotechnique

73

, No.

4

,

310

–

322

,

https://doi.org/10.1680/jgeot.21.00084

.

Google Scholar

Crossref

Simo

,

J. C.

&

Hughes

,

T. J.

(

1998

).

Computational inelasticity

.

New York, NY, USA

:

Springer Science & Business Media

.

Google Scholar

Sloan

,

S. W.

(

1987

).

Substepping schemes for the numerical integration of elastoplastic stress–strain relations

.

Int. J. Numer. Methods Engng

24

, No.

5

,

893

–

911

.

Google Scholar

Crossref

Stallebrass

,

S. E.

&

Taylor

,

R. N.

(

1997

).

The development and evaluation of a constitutive model for the prediction of ground movements in overconsolidated clay

.

Géotechnique

47

, No.

2

,

235

–

253

,

https://doi.org/10.1680/geot.1997.47.2.235

.

Google Scholar

Crossref

Šuklje

,

L.

(

1957

).

The analysis of the consolidation process by the isotache method

. In

Proceedings of the 4th international conference on soil mechanics and foundation engineering

, pp.

200

–

206

.

London, UK

:

Butterworths

.

Google Scholar

Tafili

,

M.

&

Triantafyllidis

,

T.

(

2020

).

AVISA: anisotropic visco-ISA model and its performance at cyclic loading

.

Acta Geotech.

15

, No.

9

,

2395

–

2413

.

Google Scholar

Crossref

Tafili

,

M.

,

Wichtmann

,

T.

&

Triantafyllidis

,

T.

(

2021

).

Experimental investigation and constitutive modeling of the behaviour of highly plastic lower Rhine clay under monotonic and cyclic loading

.

Can. Geotech. J.

58

, No.

9

,

1396

–

1410

.

Google Scholar

Crossref

Vaid

,

Y. P.

&

Campanella

,

R. G.

(

1977

).

Time-dependent behavior of undisturbed clay

.

J. Geotech. Engng Div.

103

, No.

7

,

693

–

709

.

Google Scholar

Crossref

Vermeer

,

P.

&

Neher

,

H.

(

1999

). A soft soil model that accounts for creep. In

Beyond 2000 in computational geotechnics

(ed.

R. B. J.

Brinkgreve

), pp.

249

–

261

.

Abingdon, UK

:

Routledge

.

Google Scholar

Wroth

,

C. P.

&

Houlsby

,

G. T.

(

1985

).

Soil mechanics: property characterization and analysis procedures

. In

Proceedings of the 11th international conference on soil mechanics and foundation engineering

,

San Francisco, CA, USA

, vol.

1

, pp.

1

–

55

.

Rotterdam, Netherlands

:

A.A. Balkema

.

Google Scholar

Yang

,

C.

,

Carter

,

J. P.

,

Sheng

,

D.

&

Sloan

,

S. W.

(

2016

).

An isotach elastoplastic constitutive model for natural soft clays

.

Comput. Geotech.

77

,

134

–

155

.

Google Scholar

Crossref

Yao

,

Y. P.

,

Kong

,

L. M.

,

Zhou

,

A. N.

&

Yin

,

J. H.

(

2015

).

Time-dependent unified hardening model: three-dimensional elastoviscoplastic constitutive model for clays

.

J. Engng Mech.

141

, No.

6

,

04014162

.

Google Scholar

Yin

,

J. H.

&

Graham

,

J.

(

1999

).

Elastic viscoplastic modelling of the time-dependent stress–strain behaviour of soils

.

Can. Geotech. J.

36

, No.

4

,

736

–

745

.

Google Scholar

Crossref

Yin

,

J. H.

&

Zhu

,

J. G.

(

1999

).

Measured and predicted time-dependent stress–strain behaviour of Hong Kong marine deposits

.

Can. Geotech. J.

36

, No.

4

,

760

–

766

.

Google Scholar

Crossref

Yin

,

J. H.

,

Zhu

,

J. G.

&

Graham

,

J.

(

2002

).

A new elastic viscoplastic model for time-dependent behaviour of normally and overconsolidated clays: theory and verification

.

Can. Geotech. J.

39

, No.

1

,

157

–

173

.

Google Scholar

Crossref

Yu

,

H. S.

(

1998

).

CASM: a unified state parameter model for clay and sand

.

Int. J. Numer. Analyt. Methods Geomech.

22

, No.

8

,

621

–

653

.

Google Scholar

Crossref

Yuan

,

Y.

&

Whittle

,

A. J.

(

2021

).

Formulation of a new elastoviscoplastic model for time-dependent behavior of clay

.

Int. J. Numer. Analyt. Methods Geomech.

45

, No.

6

,

843

–

864

.

Google Scholar

Crossref

Zhu

,

J. G.

(

2000

).

Experimental study and elastic visco-plastic modelling of the time-dependent stress–strain behaviour of Hong Kong marine deposits

.

PhD thesis

,

Hong Kong Polytechnic University

,

Hong Kong SAR, China

.

Ziegler

,

H.

(

1972

).

Systems with internal parameters obeying the orthogonality condition

.

Z. Angew. Math. Phys. ZAMP

23

, No.

4

,

553

–

566

.

Google Scholar

Crossref

Ziegler

,

H.

(

1977

).

An Introduction to thermodynamics

.

Amsterdam, The Netherlands

:

North-Holland

.

Google Scholar

Ziegler

,

H.

(

1981

).

Discussion of some objections to thermomechanical orthogonality

.

Ing.-Arch.

50

, No.

3

,

149

–

164

.

Google Scholar

Crossref

Zytynski

,

M.

,

Randolph

,

M. F.

,

Nova

,

R.

&

Wroth

,

C. P.

(

1978

).

On modelling the unloading–reloading behaviour of soils

.

Int. J. Numer. Analyt. Methods Geomech.

2

, No.

1

,

87

–

93

.

Google Scholar

Crossref

Discussion on this paper closes 1 May 2026; for further details see p. ii.

2025

Emerald Publishing Limited

Published by Emerald Publishing Limited. This article is published under the Creative Commons Attribution (CC BY 4.0) licence. Anyone may reproduce, distribute, translate and create derivative works of this article (for both commercial and non-commercial purposes), subject to full attribution to the original publication and authors. The full terms of this licence may be seen at the terms of the CC BY 4.0 licenceLink to the terms of the CC BY 4.0 licence.

Modified Cam Clay bounding surface hyper-viscoplastic model

INTRODUCTION

FORCE, FLOW AND FREE ENERGY POTENTIALS

INCREMENTAL FORMULATION

GENERALISED FORCE AND FLOW POTENTIALS

MODEL PARAMETERS

EVALUATION OF MODEL

CONCLUSION

Acknowledgements

NOTATION

REFERENCES

Email Alerts

Cited By

Modified Cam Clay bounding surface hyper-viscoplastic model Open Access

INTRODUCTION

FORCE, FLOW AND FREE ENERGY POTENTIALS

INCREMENTAL FORMULATION

GENERALISED FORCE AND FLOW POTENTIALS

MODEL PARAMETERS

EVALUATION OF MODEL

CONCLUSION

Acknowledgements

NOTATION

REFERENCES

Email Alerts

Suggested Reading

Related Chapters

Recommended for you

Recommended for you

Cited By

Modified Cam Clay bounding surface hyper-viscoplastic model