<figure><img src="image_1.png" /></figure>

HAL
open science

# A simple critical state based double-yield-surface model for clay behavior under complex loading

Zhenyu Yin, Qiang Xu, Pierre-Yves Hicher

► To cite this version:

Zhenyu Yin, Qiang Xu, Pierre-Yves Hicher. A simple critical state based double-yield-surface model for clay behavior under complex loading. Acta Geotechnica, 2013, 8 (5), pp.509-523. 10.1007/s11440-013-0206-y. hal-01004941

HAL Id: hal-01004941

https://hal.science/hal-01004941

Submitted on 8 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# A simple critical-state-based double-yield-surface model for clay behavior under complex loading

Zhen-Yu Yin · Qiang Xu · Pierre-Yves Hicher

Z.-Y. Yin
Department of Civil Engineering,
Shanghai Jiao Tong University,
Shanghai 200240, People's Republic of China

Z.-Y. Yin · Q. Xu
State Key Laboratory of Geohazard Prevention
and Geoenvironment Protection, Chengdu University
of Technology, Chengdu 610059, People's Republic of China

Z.-Y. Yin · P.-Y. Hicher
Research Institute in Civil and Mechanical Engineering,
UMR CNRS 6183, Ecole Centrale de Nantes,
44300 Nantes, France

**Abstract** The critical state concept has been widely used in soil mechanics. The purpose of this study is to apply this concept in the framework of multi-mechanism elastoplasticity. The developed model has two yield surfaces: one for shear sliding and one for compression. In this model, the location of the critical state line is explicitly considered and related to the actual material density to control the peak strength and the phase transformation characteristics. The stress reversal technique is incorporated into the model for describing clay behavior under complex loading including changes of stress direction. The determination of the model parameters is discussed; it requires only one drained or undrained triaxial test up to failure with an initial isotropic consolidation stage. The model is used to simulate drained and undrained tests under monotonic loading with different over-consolidation ratios on various remolded and natural clays, including true triaxial tests with different Lode's angles. Drained and undrained tests under cyclic loadings are also simulated by using the set of parameters determined from monotonic tests. The comparison between experimental results and numerical simulations demonstrate a good predictive ability of this new simple model.

**Keywords** Clay · Constitutive modeling · Critical state · Cyclic loading · Elastoplasticity

## 1 Introduction

The mechanical behavior of clay has been widely investigated by numerous researchers during the last decades [3, 8, 11, 18, 20, 23, 33, 42]. A clayey material is contractive (decreasing void ratio in drained condition in Fig. 1a and decreasing mean effective stress in undrained condition in Fig. 1b) in a normally or lightly over-consolidated state and becomes dilative (increasing void ratio in drained condition in Fig. 1a and increasing mean effective stress in undrained condition in Fig. 1b) when the over-consolidation ratio increases, with the stress state and the void ratio converging to the critical state (CSL in $p'-q$ plane and CSL in $e-\log p'$ ← plane, shown in Fig. 1a, b). The location of the CSL in the $e-\log p'$ ← plane is a key factor for determining the magnitude of contraction or dilation (see difference between B and $B'$, C and $C'$ in Fig. 1c) from an initial given void ratio.

Based on experimental observations, conventional elastoplastic models for clay have been developed, and generally can be classified into two categories:

(1) Critical-state-based elastoplastic models [4, 9, 20, 21, 23, 25–27, 30, 33, 34, 36, 41]. In these models, the slope of CSL in the $e-\log p'$ ← plane — $\lambda$ — $\psi$ is explicitly used as a parameter, and the location of this CSL is implicitly governed by the intersection point of the yield surface with the CSL in $p'-q$ plane.
<figure><img src="image_2.png" /><figcaption>Fig. 1 Schematic plot for the mechanical behavior of clay a drained behavior; b undrained behavior; c influence of critical state line</figcaption></figure>

(2) Multi-mechanism models [13, 14, 16, 19, 24]. The model by Hujeux [16] used the location of CSL as presented in the first category, and twelve parameters are needed. In other models, the CSL in the e-logp' ← plane has not been used. As a result, the contraction and dilation up to critical states during compression and shearing are not guaranteed. Thus, these models are difficult to simulate clay behavior along various stress paths and loading conditions.

Different from them, other models have also been
developed for clay. For instance, Mašin [22] and Huang
et al. [15] proposed clay models based on the hypoplas-
ticity framework; Calladine [5] and Schweiger et al. [28]
proposed clay models based on the concept of multi-lam-
inate modeling. In these models, the CSL in the e-logp' ←
plane has not been used, thus the contraction and dilation
up to critical states cannot be guaranteed.

In sand modeling, the critical state concept has been
explicitly used [2, 10, 31]. This explicit method has also
been adopted by the authors to account for the interlocking
effect and neighboring effect of sand grains and clay
aggregates at inter-particle contact level in different ver-
sions of the micromechanics-based model [6, 7, 37–40].
Good predictive performance has been achieved for both
sand and clay, and note that the models are physically
attractive on describing the behavior of anisotropy. How-
ever, the micromechanics-based models are usually com-
putational demand. Therefore, this paper is an attempt to
explicitly use the CSL in the double-yield-surface-model-
ing method in order to propose a model simple to use and
at the same time well adapted to the modeling of clay
behavior.

In the following sections, the development of the dou-
ble-yield-surface model using the critical state concept for
clay is firstly presented. The implementation of the stress

reversal technique into the model is then presented in order
to be able to simulate tests with shear stress reversal.
Finally, the model is validated by comparing experimental
results and numerical simulations of drained and undrained
tests on isotropically and anisotropically consolidated clay
samples with different OCRs (over-consolidation ratio)
under both monotonic and cyclic loadings.

## 2 Constitutive model

According to elastoplasticity theory, the total strain rate is
additively composed of the elastic strain rates and plastic
strain rates

$$
\dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{e} + \dot{\epsilon}_{ij}^{p} \qquad (1)
$$

where $\dot{\epsilon}_{ij}$ denotes the $(i,j)$ component of the total strain rate tensor, and the superscripts $e$ and $p$ stand, respectively, for the elastic and plastic components.

2.1 Elastic behavior

The elastic behavior is assumed to be isotropic, as follows:

$$
\dot{\epsilon}_{ij}^{e} = \frac{1+v}{E} \sigma'_{ij} - \frac{v}{E} \sigma'_{kk} \delta_{ij} \quad (2)
$$

where $v$ and $E$ are Poisson's ratio and Young's modulus; $\sigma'_{ij}$ is the effective stress tensor; $\delta_{ij}$ is Kronecker's delta.

$E$ can also be replaced by the elastic bulk modulus $K$ by
$E = 3K(1 - 2v)$. For clay, the elastic bulk modulus $K$ can
be determined by the relation $K = p'(1 + e_0)/\kappa$, with the
mean effective stress $p' = \sigma'_{kk}/3$ (see Roscoe and Burland
[26]), the slope of the swelling line $\kappa$ and the initial void
ratio $e_0$. Assuming a common value of $v = 0.25$ for clay,
only the parameter $\kappa$ is then needed for the elastic
behavior.
## 2.2 Plastic behavior

Based on conventional elastoplasticity theory, the plastic strain can be obtained, as follows:

$$ \dot{\epsilon}_{ij}^p = d\lambda \frac{\partial g}{\partial \sigma'_{ij}} \qquad (3) $$

where $d\lambda$ is plastic multiplier depending on the stress rate and the plastic hardening law; $g$ is plastic potential.

The proposed model uses two yield surfaces: one for shear sliding ($f_S$) and one for compression ($f_C$) (see Fig. 2). Equation (3) can be expressed as follows:

$$ \dot{\epsilon}_{ij}^p = \dot{\epsilon}_{ij}^{pS} + \dot{\epsilon}_{ij}^{pC} = d\lambda^S \frac{\partial g^S}{\partial \sigma'_{ij}} + d\lambda^C \frac{\partial g^C}{\partial \sigma'_{ij}} \qquad (4) $$

where the superscripts S and C represent the shear sliding and compression components, respectively. For one component with $f < 0$, the corresponding $d\lambda$ is taken equal to zero.

### 2.2.1 CSL-related density state

One of the important elements to be considered in soil modeling is the critical state concept. At critical state, the clay material remains at a constant volume while it is subjected to a continuous distortion. The void ratio corresponding to this state is $e_c$. The critical void ratio $e_c$ is a function of the mean effective stress $p'$. The relationship has traditionally been written as follows:

$$ e_c = e_{cr0} - \lambda \ln \left( \frac{p'}{p_{cr0}} \right) \qquad (5) $$

The critical state line in the e-logp' plane is explicitly located by three parameters: $e_{cr0}$ and $p_{cr0}$ determine a reference critical state point in the e-logp' plane and $\lambda$ determines the slope of critical state line. For convenience,

<figure><img src="image_3.png" /><figcaption>Fig. 2 Principle of the critical state double-yield-surface model for clay</figcaption></figure>

the value of $p_{cr0}$ is taken equal to 100 kPa. The critical state line can then be defined by the two parameters $e_{cr0}$ and $\lambda$.

Using the critical state concept, the density state of the soil is defined as the ratio $e_c/e$, where $e$ is the actual void ratio and $e_c$ is the critical void ratio obtained by Eq. (5) for the value of $p'$ at the given stress state (see Fig. 2 for the state A). This implies: $e_c/e < 1$ for normally consolidated and slightly over-consolidated clay, $e_c/e > 1$ for heavily over-consolidated clay, and $e_c/e = 1$ when the initial state of clay lies on the CSL in the e-logp' plane.

### 2.2.2 Shear-sliding criterion

In several models for sand and silt [10, 17, 32], the shape of the yield surface for the plastic shear component is linear in the $p'-q$ plane. Here, we adopt the same shear criterion for clay, written as follows:

$$ f_S = \sqrt{\frac{3}{2} r_{ij} r_{ij} - H} \qquad (6) $$

where $r_{ij} = s_{ij}/p'$ with $s_{ij} = \sigma'_{ij} - p'\delta_{ij}$; $H$ is the hardening parameter defined by a hyperbolic function in the $H - \epsilon_d^p$ plane given by

$$ H = \frac{M_p \epsilon_d^p}{1/G_p + \epsilon_d^p} \qquad (7) $$

where $G_p$ controls the initial slope of the hyperbolic curve $\eta - \epsilon_d^p(\eta = q/p')$. Equations (6) and (7) guarantee that the stress ratio $\eta$ reaches the peak stress ratio $M_p$.

According to Biarez and Hicher [3], the peak friction angle $\phi_p$ (related to $M_{pc} = 6 \sin \phi_p / (3 - \sin \phi_p)$) depends on the internal friction angle $\phi_\mu$ (related to the critical state value $M_c = 6 \sin \phi_\mu / (3 - \sin \phi_\mu)$) and the density state of soil ($e_c/e$).

$$ \tan \phi_p = \frac{e_c}{e} \tan \phi_\mu \qquad (8) $$

Equation (8) shows that, in a loose structure, the peak friction angle $\phi_p$ is smaller than $\phi_\mu$. On the other hand, a dense structure provides higher degree of interlocking. Thus, the peak friction angle $\phi_p$ is greater than $\phi_\mu$. When the loading stress reaches the peak friction angle $\phi_p$, the dense structure dilates and the degree of interlocking relaxes. As a consequence, the peak friction angle is reduced, which results in a strain-softening phenomenon.

Adopting the idea of the OCR effect on the shear modulus proposed by Hardin [12], $G_p$ is considered as a function of OCR as follows:

$$ G_p = G_{p0} \text{OCR}^2 \qquad (9) $$

where $G_{p0}$ is an input parameter instead of $G_p$; OCR can be computed based on the size of the compression yield surface as shown in the next section.
In order to take into account dilation or contraction during shear sliding, a non-associated flow rule is introduced. We propose an explicit derivation of the potential surface, given by

$$ \frac{\partial g}{\partial \sigma'_{ij}} = \frac{\partial g}{\partial p'} \frac{\partial p'}{\partial \sigma'_{ij}} + \frac{\partial g}{\partial s_{ij}} \frac{\partial s_{ij}}{\partial \sigma'_{ij}} \quad \text{with} \\ \frac{\partial g}{\partial p'} = D \left( M_{\text{pt}} - \sqrt{\frac{3}{2} r_{ij} r_{ij}} \right); \quad \frac{\partial g}{\partial s_{ij}} = \sqrt{\frac{3}{2} \frac{r_{ij}}{\sqrt{r_{ij} r_{ij}}}} \quad (10) $$

which implies the following stress-dilatancy relationship

$$ \frac{d e_v^p}{d e_d^p} = D(M_{\text{pt}} - \eta) \quad (11) $$

where the contractancy or dilatancy depends on the sign of “$M_{\text{pt}} - \eta$”; and $D$ is a material constant controlling the magnitude and evolution of the contractancy or dilatancy. $M_{\text{pt}}$ is the slope of the phase transformation line for clay which can be derived from the friction angle $\phi_\mu$ by assuming a formula similar to the peak friction angle

$$ \tan \phi_{\text{pt}} = \left(\frac{e_c}{e}\right)^{-1} \tan \phi_{\mu} \quad (12) $$

Equation (10) combined with Eq. (12) implies that in a loose structure with $e > e_c$, the phase transformation angle $\phi_{\text{pt}}$ (related to $M_{\text{ptc}} = 6 \sin \phi_{\text{pt}} / (3 - \sin \phi_{\text{pt}})$) is bigger than $\phi_\mu$; in a dense structure with $e < e_c$, the phase transformation angle $\phi_{\text{pt}}$ is smaller than $\phi_\mu$ which allows the dense structure to be first contractive and then dilative during deviatoric loading. For both loose and dense structures, when the stress state reaches the critical state line, the void ratio $e$ becomes equal to the critical void ratio $e_c$, then zero dilation or contraction takes place. Thus, the constitutive equations guarantee that stresses and void ratio reach simultaneously the critical state in the $p'-q-e$ space.

In order to interpolate $M_p$ and $M_{pt}$ between their values for compression and those for extension by means of the Lode's angle $\theta$ (see Sheng et al. [29]), the following expressions of $M_p$ and $M_{pt}$ are proposed as follows:

$$ M_p = M_{pc} \left[ \frac{2c_1^4}{1 + c_1^4 + (1 - c_1^4) \sin 3\theta} \right]^{\frac{1}{4}}; \\ M_{pt} = M_{ptc} \left[ \frac{2c_2^4}{1 + c_2^4 + (1 - c_2^4) \sin 3\theta} \right]^{\frac{1}{4}} \quad (13) $$

where $c_1 = (3 - \sin \phi_p)/(3 + \sin \phi_p)$ and $c_2 = (3 - \sin \phi_{pt})/(3 + \sin \phi_{pt})$, assuming the same peak friction angle and phase transformation angle for different Lode's angle; the Lode's angle is expressed as $\frac{-\pi}{6} \le \theta = \frac{1}{3} \sin^{-1} \left( \frac{-3\sqrt{3}J_3}{2J_2^{3/2}} \right) \le \frac{\pi}{6}$.

$J_2$ and $J_3$ are the second and the third invariants of the deviatoric stress tensor, given by $J_2 = s_{ij}s_{ij}/2, J_3 = s_{ij}s_{jk}s_{ki}/3$.

### 2.2.3 Compression criterion

In order to describe the compressible behavior of clay, a second yield surface is added. The second yield function is assumed to be as follows:

$$ f_c = p'^2 + \frac{3}{2} \frac{s_{ij}s_{ij}}{R^2} - p_c^2 \quad (14) $$

where the constant $R$ controls the shape of the yield surface and $p_c$ is the hardening parameter controlling the size of the yield surface. The yield surface expands with the plastic volumetric strain. The hardening rule of the Cam-Clay model is adopted:

$$ d p_c = p_c \frac{1 + e_0}{\lambda - \kappa} d e_v^p \quad (15) $$

An associated flow rule is adopted for the compression behavior.

Parameter $R$ can be derived as follows:

One-dimensional straining corresponds to

$$ \left( \frac{d\varepsilon_d}{d\varepsilon_v} \right)_{K0} = \frac{2}{3} \quad (16) $$

Assuming that elastic strains are much smaller than plastic strains, Eq. (13) can be approximated by

$$ \left( \frac{d e_d^p}{d e_v^p} \right)_{K0} = \frac{2}{3} \quad (17) $$

Combining Eq. (14) with the flow rule derived from Eq. (11), for the condition of one-dimensional consolidation, $R$ can be expressed by

$$ R = \sqrt{\frac{3}{2}} \eta_{K0} \quad (18) $$

where $\eta_{K0} = 3M_c/(6 - M_c)$ for normally consolidated clay adopting Jacky's formula ($K_0 = 1 - \sin \phi_\mu$). Therefore, $R$ is not needed for input.

## 2.3 Incorporation of shear stress reversal technique

The above constitutive equations were developed based on experimental clay behavior, considering isotropically consolidated samples under monotonic loading. When the direction of the shear stress changes during the loading, the stress reversal effect on the stress-strain relationship needs to be incorporated.

Let us denote the stress state and plastic strain state at the moment of shear reversal to be the residual stresses
(marked as $\sigma'_{ij}^R$) and the residual plastic strains (marked as $\epsilon'_{ij}^R$). These two mechanical state variables have a significant influence on the subsequent shear behavior (see Fig. 3). As a consequence, the yield function (see Eq. (6)), the hardening rule (see Eq. (7)) and the flow rule (see Eq. (10)) have been revised as follows:

$$f_S = \sqrt{\frac{3}{2}(r_{ij} - r_{ij}^R)(r_{ij} - r_{ij}^R)} - H \quad (19)$$

$$H = \frac{M_p^* \epsilon_d^{p*}}{G_p + \epsilon_d^{p*}} \quad \text{with} \quad \epsilon_d^{p*} = \sqrt{\frac{2}{3}(e_{ij}^p - e_{ij}^{pR})(e_{ij}^p - e_{ij}^{pR})} \quad (20)$$

$$\frac{\partial g}{\partial p'} = D \quad M_{pt}^* - \sqrt{\frac{3}{2}(r_{ij} - r_{ij}^R)(r_{ij} - r_{ij}^R)};$$

$$\frac{\partial g}{\partial s_{ij}} = \sqrt{\frac{3}{2}n_{ij}} \quad \text{with} \quad n_{ij} = \frac{(r_{ij} - r_{ij}^R)}{\sqrt{(r_{ij} - r_{ij}^R)(r_{ij} - r_{ij}^R)}} \quad (21)$$

$$M_p^* = \sqrt{\frac{3}{2}(M_p n_{ij} - r_{ij}^R)(M_p n_{ij} - r_{ij}^R)}; \quad (22)$$

$$M_{pt}^* = \sqrt{\frac{3}{2}(M_{pt} n_{ij} - r_{ij}^R)(M_{pt} n_{ij} - r_{ij}^R)}$$

<figure><img src="image_4.png" /><figcaption>Fig. 3 Principle of stress reversal in loading and unloading</figcaption></figure>

where $r_{ij}^R = \sigma'_{ij}^R/p'^R - \delta_{ij}$ (with $p'^R = \sigma'_{ij}^R/3$) is the mobilized stress ratio at the moment of stress reversal; $e_{ij}^{pR} = \epsilon_{ij}^{pR} - \epsilon_{kk}^{pR} \delta_{ij}/3$ is the plastic deviatoric strain tensor at the moment of stress reversal. The estimation of Lode's angle to calculate $M_p$ and $M_{pt}$ is modified as $\theta = \frac{1}{3} \sin^{-1} \left( \frac{-3\sqrt{3}J_3'}{2J_2^{3/2}} \right)$ with $J_2' = r'_{ij} r'_{ij}/2, J_3' = r'_{ij} r'_{jk} r'_{ki}/3$ where $r'_{ij} = r_{ij} - r_{ij}^R$.

Equation (20) indicates that the same form of hardening rule can be used for both loading and unloading, but that some scaling process on the values of $M_p$ is required. Eq. (21) implies that the amount of dilation is different upon shear reversal. This concept is similar to that proposed by Balendran and Nemat-Nasser [1] and Gajo and Muir Wood [10]. This stress reversal technique is similar to that used in Masing's rule, in bounding surface plasticity (Taiebat and Dafalias [31]) and in subloading surface plasticity (Yamakawa et al. [35]).

## 2.4 Determination of model parameters

The model contains six material parameters divided into four groups (see Table 1). Based on their physical meanings, these parameters can be determined as follows:

(1) the slope of the swelling line $\kappa$ can be measured from the unloading curve of an isotropic compression test;

(2) the internal friction angle $\phi_{\mu}$ can be obtained from $M_c$ measured from drained or undrained triaxial tests up to failure ($\phi_{\mu} = \arcsin[3M_c/(6+M_c)]$);

(3) the slope of the CSL $\lambda$ is assumed to be the same as the slope of the isotropic compression line and can, therefore, be measured from an isotropic compression test;

(4) the position of the critical state in the e-logp' plane can be obtained from drained or undrained triaxial tests up to failure; then, the reference critical void ratio $e_{cr0}$ corresponding to $p' = 100$ kPa can be measured in the e-logp' plane;

(5) the plastic stiffness $G_{p0}$ can be obtained by fitting from the deviatoric stress-strain curve at small strain level (see Fig. 4a for the effect of $G_{p0}$ on the

<table><caption>Table 1 Values of model parameters for simulated clays</caption><thead><tr><th>Group</th><th>Parameter</th><th>Definition</th><th>Black</th><th>White</th><th>Fujinomori</th><th>Mixture</th><th>BBC</th><th>LCT</th></tr></thead><tbody><tr><td>Elasticity</td><td>κ</td><td>Slope of swelling line</td><td>0.079</td><td>0.034</td><td>0.02</td><td>0.034</td><td>0.023</td><td>0.008</td></tr><tr><td rowspan="3">Critical state</td><td>λ</td><td>Slope of compression line</td><td>0.244</td><td>0.089</td><td>0.093</td><td>0.173</td><td>0.171</td><td>0.066</td></tr><tr><td>e<sub>cr0</sub></td><td>Reference critical void ratio</td><td>1.33</td><td>0.653</td><td>0.77</td><td>0.977</td><td>0.986</td><td>0.435</td></tr><tr><td>φ<sub>c</sub></td><td>Critical friction angle</td><td>21.4°</td><td>23°</td><td>34°</td><td>19.2°</td><td>35°</td><td>30°</td></tr><tr><td>Plastic stiffness</td><td>G<sub>p0</sub></td><td>Plastic hardening modulus</td><td>50</td><td>70</td><td>115</td><td>40</td><td>155</td><td>170</td></tr><tr><td>Dilatancy</td><td>D</td><td>Dilatancy constant</td><td>1</td><td>1</td><td>1.2</td><td>0.3</td><td>1</td><td>0.8</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Parametric study of Gp0 effect on the deviatoric stress versus axial strain for undrained triaxial test on NC clay; b D effect on the effective stress path for undrained triaxial test on NC clay; c D effect on the volumetric strain versus axial strain for drained triaxial test on NC clay</figcaption></figure>

deviatoric stress–strain curve, keeping other param-
eters (White clay in Table 1) constant);

(6) the dilatancy constant *D* can be finally obtained by
curve fitting from the effective stress path of an
undrained test (see Fig. 4b for the effect of *D* on the
effective stress path, keeping other parameters (White
clay in Table 1) constant) or, alternatively, from the
evolution of the volumetric strain during a drained
test (see Fig. 4c for the effect of *D* on the evolution of
the volumetric strain).

Additionally, the model involves two state variables (the void ratio *e* and the size of the compression yield surface *p*<sub>c</sub>) which require their initial values to be determined. *e*<sub>0</sub> can be measured from the tested sample. The initial size *p*<sub>c0</sub> is determined by the consolidation history of the tested sample (i.e., the consolidated stress for reconstituted clay and the preconsolidation stress for natural clay).

In summary, all values of the model parameters and
state variables can be easily determined based on one
drained or undrained triaxial test up to failure with an
isotropic consolidation stage.

**3 Experimental verification**

The experimental verification is presented herein with
reference to drained and undrained tests results under
monotonic loading on various remolded and natural clays
(Black clay by Zervoyanis [42]; White clay by Biarez and
Hicher [3]; Fujinomori clay by Nakai and Hinokio [23];
Boston blue clay by Ladd and Varallyay [18]; Lower
Cromer Till by Gens [11]; and clay mixture by Li and
Meissner [20]). These tests were performed on clays of
different mineral contents and Atterberg limits. Figure 5
shows the classification of these clays using Casagrande’s
plasticity chart. According to this chart, the selected
experimental results consist of both low and high plastic
inorganic clays, as indicated in Fig. 5.

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Classification of clays by liquid limit $w_L$ and plasticity index $I_p$</figcaption></figure>

Drained tests on Fujinomori clay and undrained tests on
clay mixture under cyclic loadings were also simulated for
the validation.

3.1 Different over-consolidation ratios

3.1.1 Black clay

Drained triaxial tests on black kaolinite clay samples were
performed by Zervoyanis [42]. Four tests began with an
isotropic consolidation up to 800 kPa, then three of them
were unloaded to 400, 200, and 100 kPa, respectively
(OCR = 1, 2, 4, and 8), then followed by axial loading up to
failure under drained condition. The tested Black clay is a
remolded clay prepared from a slurry obtained by mixing
clay powder and water at a water content equal to two times
the liquid limit. The parameters presented in Table 1 were
calibrated from one triaxial compression test on the normally
consolidated specimen with an isotropic consolidation stage.

In Fig. 6a, the stress–strain curves for OCR = 1 and 2
show a continuous strain hardening due to a continuous
contractancy corresponding to the paths in Fig. 6c above
the critical state line, and the stress–strain curves for
OCR = 4 and 8 show strain hardening followed by strain
softening corresponding to the paths in Fig. 6c below the
critical state line, as respected by using Eq. (8). The pre-
dicted void ratio changes in Fig. 6b, c show a contractive
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Comparison between experimental results and model predictions for drained triaxial tests on *Black clay* with different OCRs: a deviatoric stress versus axial strain; b void ratio versus axial strain; c void ratio versus mean effective stress</figcaption></figure>

behavior for OCR = 1 and 2, and a dilutive behavior for OCR = 4 and 8, as respected by using Eq. (12). The paths in the *e*-log*p'* plane in Fig. 6c show that the void ratio approaches the critical state line when the stress state approaches the CSL in the *p'*-*q* plane. An overall good agreement is observed between experimental and numerical results for different OCRs.

### 3.1.2 White clay

Two drained triaxial tests on normally consolidated White clay and three undrained triaxial tests on normally and

over-consolidated White clay were reported by Biarez and Hicher [3]. In the case of drained tests, specimens were isotropically consolidated up to 400 and 800 kPa, respectively, and then loaded to failure keeping the confining stresses. In the case of undrained tests, three samples were isotropically consolidated up to 800 kPa, and two of them were unloaded to 400 and 67 kPa (OCR = 1, 2, 12). The parameters presented in Table 1 were calibrated from one drained compression test at the confining stress $\sigma'_{3c}$ = 800 kPa with an isotropic consolidation stage.

Using the calibrated parameters, the predicted test results are plotted in Fig. 7. For drained tests, the stress-

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Comparison between experimental results and model predictions for triaxial tests on *White clay* a deviatoric stress versus axial strain for drained tests with OCR = 1; b void ratio versus axial strain for drained tests with OCR = 1; c void ratio versus mean effective stress for drained tests; d deviatoric stress versus axial strain for undrained tests with different OCRs; e effective stress path for undrained tests with different OCRs; f void ratio versus mean effective stress for undrained tests</figcaption></figure>
strain curves in Fig. 7a and the void ratio change in Fig. 7b, c shows good agreement between experimental and numerical curves. For undrained tests, the computed and measured stress–strain curves are also in good agreement (Figs. 7d–f). The effective stress paths indicate that, for the normally consolidated and slightly over-consolidated samples, the stress paths do not overpass the critical state line; whereas, for the strongly over-consolidated specimen, the stress path goes above the critical state line, at which dilation occurs, leading to an increase in the mean effective stress, before converging toward the critical state line. Overall, by using parameters calibrated from one drained test, the model is capable of predicting the stress–strain behavior of undrained tests for specimens with different OCRs.

### 3.1.3 Natural clays

Undrained triaxial tests on isotropically consolidated samples of clay mixture with different OCRs (OCR = 1.0, 1.6, and 4.0) were performed by Li and Meissner [20]. The main minerals of the clay mixture are kaolinite (60 %), illite (5 %), and quartz. The parameters presented in Table 1 were calibrated from the undrained triaxial test performed on the over-consolidated specimen (OCR = 4) with its consolidation stage. The determined parameters were used to simulate undrained triaxial tests on samples with OCR = 1 and 1.6. As shown in Fig. 8, the model predictions are in good agreement with the experimental results.

Numerical simulations of tests performed on samples of Boston blue clay (BBC) by Ladd and Varallyay [18] were also undertaken, as well as tests performed on samples of Lower Cromer Till (LCT) by Gens [11]. The database for both clays includes undrained triaxial tests on both isotropically and anisotropically consolidated specimens with different OCRs. Figure 9a, b shows the computed results of anisotropic compression and oedometer tests compared to

experimental data on BBC. Figures 9c, d and 10a, g, j show the comparisons between measured and predicted results of undrained triaxial tests on isotropically consolidated samples of LCT and BBC. For both natural clays, the present model gives good prediction for tests on lightly and heavily over-consolidated samples using the set of parameters determined from the tests on normally consolidated clay samples. The predictions for tests on anisotropically consolidated samples will be discussed later.

## 3.2 Influence of Lode's angle

### 3.2.1 Fujinomori clay

Drained triaxial tests on Fujinomori clay were performed at constant $p'$ in compression and extension on isotropically consolidated samples with different OCRs (OCR = 1, 2, 4, 8) by Nakai and Hinokio [23]. Drained true triaxial tests on normally consolidated Fujinomori clay samples were also performed at constant $p'$ with different Lode's angles (0°, 15°, 30°, 45°). The parameters presented in Table 1 were calibrated from one drained test in compression on a normally consolidated specimen with an isotropic consolidation stage.

Comparisons between experimental and numerical results for triaxial tests in compression and extension are presented in Fig. 11. One can see that the peak stress ratio of the over-consolidated clay samples and the amount of dilatancy increase with increased OCR. The model was able to capture the trend of the stress–strain behavior of clay with different OCRs under both compression and extension conditions.

The observed and predicted responses of normally consolidated Fujinomori clay in true triaxial tests with different Lode's angles are shown in Fig. 12. The stress ratio $\sqrt{3J_2}/p'$ is plotted against the principal strains (major strain in Fig. 12a, intermediate strain in Fig. 12b, minor strain in Fig. 12c, respectively). The predictions agree

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Comparison between experimental results and model predictions for undrained tests on clay mixture a void ratio versus mean effective stress for isotropic compression test, b deviatoric stress versus axial strain; c effective stress path</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 9** Comparison between experimental results and model predictions for tests on Boston blue clay a void ratio versus mean effective stress for anisotropic compression test, b void ratio versus vertical effective stress for oedometer test, c, e deviatoric stress versus axial strain for undrained tests; d, f effective stress path for undrained tests</figcaption></figure>

reasonably well with the experimental data in all cases. The comparisons demonstrate that the proposed model can take into account the influence of the intermediate principal stress on the stress–strain relationship.

### 3.2.2 Other natural clays

Undrained triaxial tests in extension were also carried out on Boston blue clay by Ladd and Varallyay [18] and on Lower Cromer Till by Gens [11]. Figures 9e, f and 10j, k show the comparisons between measured and predicted results of undrained triaxial tests on BBC and LCT. For both natural clays, the present model gives good prediction for tests in extension ($\theta = 60^\circ$) using parameters determined from the test in compression ($\theta = 0^\circ$).

## 3.3 Influence of consolidation stress ratio

### 3.3.1 Boston blue clay

Two undrained triaxial tests on anisotropically consolidated sample of BBC under compression and extension were also performed by Ladd and Varallyay [18]. The model parameters determined from tests on isotropically consolidated samples were used for the predictions.

Figure 9e, f shows a good agreement between the numerical and the experimental results for these two tests. It is interesting to point out that the undrained softening response following an anisotropic consolidation is reproduced by the model.

### 3.3.2 Lower Cromer Till

Drained and undrained triaxial tests in compression and extension with different OCRs on anisotropically consolidated samples of LCT were conducted by Gens [11].

Simulations were carried out to evaluate the model's performance for predicting the compression and extension tests on anisotropically consolidated samples (Fig. 10b, c). The samples were first anisotropically consolidated under $K_0 = 0.5$ up to $\sigma'_a = 350$ kPa. Then, they were unloaded along a different stress path to four different over-consolidation ratios (OCR = 1, 2, 4, 7) before being sheared in compression and extension under undrained condition. Figure 10b, c shows the comparisons between the experimental data and the model predictions. The comparisons indicate a good agreement in the major features of the undrained behavior for anisotropically consolidated samples.
<figure><img src="image_11.png" /><figcaption>**Fig. 10** Comparison between experimental results and model predictions for tests on LCT clay **a** void ratio versus mean effective stress for compression test, **b**, **c** undrained tests on anisotropically consolidated samples with different OCRs; **d–f** drained tests on anisotropically consolidated samples with different OCRs; **g–i** drained tests on anisotropically consolidated samples with OCR = 1; **j–l** undrained tests on anisotropically consolidated samples with OCR = 1</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11** Comparison between experimental results and model predictions for drained triaxial tests on Fujinomori clay with different OCRs
**a** stress ratio *q*/*p*' versus deviatoric strain under compression; **b** volumetric strain versus deviatoric strain under compression; **c** stress ratio *q*/*p*' versus deviatoric strain under extension; **d** volumetric strain versus deviatoric strain under extension</figcaption></figure>

Figure 10d–f shows the comparisons between experi-
mental data and model predictions of drained tests on
anisotropically consolidated specimens with different
OCRs. Similar to the undrained tests, the samples were first
anisotropically consolidated under K₀ = 0.5, then unloa-
ded to different values of OCR (1, 1.5, 2, 4, and 7) with

different values of K₀ (0.4, 0.5, 0.67, 1.0) along a different
stress path. Then, the samples were subjected to an axial
compression until a vertical strain of 15 %. Maximum
shear strength and volume change for samples with dif-
ferent values of OCR are well captured by the present
model using the same set of parameters as for undrained
<figure><img src="image_13.png" /><figcaption>**Fig. 12** Comparison between experimental results and model predictions for drained true triaxial tests on normally consolidated Fujinomori clay with different Lode's angle $a$ stress ratio $lql/p'$ versus major principle strain; $b$ stress ratio $lql/p'$ versus intermediate principle strain; $c$ stress ratio $lql/p'$ versus minor principle strain</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13** Determination of $c$ for the value of plastic hardening modulus according to the plasticity index</figcaption></figure>

tests. However, some discrepancies between experiments and simulations can be observed in the evolution of the deviatoric stress and volume change with the axial strain.

Four drained triaxial compression tests on normally consolidated samples were also selected for simulation. After being consolidated to different values of $K_0$ (0.4, 0.5, 0.67, 1.0), the samples were loaded to failure in drained condition. Figure 10g-i shows a good agreement between the numerical and the experimental results of the drained triaxial tests, using the set of parameters determined from an undrained test (Table 1). The measured volumetric strain increases with the value of $K_0$, which can be explained by the position of the compression lines in the $e-\log p'$ plane (Fig. 10i).

Undrained triaxial tests on normally consolidated samples with different consolidation stress ratios were also simulated. The samples were first anisotropically consolidated with four different consolidation stress ratios $K_0$ (i.e., the ratio of radial to axial stress $\sigma'_r/\sigma'_a$: 0.4, 0.5, 0.67 and 0.8. Then, for each $K_0$, two subsequent undrained shearing tests were conducted: one in compression (with an increase in the axial strain) and the other in extension (with a decrease in the axial strain). For all the eight loading paths mentioned above, Fig. 10j-k shows a good agreement

between the numerical and the experimental results using the set of parameters given in Table 1.

### 3.3.3 Softening response under undrained compression following anisotropic consolidation

A peculiar behavior worth to be noted is the softening response in undrained compression for the two cases with $K_0$ consolidation (see $K_0 = 0.4$ and 0.5 in Fig. 10). The same form of softening response has been observed in other types of clay (e.g., Fig. 9e, f on Boston blue clay by Ladd and Varallyay [18]). The measured softening response cannot be attributed to the destructuration process since the tested clays were reconstituted in the laboratory. This form of softening response is difficult to model by the conventional methods using kinematic hardening of the yield surface (e.g., Ling et al. [21]; Wheeler et al. [33]). However, using specific rotational kinematic hardening rules and yield-surface shapes, Pestana et al. [25] and Dafalias et al. [9] have managed to simulate the softening response after $K_0$ consolidation.

Different from the approach via kinematic hardening of the yield surface, the proposed approach employs the density state variable $e_c/e$. This assures the void ratio as well as the stress state to approach the critical state simultaneously, for any loading path. At large shear strains, the stress state converges toward the critical state. Thus, the magnitude of the mean effective stress $p'$ is governed by the location of the critical state line (in the $e-\log p'$ plane). The void ratio also approaches the critical state, thus the shear strength $q$ (on the $p'-q$ plane) is determined from $p'$ corresponding to the critical state $e_c$ (see the schematic plot in Fig. 10l). As a consequence, the softening response will occur when the deviatoric stress at the end of the $K_0$ consolidation is higher than the undrained shear strength determined from the critical state.

Overall, for all the examples selected in this study, the numerical simulations are in agreement with the
experimental results. For two parameters ($G_{p0}$ and $D$) dif-
ferent from Cam-Clay models, the value of $D$ varies from
0.3 to 1.2 and mostly around 1 (see Table 1); and the value
of $G_{p0}$ is suggested to be determined from the plasticity
index based on widely selected clays (see Fig. 13), thus is
not needed as input. Therefore, the proposed model
incorporating the density state controlling explicitly the
location of the critical state is simple in terms of parame-
ters determination, and is able to describe drained and
undrained behaviors of clay subjected to monotonic load-
ing after isotropic or anisotropic consolidation.

3.4 Cyclic loading

3.4.1 Fujinomori clay

The parameters determined from the drained compression
test under monotonic loading (see Table 1; Fig. 11) were
used to simulate three drained triaxial tests under cyclic
loading on Fujinomori clay performed by Nakai and
Hinokio [23]. The model predictions of the clay response
are compared with the experimental results: Fig. 14a, b
presents the results of a drained cyclic test under con-
stant confining stress; Fig. 14c, d presents the results of a
varying-amplitude cyclic test under constant mean

effective stress, in which the stress ratio increases with
the number of cycles; Fig. 14e, f presents the result of a
constant-amplitude cyclic test under constant mean
effective stress. For each test, the curves giving the
stress ratio versus the deviatoric strain and the stress
ratio versus the volumetric strain are plotted. All com-
parisons between the test results and the numerical
simulations demonstrate that the model can reasonably
well describe the cyclic behavior of clay in drained tri-
axial tests.

3.4.2 Clay mixture

Undrained triaxial tests were conducted on normally con-
solidated samples under cyclic loading by Li and Meissner
[20]. The cyclic loading program involves one- and two-
way cyclic tests. All the cyclic tests were stress controlled,
with a sinusoidal wave form at frequency of 0.1 Hz. The
cyclic stress ratio, defined as the ratio of the applied cyclic
shear stress to the monotonic shearing strength in com-
pression, ranges from 0.50 to 0.80.

The set of parameters determined from the monotonic
tests were used to simulate the undrained cyclic tests. Note
that the initial slopes of the $q - \varepsilon_d$ curve for monotonic,
one- and two-way cyclic tests are different from each

<figure><img src="image_15.png" /><figcaption>**Fig. 14** Comparison between experimental results and model predictions for drained triaxial tests under cyclic loading on normally consolidated Fujinomori clay with different stress paths a, b stress ratio *q*/*p'* versus deviatoric and volumetric strains for test at constant confining stress; c, d stress ratio *q*/*p'* versus deviatoric and volumetric strains for test at constant *p'* with increasing *q*/*p'*; e, f stress ratio *q*/*p'* versus deviatoric and volumetric strains for test at constant *p'* with constant *q*/*p'*</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 15** Comparison between experimental results and model predictions for undrained tests on clay mixture a deviatoric stress versus deviatoric strain for one-way cyclic loading; b pore water pressure versus number of cycles for one-way cyclic loading; c deviatoric stress versus deviatoric strain for two-way cyclic loading; d pore water pressure versus number of cycles for two-way cyclic loading</figcaption></figure>
others due to different loading rates. Since the rate-dependency behavior of clay was not considered in the elastoplastic model, different values of $G_{p0}$ were selected for each case ($G_{p0} = 670$ for one-way test, $G_{p0} = 5,000$ for two-way test), as done by Li and Meissner [20] and Yu et al. [41]. Simulations are given in Fig. 15a, b for one-way cyclic tests and in Fig. 15c, d for two-way cyclic tests and compared with experimental results. The excess pore pressure ($\Delta u$) and the deviatoric stress are plotted as a function of the number of cycles and deviatoric strain, respectively. One can see that the proposed model can capture with reasonable accuracy, the undrained behavior of clay subjected to one- and two-way cyclic loading.

## 4 Conclusions

A simple critical-state-based double-yield-surface model was developed for describing the mechanical behavior of clay. The model has two yield surfaces: one for shear sliding and one for compression. A CSL-related density state was defined to link the peak strength and the phase transformation characteristics to the material's void ratio. Therefore, the strain-hardening behavior with contraction for normally consolidated clay and the strain-softening behavior with dilation for over-consolidated clay can be modeled. The model guarantees also that stresses and void ratio reach simultaneously the critical state line in the $p'-q-e$ space. The stress reversal technique was incorporated into the model for describing the mechanical behavior of clay under loading with changes in stress direction (for instance, anisotropic consolidation followed by monotonic loading/unloading, cyclic loadings). The model has 6 material parameters and 2 state variables which can be easily determined based on one drained or undrained triaxial test up to failure with an isotropic consolidation stage.

The capability of the model to reproduce the main features of clay behavior was examined by comparing experimental results and numerical simulations of drained and undrained triaxial tests under monotonic loading with different conditions (different OCRs, different Lode's angles, different consolidation stress ratios) on various remolded clays (Black clay and White clay) and natural clays (Fujinomori clay, Boston blue clay, Lower Cromer Till, and clay mixture). Drained and undrained tests under cyclic loadings were also simulated by using the sets of parameters determined from monotonic tests.

All comparisons between experimental results and numerical simulations demonstrate that the proposed model is capable of reproducing the behavior of clays with different stress histories, different stress paths, different drainage conditions and different loading conditions.

**Acknowledgments** This research was financially supported by the opening project of the State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (Grant No. SKLGP2013K025), the National Natural Science Foundation of China (Grant No. 41240024), the Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20110073120012), and the Shanghai Pujiang Talent Plan (Grant No. 11PJ1405700).

## References

1. Balendran B, Nemat-Nasser S (1993) Double sliding model for cyclic deformation of granular materials including dilatancy effects. J Mech Phys Solids 41(3):573-612

2. Been K, Jefferies MG (1985) A state parameter for sands. Geotechnique 35(2):99-112

3. Biarez J, Hicher PY (1994) Elementary mechanics of soil behaviour, Balkema

4. Bryson LS, Salehian A (2011) Performance of constitutive models in predicting behavior of remolded clay. Acta Geotech 6:143-154

5. Calladine CR (1971) Microstructural view of the mechanical properties of saturated clay. Geotechnique 21(4):391-415

6. Chang CS, Yin Z-Y (2010) Modeling stress-dilatancy for sand under compression and extension loading conditions. ASCE J Eng Mech 136(6):777-786

7. Chang CS, Yin Z-Y (2010) Micromechanical modelling for inherent anisotropy in granular materials. ASCE J Eng Mech 136(7):830-839

8. Chowdhury EQ, Nakai T (1998) Consequence of the tij-concept and a new modelling approach. Comput Geotech 23(4):131-164

9. Dafalias YF, Manzari MT, Papadimitriou AG (2006) SANI-CLAY: simple anisotropic clay plasticity model. Int J Numer Anal Methods Geomech 30(12):1231-1257

10. Gajo A, Muir Wood D (1999) Severn-Trent sand: a kinematic-hardening constitutive model: the $q-p$ formulation. Geotechnique 49(5):595-614

11. Gens A (1982) Stress-strain and strength of a low plasticity clay. Ph.D. Thesis at Imperial College, London University

12. Hardin BO (1978) The nature of stress-strain behaviour of soils. In: Proceedings of the earthquake engineering and soil dynamics, vol 1, Pasadena, USA, pp 3-90

13. Hirayama H (1987) Interpretation of the cam-clay model as a simplified double-yield-surface model. Soils Found 27(3):105-111

14. Hsieh HS, Kavazanjian JE, Borja RI (1990) Double-yield-surface Cam-clay plasticity model. I: theory. J Geotech Eng 116(9): 1381-1401

15. Huang WX, Wu W, Sun DA, Scott S (2006) A simple hypoplastic model for normally consolidated clay. Acta Geotech 1:15-27

16. Hujeux JC (1985) Une loi de comportement pour le chargement cyclique des sols. In: Davidovici V (ed) Génie Parasismique. Presses ENPC, France, pp 278-302

17. Jefferies MG (1993) NorSand: a simple critical state model for sand. Geotechnique 43(1):91-103

18. Ladd CC, Varallyay J (1965) The influence of the stress system on the behaviour of saturated clays during undrained shear. Research Rep. No. R65-11, Department of Civil Engineering, MIT, Cambridge, MA

19. Lade PV (2007) Modeling failure in cross-anisotropic frictional materials. Int J Solids Struct 44(16):5146-5162

20. Li T, Meissner H (2002) Two-surface plasticity model for cyclic undrained behavior of clays. ASCE J Geotech Geoenviron Eng 128(7):613-626
21. Ling HI, Yue D, Kaliakin VN (2002) Anisotropic elastoplastic bounding surface model for cohesive soils. ASCE J Eng Mech 128(7):748-758

22. Mašin D (2005) A hypoplastic constitutive model for clays. Int J Numer Anal Methods Geomech 29(4):311-336

23. Nakai T, Hinokio M (2004) A simple elastoplastic model for normally and overconsolidated soils with unified material parameters. Soils Found 44(2):53-70

24. Ohmaki S (1979) A mechanical model for the stress-strain behaviour of normally consolidated cohesive soil. Soils Found 19(3):29-44

25. Pestana JM, Whittle AJ, Gens A (2002) Evaluation of a constitutive model for clays and sands: Part II—Clay behaviour. Int J Numer Anal Methods Geomech 26(11):1123-1146

26. Roscoe KH, Burland JB (1968) On the generalized stress-strain behavior of 'wet' clay. Engineering Plasticity. Cambridge University Press, Cambridge, pp 553-609

27. Schofield AN, Wroth CP (1968) Critical state soil mechanics. Mcgraw Hill, Maidenhead

28. Schweiger H, Wiltafsky C, Scharinger F, Galavi V (2009) A multilamine framework for modelling induced and inherent anisotropy of soils. Geotechnique 59(2):87-101

29. Sheng D, Sloan SW, Yu HS (2000) Aspects of finite element implementation of critical state models. Comput Mech 26:185-196

30. Stallebrass SE, Taylor RN (1997) The development and evaluation of a constitutive model for the prediction of ground movements in overconsolidated clay. Geotechnique 47(2):235-253

31. Taiebat M, Dafalias YF (2008) SANISAND: simple anisotropic sand plasticity model. Int J Numer Anal Methods Geomech 32(8):915-948

32. Vermeer PA (1978) A double hardening model for sand. Geotechnique 28(4):414-433

33. Wheeler SJ, Näätänen A, Karstunen M, Lojander M (2003) An anisotropic elasto-plastic model for soft clays. Can Geotech J 40:403-418

34. Whittle AJ, Kavvadas MJ (1994) Formulation of MIT-E3 constitutive model for overconsolidated clays. J Geotech Eng 120(1):173-198

35. Yamakawa Y, Hashiguchi K, Ikeda K (2010) Implicit stress-update algorithm for isotropic Cam-clay model based on the subloading surface concept at finite strains. Int J Plast 26(5):634-658

36. Yao YP, Hou W, Zhou AN (2009) UH model: three-dimensional unified hardening model for overconsolidated clays. Geotechnique 59(5):451-469

37. Yin Z-Y, Chang CS (2009) Microstructural modelling of stress-dependent behaviour of clay. Int J Solids Struct 46(6):1373-1388

38. Yin Z-Y, Chang CS (2013) Stress-dilatancy for sand under loading and unloading conditions. Int J Num Anal Methods Geomech. doi:10.1002/nag.1125

39. Yin Z-Y, Chang CS, Hicher PY, Karstunen M (2009) Micromechanical analysis of kinematic hardening in natural clay. Int J Plast 25(8):1413-1435

40. Yin Z-Y, Chang CS, Hicher PY (2010) Micromechanical modelling for effect of inherent anisotropy on cyclic behaviour of sand. Int J Solids Struct 47(14-15):1933-1951

41. Yu H-S, Khong C, Wang J (2007) A unified plasticity model for cyclic behaviour of clay and sand. Mech Res Commun 34:97-114

42. Zervoyannis C (1982) Etude synthetique des proprietes mecaniques des argiles et des sables sur chemins oedometrique et triaxial de revolution. Ph.D. thesis, Ecole Centrale de Paris