A state-dependent soil model and its application to principal stress rotation simulations

Abstract

The plastic strain caused by principal stress rotation is one of the most important factors contributing to substantial deformation under earthquake, wave or traffic loading. The original Pastor–Zienkiewicz Mark III model, a well-known model for the analysis of the dynamic response under cyclic loading, is unable to consider the effects of principal stress orientation as well as state-dependent dilatancy. In this article, a new constitutive model for sand is developed to consider both aforementioned effects based on the original Pastor–Zienkiewicz Mark III model. There are 14 model parameters in total for the static condition and three extra parameters for cyclic loading, and a corresponding calibration method of model parameters is proposed. The predictive capability of the proposed model is verified with the results of a series of experiments on sand, including undrained monotonic tests in different fixed principal stress orientations and undrained cyclic rotational shear tests. The comparisons indicate that the proposed model can effectively incorporate the effects of principal stress orientation and state-dependent dilatancy.

Introduction

The phenomenon of principal stress rotation is very common in sea-floor sediments under wave loading and foundations under earthquake loading or traffic loading.^1,2 However, models formulated under the traditional plasticity theory in principal stress space cannot reflect the principal stress rotation effect. For example, such models indicate that plastic strain will not be produced by the pure principal stress rotation when the magnitude of the cyclic deviatoric stress remains unchanged, but the direction of principal stress rotates progressively. Many experimental results^3–5 have shown that either the principal stress rotation in cyclic rotational shear tests or the fixed principal stress orientation variation in monotonic loading tests has a significant impact on stress–strain behaviour of sand. These effects are collectively referred to as the effects of principal stress orientation.

Only a few soil constitutive models in the literature reflect the effects of principal stress orientation. Sassa and Sekiguchi ⁶ developed a new model (the Sassa model) to consider principal stress rotation under two-dimensional (2D) plane strain conditions based on the Pastor–Zienkiewicz Mark III (PZ3) model by introducing the principal stress angle. In recent years, new models and methods^7–9 have been proposed to investigate the principal stress rotation and the property of anisotropy that are considered to be the most likely reason for the effects of principal stress orientation. In general, the Sassa model is simpler in practice and employs fewer equations and parameters that are relatively easily formulated and calibrated. However, there is still room for improvements and enhancements. First, the Sassa model can only be used in simple 2D plane strain problems, as the out-of-plane stress, which could be important in determining the plastic flow condition, was neglected. Second, the effects of unloading and reloading are not considered for principal stress rotation. Finally, the Sassa model is not able to consider the effects of state-dependent dilatancy directly because it is formulated for a single initial void ratio. Thus, the Sassa model requires different sets of model parameters for the same type of sand consisting of different relative densities.

The theory of state-dependent dilatancy was proposed by Wood et al. ¹⁰ and expanded upon by Manzari and Dafalias, ¹¹ Li and Dafalias ¹² and Li. ¹³ The theory states that the constitutive behaviour of sand is closely related to state parameter, which depends on the current physical state such as void ratio and stress states, including confining pressure. Although the theory of state-dependent dilatancy has been recognized and some state-dependent constitutive models have been proposed, most of them have not considered the effects of the principal stress orientation.

In this article, the PZ3 model was used as the base model to develop a new constitutive model for sand following the hierarchical approach proposed by Desai. ¹⁴ The PZ3 model developed by Pastor et al.^15,16 based on the generalized plasticity theory ¹⁷ is a well-known model for the analysis of the dynamic response under cyclic loading such as that due to earthquakes and wave loading. In the proposed model, additional features are included hierarchically by introducing further hardening rules into the PZ3 model while preserving all the features of the base model such as many validated, relevant stress paths when the new feature is not active. This approach avoids a complete reformulation of the constitutive models and only requires the additional features to be validated. This method has proven to be successful and effective by Manzanal et al., ¹⁸ who developed a new model that can consider the dilatancy of sand based on the PZ3 model.

Following this course of studies and within the framework of generalized plasticity, a new constitutive model to reflect both the effects of principal stress orientation and state-parameter dependency is proposed and formulated in section ‘Model description’. The calibration method for model parameters is proposed. The predictive capacity of this new model is verified through comparison with the results of fixed principal stress orientation monotonic tests ³ and cyclic rotational shear tests^5,19 and reported in section ‘Model evaluation’ together with further discussion. The key conclusions of this article are then given in section ‘Discussion and conclusion’.

Model description

Before establishing the constitutive model, the principal stress orientation and the stress variables in the proposed model are defined in sections ‘Definition of the principal stress orientation’ and ‘Inclusion of three stress invariants’, respectively. Then, the theory of state-dependent dilatancy is introduced in section ‘Introduction of the state-dependent dilatancy’. The basic equations of the model are given in section ‘Basic equations’, and the calibration method is illustrated in section ‘Model calibration’.

Definition of the principal stress orientation

The principal stress orientation can be defined in the 2D plane $[\frac{{\sigma \prime }_z-{\sigma \prime }_x}{2}-{\tau }_{\mathrm{zx}}]$ using the angle of the major principal stress orientation from the vertical axis as shown in Figure 1, which can be expressed as

$$\tan \left(2\psi \right)=\frac{2{\tau }_{\mathrm{zx}}}{{\sigma \prime }_z-{\sigma \prime }_x}$$

(1)

Figure 1.

Principal stress angle $\psi$ .

Inclusion of three stress invariants

The intermediate principal stress was not considered in the Sassa model. However, the coefficient of intermediate principal stress $b=({\sigma }_2-{\sigma }_3)/({\sigma }_1-{\sigma }_3)$ is a very important control factor in the experiments; therefore, the influence of the Lode angle should not be neglected. In this article, the mean effective stress $p\prime$ and generalized shear stress q are used in the new model. Meanwhile, the Lode angle $\theta$ is also considered; therefore

$$p^{\prime }=\frac{1}{3}({{\sigma }^{\prime }}_x+{{\sigma }^{\prime }}_y+{{\sigma }^{\prime }}_z)=\frac{1}{3}{I}_1$$

(2)

$$\begin{array}{l}q=\sqrt{3J_2}=\frac{1}{\sqrt{2}}({({{\sigma }^{\prime }}_x-{{\sigma }^{\prime }}_y)}^2+{({{\sigma }^{\prime }}_x-{{\sigma }^{\prime }}_y)}^2\\ +{({{\sigma }^{\prime }}_x-{{\sigma }^{\prime }}_y)}^2+6({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2){)}^{\frac{1}{2}}\end{array}$$

(3)

$$\sin 3\theta =-\frac{3\sqrt{3}}{2}\frac{J_3}{{\left(J_2\right)}^{\frac{3}{2}}}=-\frac{27J_3}{2q^3}$$

(4)

where $J_3=({{\sigma }^{\prime }}_x-p^{\prime })·({{\sigma }^{\prime }}_y-p^{\prime })·({{\sigma }^{\prime }}_z-p^{\prime })+2{\tau }_{xy}{\tau }_{yz}{\tau }_{zx}$ $-(\sigma {\prime }_x-p\prime ){\tau }_{\mathrm{yz}}^2-(\sigma {\prime }_y-p\prime ){\tau }_{\mathrm{zx}}^2-(\sigma {\prime }_z-p\prime ){\tau }_{\mathrm{xy}}^2$ , here, $I_1$ , $J_2$ and $J_3$ are the three stress invariants. $p\prime$ , q and $\theta$ are the functions of $I_1$ , $J_2$ and $J_3$ , respectively; therefore, they are also stress invariants.

Introduction of the state-dependent dilatancy

The theory of state-dependent dilatancy has been successfully used in the modelling of sand behaviour. In the theory of constitutive model for sands, one of the fundamental issues is to describe dilatancy d correctly. Rowe ²⁰ suggested that the dilatancy d was simply a unique function of the stress ratio

$$d=d(\eta )=M-\eta$$

(5)

where M is the critical stress ratio. However, it was soon found that the dilatancy was not only related with stress ratio $\eta =q/p\prime$ but also with the material internal state. Based on the aforementioned observations and the critical state constitutive framework, Li and Dafalias ¹² proposed a new expression for dilatancy, which was the function of the stress ratio $\eta$ and state parameter $\phi$ , expressed as

$$d=d(\eta ,\phi )=d_0(e^{m\phi }-\frac{\eta }{M})$$

(6)

where $\phi$ is a state parameter to reflect the current state of sand and m is a material parameter. State parameter $\phi$ was defined by Been and Jefferies. ²¹ Based on the linear representations of the steady-state line for sand by Li and Wang ²² shown in Figure 2, state parameter $\phi$ can be expressed as

$$\phi =e-e_{\mathrm{c}}=e-(e_{\Gamma }-{\lambda }_{\mathrm{c}}{\left(\frac{p\prime }{p_{\mathrm{a}}}\right)}^{\kappa })$$

(7)

where e is the current void ratio; $e_{\mathrm{c}}$ is the critical void ratio, which is a function of the confining pressure $p\prime$ ; and model parameters $e_{\Gamma }$ , ${\lambda }_{\mathrm{c}}$ , $\kappa$ and $p_{\mathrm{a}}$ represent atmospheric pressure.

Figure 2.

Definition of state parameter.

Basic equations

Within the framework of generalized plasticity, the yield surface and plastic potential need not be explicitly defined. Instead, the loading direction vector $\mathbf{n}$ and the plastic flow direction vector ${\mathbf{n}}_{\mathbf{g}}$ are used. With appropriate laws, which can be obtained by tests for the plastic flow, loading direction and plastic moduli, some salient or special behaviours of sand can be described.

In the original PZ3 model, the loading direction vector $\mathbf{n}=(n_{\mathrm{v}},n_{\mathrm{s}}{)}^{\mathrm{T}}$ and the plastic flow ${\mathbf{n}}_{\mathbf{g}}=(n_{\mathrm{gv}},n_{\mathrm{gs}}{)}^{\mathrm{T}}$ are defined. Based on the PZ3 model, the vectors $\mathbf{n}$ and ${\mathbf{n}}_{\mathbf{g}}$ are developed and expanded to reflect the influence of the Lode angle $\theta$ , the orientation of principal stress $\psi$ and the state parameter $\phi$ on the basis of experimental regularities.

The loading direction vector $\mathbf{n}$ can be defined as follows

$$\mathbf{n}={(n_{\mathrm{v}},n_{\mathrm{s}},n_{\theta },n_{\psi },n_{\phi })}^{\mathrm{T}}$$

(8)

The plastic flow direction vector ${\mathbf{n}}_{\mathbf{g}}$ can be defined as follows

$${\mathbf{n}}_g={(n_{\mathrm{gv}},n_{\mathrm{gs}},n_{g\theta },n_{g\psi },n_{g\phi })}^{\mathrm{T}}$$

(9)

The incremental stress tensor can be expressed as

$$\mathrm{d}\sigma {\prime }_{\mathrm{ij}}=({\mathbf{D}}_{\mathbf{ijkl}}^{\mathbf{e}}-\frac{{\mathbf{D}}_{\mathbf{ijmn}}^{\mathbf{e}}:{\mathbf{n}}_{\mathbf{g}}\otimes \mathbf{n}:{\mathbf{D}}_{\mathbf{stkl}}^{\mathbf{e}}}{H_{\mathrm{L},\mathrm{U}}+\mathbf{n}:{\mathbf{D}}_{\mathbf{stkl}}^{\mathbf{e}}:{\mathbf{n}}_{\mathbf{g}}}):\mathrm{d}{\epsilon }_{\mathrm{kl}}$$

(10)

where $D_{\mathrm{ijkl}}^e$ is the elastic stiffness tensor (detailed in Appendix 1).

In the original PZ3 model, $n_{\mathrm{v}}$ and $n_{\mathrm{s}}$ are assumed to take the following forms

$$n_{\mathrm{v}}=(1+\alpha )(M_{\mathrm{f}}-\eta );n_{\mathrm{s}}=1$$

(11)

and $n_{\mathrm{gv}}$ and $n_{\mathrm{gs}}$ are assumed to take the following forms

$$n_{\mathrm{gv}}=(1+\alpha )(M_{\mathrm{g}}-\eta );n_{\mathrm{gs}}=1$$

(12)

here, the critical stress ratio $M_{\mathrm{g}}$ can be obtained by an undrained triaxial test, and the model parameter $M_{\mathrm{f}}$ is calculated by a simple approximate relation where $M_{\mathrm{f}}=D_{\mathrm{r}}·M_{\mathrm{g}}$ . $D_{\mathrm{r}}$ is the relative density. This assumption leads to the limitation that the PZ3 model regards sands of the same type with different initial relative densities and confining pressures as different materials, which means that the PZ3 model requires different parameters for a single type of sand with different relative density under different confining pressures. Furthermore, the principal stress orientation is not considered.

In this study, the state parameter $\phi$ and principal stress orientation $\psi$ are both introduced into the proposed model based on the following experimental evidence and the theory of state-parameter dependency:

1. Contractive soil behaviour continues until the state of phase transformation becomes more marked with increasing $\phi$ .

2. Contractive soil behaviour continues until the state of phase transformation becomes more marked with increasing $\psi$ .

3. Slope of the phase transformation line decreases with increasing $\phi$ .

4. Slope of the phase transformation line decreases with increasing $\psi$ .

5. Critical state is unique.

Based on these experimental results, $n_{\mathrm{gv}}$ and $n_{\mathrm{gs}}$ are redefined in the following forms

$$n_{\mathrm{gv}}=(1+\alpha \left(\psi \right))(M_{\mathrm{g}}(\theta ,\psi ,\phi )-q/p\prime );n_{\mathrm{gs}}=1$$

(13)

where

$$M_{\mathrm{g}}(\theta ,\psi ,\phi )=M_{\mathrm{g}0}(\theta ,\phi )-U\left(\psi \right)·a·M_{\mathrm{g}0}(\theta ,\phi )$$

(14)

$$M_{\mathrm{g}0}(\theta ,\phi )=\frac{18}{18+3(1-\sin 3\theta )}M{e}^{m\phi }$$

(15)

$$\alpha \left(\psi \right)={\alpha }_0+g·U\left(\psi \right)$$

(16)

$$U\left(\psi \right)=\frac{1-\cos (2\psi )}{2}$$

(17)

here, ${\alpha }_0$ , a, g and m are the model parameters. M is the critical stress ratio.

$n_{\mathrm{v}}$ and $n_{\mathrm{s}}$ are redefined in the following forms

$$n_{\mathrm{v}}=(1+\alpha \left(\psi \right))(M_{\mathrm{f}}(\theta ,\psi ,\phi )-\frac{q}{p\prime });n_{\mathrm{s}}=1$$

(18)

where

$$M_{\mathrm{f}}(\theta ,\psi ,\phi )=M_{\mathrm{f}0}(\theta ,\phi )-U\left(\psi \right)·a·M_{\mathrm{f}0}(\theta ,\phi )$$

(19)

$$M_{\mathrm{f}0}(\theta ,\phi )=\frac{18}{18+3(1-\sin 3\theta )}M{e}^{-n\phi }$$

(20)

here, n is a model parameter.

Then, the loading direction tensor $\mathbf{n}$ can be obtained as follows

$$\{\begin{array}{c}{n}_{\mathrm{v}}=\left(1+\alpha \left(\psi \right)\right)\left(M_{\mathrm{f}}\left(\theta ,\psi ,\phi \right)-q/p^{\prime }\right)\\ n_s=1\\ n_{\theta }=-\frac{q}{M_{\mathrm{f}0}\left(\theta ,\phi \right)}\frac{\partial M_{\mathrm{f}0}\left(\theta ,\phi \right)}{\partial \theta }\\ n_{\psi }=U^{\prime }(\psi )·\\ \left\{a·M_{\mathrm{f}0}\left(\theta ,\phi \right)·p^{\prime }·\left(1+\frac{1}{\alpha \left(\psi \right)}\right)+\frac{g}{\alpha \left(\psi \right)}\left(\frac{q}{1+\alpha \left(\psi \right)}+X_2·\ln \frac{X_2}{X_1}\right)\right\}\\ n_{\phi }=-\frac{q}{M_{\mathrm{f}0}\left(\theta ,\phi \right)}\frac{\partial M_{\mathrm{f}0}\left(\theta ,\phi \right)}{\partial \phi }\end{array}$$

(21)

The details of $X_1$ and $X_2$ can be found in equations (41) and (42) of Appendix 1.

The plastic flow direction vector ${\mathbf{n}}_{\mathbf{g}}$ can be obtained as follows

$$\{\begin{array}{c}{n}_{\mathrm{gv}}=\left(1+\alpha \left(\psi \right)\right)\left(M_{\mathrm{g}}\left(\theta ,\psi ,\phi \right)-q/p^{\prime }\right)\\ n_{\mathrm{gs}}=1\\ n_{\mathrm{g}\theta }=-\frac{q}{M_{\mathrm{g}0}\left(\theta ,\phi \right)}\frac{\partial M_{\mathrm{g}0}\left(\theta ,\phi \right)}{\partial \theta }\\ n_{\mathrm{g}\psi }=U^{\prime }(\psi )\\ ·\{a{M}_{\mathrm{g}0}\left(\theta ,\phi \right)·p^{\prime }·(1+\frac{1}{\alpha \left(\psi \right)})+\frac{g}{\alpha \left(\psi \right)}(\frac{q}{1+\alpha \left(\psi \right)}+Y_2·\ln \frac{Y_2}{Y_1})\}\\ n_{\mathrm{g}\phi }=-\frac{q}{M_{\mathrm{g}0}\left(\theta ,\phi \right)}\frac{\partial M_{\mathrm{g}0}\left(\theta ,\phi \right)}{\partial \phi }\end{array}$$

(22)

The details of $Y_1$ and $Y_2$ can be found in equations (43) and (44) of Appendix 1.

The plastic modulus for loading should also depend on $\psi$ and $\phi$

$$\begin{array}{l}{H}_{\mathrm{L}}=H_0·p^{\prime }{\left(1-\frac{\eta \left(\theta ,\psi ,\phi \right)}{{\eta }_{\mathrm{f}}^{\ast }}\right)}^4\\ \left(1-\frac{q/p^{\prime }}{M_{\mathrm{g}}\left(\theta ,\psi ,\phi \right)}+{\beta }_0{\beta }_1\exp \left(-{\beta }_0\xi \right)\right)·D_{\mathrm{M}}\end{array}$$

(23)

where $H_0$ , ${\beta }_0$ and ${\beta }_1$ are the model parameters and

$$\eta (\theta ,\psi ,\phi )=\frac{q}{p\prime }+(1-U\left(\psi \right))·a·M_{g0}(\theta ,\phi )$$

(24)

$${\eta }_{\mathrm{f}}^{*}=(M_{\mathrm{f}0}\left(\phi \right)-a·M_{\mathrm{g}0}\left(\phi \right))(1+\frac{1}{a_0+g})$$

(25)

$$\xi =\int \mathrm{d}\xi =\int \left|\mathrm{d}{\epsilon }_{\mathrm{q}}^{\mathrm{p}}\right|$$

(26)

here, $\xi$ is the accumulated deviatoric plastic strain.

Regarding the reloading process, it is necessary to consider the history of past events. Thus, a discrete memory factor $D_{\mathrm{M}}$ was introduced as follows

$$D_{\mathrm{M}}={\left(\frac{{\zeta }_{max}}{\zeta }\right)}^{{\gamma }_{\mathrm{d}}}$$

(27)

where $\zeta$ is defined as

$$\zeta =p\prime {(1-\left(\frac{1+\alpha \left(\psi \right)}{\alpha \left(\psi \right)}\right)·\frac{q/p\prime }{M_{\mathrm{g}}(\psi ,\phi )})}^{1/\alpha \left(\psi \right)}$$

(28)

and ${\gamma }_{\mathrm{d}}$ is the model parameter.

The plastic modulus $H_{\mathrm{U}}$ for unloading is assumed to be

$$H_{\mathrm{U}}=\{\begin{array}{c}{H}_{\mathrm{U}0}{\left(\frac{M_{\mathrm{g}}(\psi ,\phi )}{q/p\prime }\right)}^{{\gamma }_u}\mathrm{for}\left|\frac{M_{\mathrm{g}}(\psi ,\phi )}{q/p\prime }\right|>1\hfill \\ \multicolumn{1}{c}{H_{\mathrm{U}0}\mathrm{for}\left|\frac{M_{\mathrm{g}}(\psi ,\phi )}{q/p\prime }\right|\le 1}\end{array}$$

(29)

where $H_{\mathrm{U}0}$ and ${\gamma }_{\mathrm{u}}$ are the original model parameters.

Model calibration

There are 14 model parameters in total for static conditions and 3 extra parameters for cyclic conditions with a, g, m and n introduced on top of the parameters of the original PZ3 model. It should be noted that when $a=g=0$ , there will be no effect due to principal stress rotation and when $n=m=0$ , there will be no effect due to state-dependent dilatancy. Therefore, when those four parameters are equal to zero at the same time, the current model will reduce to the original PZ3 model, thus confirming the hierarchical approach.

a, g, m, n, $e_{\Gamma }$ , ${\lambda }_{\mathrm{c}}$ , $\kappa$ and M are eight parameters not found in the original PZ3 model. $e_{\Gamma }$ , ${\lambda }_{\mathrm{c}}$ , $\kappa$ and M can be determined directly by fitting the triaxial test data for the critical state line in the $[e-p\prime ]$ plane.

Parameter m can be calculated by the undrained monotonic loading shear test at a principal stress angle of $\psi =0^{°}$ and a Lode angle $\theta =0^{°}$ (coefficient of intermediate principal stress $b=0.5$ ). In the phase transformation state, the dilatancy d is equal to zero, as follows

$$d=\frac{\partial g/\partial p}{\partial g/\partial q}=(1+{\alpha }_0)(M_{\mathrm{g}}(0,0,{\phi }^{*})-{\eta }^{*})=0$$

(30)

here, the state parameter and stress ratio denoted by superscript * represent measured values by the corresponding experiments.

Both $a_0$ and g parameters are positive; therefore

$$M_{\mathrm{g}}(0,0,{\phi }^{*})={\eta }^{*}$$

(31)

Because the Lode angle $\theta =0^{°}$ , using equations (14) and (15), the parameter m can be obtained as

$$m=\frac{1}{{\phi }^{*}}\ln \frac{7}{6M}·{\eta }^{*}$$

(32)

a and g are the parameters of the principal stress orientation. After m is obtained, a can be calculated through the undrained monotonic loading shear test at a principal stress angle of $\psi ={45}^{°}$ and Lode angle $\theta =0$ . In the phase transformation state, the dilatancy d is equal to zero

$$d=\frac{\partial g/\partial p}{\partial g/\partial q}=(1+{\alpha }_0+g)(M_{\mathrm{g}}(0,{45}^{°},{\phi }^{*})-{\eta }^{*})=0$$

(33)

Hence

$$M_{\mathrm{g}}(0,{45}^{°},{\phi }^{*})={\eta }^{*}$$

(34)

After m is obtained, the parameter a can be obtained using equation (14) as follows

$$a=2(1-\frac{7{\eta }^{*}}{6M{e}^{m{\phi }^{*}}})$$

(35)

${\alpha }_0$ and g can be calculated through drained monotonic loading shear tests. Ignoring the small elastic deformations, we get

$$\frac{\mathrm{d}{\epsilon }_{\mathrm{v}}}{\mathrm{d}{\epsilon }_{\mathrm{q}}}\approx \frac{\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{p}}}{\mathrm{d}{\epsilon }_{\mathrm{q}}^{\mathrm{p}}}=d=\frac{\partial g/\partial p}{\partial g/\partial q}=(1+\alpha (\psi ))(M_{\mathrm{g}}(\theta ,\psi ,\phi )-\frac{q}{p\prime })$$

(36)

Therefore, from equation (30), the parameter ${\alpha }_0$ can be identified using the dilatancy d at principal stress angle $\psi =0^{°}$ and Lode angle $\theta =0^{°}$ . For example, choosing the dilatancy $d^{\mathrm{p}}$ at the peak stress ratio, equation (30) yields

$${\alpha }_0=\frac{d^{\mathrm{p}}}{M_{\mathrm{g}}(0,0,{\phi }^{*})-{\eta }^{*}}-1$$

(37)

After ${\alpha }_0$ is calculated, the parameter g can be obtained using the same method, but at $\psi ={45}^{°}$ , equation (33) yields

$$g=2·(\frac{d^{\mathrm{p}}}{M_{\mathrm{g}}(0,0,{\phi }^{*})-{\eta }^{*}}-1-{\alpha }_0)$$

(38)

The parameter n cannot be easily calibrated directly and must be determined by fitting the stress–strain curve.

In addition to ${\alpha }_0$ , there are still five model parameters that are shared with the original PZ3 model under static conditions. Ling and Liu ²³ suggested a method of calibrating these five model parameters. $H_0$ can be determined by fitting the curves of $p\prime$ and q versus the axial strain, which can also be found by matching the plot of q versus $p\prime$ for undrained tests. ${\beta }_0$ and ${\beta }_1$ can be determined by matching the stress–strain curves, but ${\beta }_0$ is generally suggested to take the value of 4.2 and ${\beta }_1$ of 0.2 by Pastor et al. ¹⁶ The elastic parameters $G_0$ and $\nu$ can be obtained by observing the initial behaviour of the stress–strain relationship.

There are extra three model parameters, $H_{\mathrm{U}0}$ , ${\gamma }_{\mathrm{u}}$ and ${\gamma }_{\mathrm{d}}$ , for cyclic conditions. $H_{\mathrm{U}0}$ is determined by fitting the results to the initial slope of the first unloading stress–strain curve; ${\gamma }_{\mathrm{u}}$ is obtained by fitting the first unloading stress–strain curve. ${\gamma }_{\mathrm{d}}$ can be obtained by matching both the hysteretic loops in the stress–strain curve and $[q-{\epsilon }_{\mathrm{v}}]$ curve.

Model evaluation

To evaluate the predictive capability of the proposed model, one element response was analysed with the new model implemented into the finite element method (FEM) platform of DIANA–SWANDYNE II, ²⁴ and a series of simulations were carried out to compare with published experimental results in different stress conditions.

Monotonic loading, different fixed principal stress orientation

Nakata et al. ³ performed a series of monotonic loading experiments using the torsional cylinder shear apparatus to investigate the undrained deformation behaviour of Toyoura sand ( $G_s=2.643,e_{\max }=0.973,e_{\min }=0.635$ ) with different relative densities subjected to the change of principal stress orientation. In this study, three sets of tests in different principal stress angles $\psi ={15}^{°},{30}^{°},{45}^{°},{60}^{°},{75}^{°}$ with different initial relative densities $D_{\mathrm{r}}=90\%,$ $60\%,30\%$ are simulated. The stress path in different principal stress angles is shown in Figure 3, displaying that the direction of principal stresses is fixed and the stress path in the $[\frac{{\sigma \prime }_z-{\sigma \prime }_x}{2}-{\tau }_{\mathrm{zx}}]$ plane is a line.

Figure 3.

Stress paths in the $[\frac{{{\sigma }^{\prime }}_z-{{\sigma }^{\prime }}_x}{2}-{\tau }_{zx}]$ plane for monotonic loading tests in different directions of principal stress.

The model parameters a, g, m and n are the most important parameters and are calibrated by using the method described in section ‘Model calibration’. For the other 10 model parameters, refer to the work by Manzanal et al. ¹⁸ Some minor adjustments are made to achieve better simulation results. The model parameters are summarized in Table 1.

Table 1.

Model parameters of the proposed model for Toyoura sand.

M	${\alpha }_0$	$G_0$	$\nu$	$H_0$	${\beta }_0$	${\beta }_1$	m	n	$e_{\Gamma }$	${\lambda }_{\mathrm{c}}$	$\kappa$	a	g	$H_{\mathrm{U}0}$	${\gamma }_{\mathrm{u}}$	${\gamma }_d$
1.25	0.45	125	0.2	125	4.2	0.2	0.3	0.8	0.934	0.019	0.7	0.05	0.1	175	2.0	1.0

Figures 4–6 show the experimental results. Figures 7–9 show the simulated results by the proposed model. Both numerical simulation and experimental results display the following: (1) for the same initial relative density, as the principal stress angle becomes larger, the behaviour clearly becomes softer and more contractive; and (2) for the same principal stress angle, the contractive behaviour of sand decreases with an increase in the initial relative density.

Figure 4.

Measured shear properties on Toyoura sand under undrained conditions in different directions of principal stress with $D_{\mathrm{r}}=90\%$ . ³

Figure 5.

Measured shear properties on Toyoura sand under undrained conditions in different directions of principal stress with $D_{\mathrm{r}}=60\%$ . ³

Figure 6.

Measured shear properties on Toyoura sand under undrained conditions in different directions of principal stress with $D_{\mathrm{r}}=30\%$ . ³

Figure 7.

Predicted shear properties on Toyoura sand under undrained conditions in different directions of principal stress with $D_{\mathrm{r}}=90\%$ .

Figure 8.

Predicted shear properties on Toyoura sand under undrained conditions in different directions of principal stress with $D_{\mathrm{r}}=60\%$ .

Figure 9.

Predicted shear properties on Toyoura sand under undrained conditions in different directions of principal stress with $D_{\mathrm{r}}=30\%$ .

To illustrate the effects of principal stress orientation and state parameter by the proposed model, the experiments of Nakata et al. ³ are simulated with $a=g=m=n=0$ and $a=g=0$ .

• Case 1: $a=g=m=n=0$

In this case, the proposed model is reduced to the original PZ3 model; therefore, both of the effects of principal stress orientation and the effects of state parameters cannot be reflected. The other parameters are identical to those in Table 1. The results shown in Figure 10 indicate that at the different principal stress orientations of $\psi ={15}^{\mathrm{o}},{30}^{\mathrm{o}},{45}^{\mathrm{o}},{60}^{\mathrm{o}},{75}^{\mathrm{o}}$ and different initial relative ratios of $D_{\mathrm{r}}=30\%,60\%,90\%$ , all of the stress paths and the stress–strain behaviour curves are identical.

• Case 2: $a=g=0$

In this case, the proposed model is reduced to a state-parameter model considering the dilatancy of sand, which is similar to the model of Manzanal et al. ¹⁸ However, the effects of principal stress orientation are not reflected. The other parameters are identical to those in Table 1. The results shown in Figure 11 indicate that (1) for the same relative density, at the different principal stress orientations of $\psi ={15}^{\mathrm{o}},{30}^{\mathrm{o}},{45}^{\mathrm{o}},{60}^{\mathrm{o}},{75}^{\mathrm{o}}$ , the stress path and the stress–strain behaviour curves are identical; and (2) the effects of relative density can be reflected through the state parameter.

Figure 10.

Simulated results when $a=g=m=n=0$ .

Figure 11.

Simulated results when $a=g=0$ .

The prediction of the proposed model agrees reasonably well with experimental result by Nakata. The comparison confirms the validity of the proposed model in undrained monotonic loading experiments in different fixed directions of principal stress and state parameters.

Undrained cyclic rotational shear

Two sets of cyclic rotational shear experiments have been simulated to display the effects of principal stress orientation by the proposed model. The corresponding stress path is a circle in the $[\frac{{\sigma \prime }_z-{\sigma \prime }_x}{2}-{\tau }_{\mathrm{zx}}]$ plane (see Figure 12) and a point in the principal stress plane, which indicates that the total principal stresses remain constant but the principal stress orientation continues rotating.

Figure 12.

Stress paths in the $[\frac{{\sigma \prime }_z-{\sigma \prime }_x}{2}-{\tau }_{\mathrm{zx}}]$ plane for the undrained cyclic rotational shear tests.

Yang et al. ⁵ conducted a series of undrained cyclic rotational shear experiments on Toyoura sand. One of his experiments, Series I, was simulated and investigated first. The proposed model parameters for Toyoura sand are summarized in Table 1 and had an initial void ratio $e_0=0.707$ , $q=34.65\mathrm{kPa}$ . For comparison, the line elastic model and the PZ3 model are also employed to simulate the same stress path. The PZ3 model parameters are identical to the proposed model without the parameters of the principal stress orientation and state parameter.

The simulated hysteretic curve (Figure 13(b)) from the proposed model agrees well with the experimental curve (Figure 13(a)). The results of the linear elastic model and the PZ3 model are shown in Figure 13(c) and (d), respectively. It is seen that when we use the linear elastic model or the PZ3 model to simulate the undrained cyclic rotational shear experiment, the stress–strain behaviour is a straight line rather than a hysteretic curve, and the magnitude of the strain is negligible.

Figure 13.

Shear strain ${\epsilon }_{z\theta }$ and shear stress ${\tau }_{z\theta }$ relationship in the undrained cyclic rotational shear tests on Toyoura sand: (a) measured test result, ⁵ (b) predicted result using the proposed model, (c) predicted result by the linear elastic model and (d) predicted result by the PZ3 model.

The effects of the principal stress orientation may not have been reflected clearly because the experiment of Yang has only 35 cycles. Therefore, the simulations have been performed on Luan et al.’s ¹⁹ experiments which are similar to Yang’s experiments but have more cycles (approximately 200 cycles) and on Fujian standard sand, which is often used as the standard experimental sand in China and has a mean diameter $d_{50}=0.34\mathrm{mm}$ , uniformity coefficient $C_{\mathrm{u}}=1.542$ , specific gravity $G_{\mathrm{s}}=2.634$ , maximum void ratio $e_{max}=0.848$ and minimum void ratio $e_{min}=0.519$ .

The selection of model parameters for Fujian standard sand refers to the work by Liu and Song, and the model parameters are summarized in Table 2 (initial void ratio $e_0=0.76$ ).

Table 2.

Model parameters of the proposed model for Fujian standard sand.

M	${\alpha }_0$	$G_0$	$\nu$	$H_0$	${\beta }_0$	${\beta }_1$	m	n	$e_{\Gamma }$	${\lambda }_{\mathrm{c}}$	$\kappa$	a	g	$H_{\mathrm{U}0}$	${\gamma }_{\mathrm{u}}$	${\gamma }_d$
1.23	0.3	134	0.25	134	4.2	0.2	1.2	2.5	0.802	0.028	0.68	0.2	0.3	28	2.0	1.0

The simulation results with the proposed model shown in Figure 14 indicate that the higher the number of cycles, the higher the plastic strain. The magnitude of the simulated plastic strain is consistent with the experimental results shown in Figure 15, though there is some difference in the shape of hysteresis loops.

Figure 14.

Predicted shear strain ${\epsilon }_{z\theta }$ and shear stress ${\tau }_{z\theta }$ relationship using the proposed model in undrained cyclic rotational shear tests on Fujian standard sand.

Figure 15.

Measured shear strain ${\epsilon }_{z\theta }$ and shear stress ${\tau }_{z\theta }$ relationship in undrained cyclic rotational shear tests on Fujian standard sand. ¹⁹

Discussion and conclusion

Using a hierarchical approach, a new constitutive model is developed in this article by introducing principal stress angle and state parameter into the PZ3 model. The feature of stress–strain behaviour due to the change of principal stress orientation can be reasonably reflected when compared with the experiment. Meanwhile, the inclusion of state parameters makes it feasible to use the same set of parameters to predict the results of undrained monotonic loading tests in different initial confining pressures or initial void ratios. The new model uses void ratio e and critical stress ratio M as independent variables, which allows for more explicit physical meanings in the determination of the model parameters. Even though this new feature of the model was validated on undrained monotonic loading tests, which do not involve the change of void ratio during the test, the inclusion of the state parameter would allow the model to successfully simulate drained triaxial tests that involve a continuous change in the void ratio. Furthermore, a systematic model parameter calibration procedure is proposed, which enhances the practicality of the new model.

The capability of the model simulation is verified with the experimental results obtained from two different types of experiment involving both monotonic and cyclic loading. In accordance with experimental results, the proposed model can reasonably reflect the effects of principal stress orientation and the state-dependent dilatancy, especially in the stress condition of the pure principal stress rotation.

The current model only allows rotation of principal stress direction in a 2D plane. Nevertheless, this would already allow the use of the model to investigate problems involving saturated soil, pore fluid interactions and wave-induced dynamic responses of level seabeds with or without structures. However, in order to investigate 3D problems, the model should be extended by introducing two more principal stress angles or using the Cartesian stress components as base variables. However, it would be difficult to perform experiments with two different rotations of principal axes for model validation. One possible future direction is to split the one direction rotation of principal axis of hollow cylinder testing into two by rotating the coordinate system, which can be used to test future models involving more than one rotation of principal axes. Finally, because natural soil is layered due to the deposition, principal stress rotation in anisotropic soil should also be considered in the future.