Parameter Evaluation of Exponential-Form Critical State Line of a State-Dependent Sand Constitutive Model

Featured Application

In this paper, the influence of two model parameters$\mathbf{a}$and$\mathbf{b}$of critical state line on sand dilatancy behaviour are discussed. However, the emphasis is still put on liquefaction behaviour of clean loose sand. Moreover, a large number of real field studies have indicated that the presence of silt is very common at the liquefied sites during earthquakes. Therefore, this study could also be extended to investigate the potential influence of fine fraction of silty sand mixture on the critical state line and its instability behaviour.

Abstract

In sand constitutive models, it is of cardinal importance to consider a state parameter to distinguish the real dilatancy for cohesionless soils (sand), which is different from cohesive soils (clay). Thus, one of the key issues in simulating the sand behaviour is the better representation and parameter calibration of critical state line (CSL) for estimating contraction in loose state and dilatancy in dense state, respectively. For this purpose, a new exponential form for CSL with two model constants a and $\mathrm{b}$ has been presented in the literature. This paper provides a valuable insight into the two model constants, controlling the shape of the critical state line by simulating a uniform quartz reference sand (Hostun RF) in loose and dense states under undrained triaxial conditions. It can be concluded that the liquefaction behaviour in loose state is fundamentally affected by even a minor variation in model constant $\mathrm{a}$, but insensitive to model constant $\mathrm{b}$. Moreover, the linear fitting calibration of CSL recommended in the literature is complicated in consideration of the non-unified critical state line. Thus, the maximum void ratio in the natural state could be considered as a comparison basis on which to evaluate the liquefaction potential as an alternative. The numerical results showed good agreement with real experimental data. However, in dense state, the dilatant behaviour of sand was found to be mainly controlled by model parameter $\mathrm{b}$. In addition, the influence of a non-unified critical state under various confining pressures on the determination of $\mathrm{b}$ should not be neglected. With the correction of $\mathrm{b}$, the numerical results were found to be consistent with the experimental data concerning Hostun RF sand.

1. Introduction

Unlike metal materials, the shear strength of natural silica sand is controlled by internal frictions. Subjected to a relatively low loading, the destruction of sand specimens is caused by relative sand particle movements, instead of particle crushing. Especially during an earthquake, much of the severe damage of geotechnical structures can be attributed to sand liquefaction, which results in large land deformations due to positive excess pore water pressure. Therefore, a better constitutive modelling of sand on various states should consider both dilatant and liquefaction behaviours. In the triaxial stress–strain space, the dilatancy $\mathrm{d}$ is commonly described as the incremental ratio of the plastic volumetric strain to that of the deviatoric strain: $\mathrm{d}={\dot{\mathsf{\epsilon }}}_{\mathrm{v}}^{\mathrm{p}}/{\dot{\mathsf{\epsilon }}}_{\mathrm{q}}^{\mathrm{p}}$, where ${\mathsf{\epsilon }}_{\mathrm{v}}={\mathsf{\epsilon }}_1+2{\mathsf{\epsilon }}_3$, ${\mathsf{\epsilon }}_{\mathrm{q}}=\frac{2}{3}({\mathsf{\epsilon }}_1-{\mathsf{\epsilon }}_3)$, superscript ‘p’ stands for plastic and overdot stands for increment. On the basis of the experimental results of clayey soils, with the assumption of incremental plastic energy in the Cam-Clay model, Roscoe and Schofield [1] proposed $\mathrm{d}=\mathrm{M}-\mathrm{q}/{\mathrm{p}}^{\prime }=\mathrm{M}-\mathsf{\eta }\mathrm{d}=\mathrm{M}-\mathsf{\eta }$, where $\mathsf{\eta }$ represents stress ratio. Later, Roscoe and Burland [2] modified the dilatancy law in order to ensure that the direction of strain flux is horizontal while $\mathsf{\eta }=0$. However, Verdugo and Ishihara [3], and Li and Dafalias [4] pointed out that the above equations for cohesive soils does not adapt to cohesionless soils because $\mathrm{d}$ is expressed as a unique function of stress ratio with a constant material parameter $\mathrm{d}=\mathrm{d}\left(\mathsf{\eta },\mathrm{M}\right)$ in both the above versions. A state parameter containing both the effect of void ratio and confining pressure should be carefully taken into account in the dilatancy law $\mathrm{d}=\mathrm{d}\left(\mathsf{\eta },\mathsf{\psi },\mathrm{M}\right)$ for sand. The simplification of the former leads to multiple sets of model parameters because the same sands in loose and dense states are treated as totally different materials. Numerous experimental results show that the stress–strain response of sand is largely dependent on the relative density and effective confining pressure. For this purpose, a dimensionless scalar $\mathsf{\psi }={\mathrm{e}}_{\mathrm{act}}-{\mathrm{e}}_{\mathrm{c}}$ was introduced by Been and Jefferies [5], which estimates the distance between the current and critical void ratios under the same effective confining pressure, ${\mathrm{p}}^{\prime }$. Thus, the sand contracts in loose state while $\mathsf{\psi }>0$, and the sand dilates in dense state while $\mathsf{\psi }<0$. To avoid the counting problem of positive and negative signs, some researchers propose using the ratio of the current and critical void ratios $\mathsf{\psi }={\mathrm{e}}_{\mathrm{act}}/{\mathrm{e}}_{\mathrm{c}}$ instead of the difference measured in $\mathrm{e}-\log {\mathrm{p}}^{\prime }$ plan.

One of the fundamental issues in sand constitutive modelling is the precise representation of critical state line in terms of the effective confining pressure ${\mathrm{e}}_{\mathrm{c}}={\mathrm{e}}_{\mathrm{c}}\left({\mathrm{p}}^{\prime }\right)$. In general, plotting these lines of cohesionless sand in $\mathrm{e}-\log {\mathrm{p}}^{\prime }$ plan makes the expression algebraically more difficult. Firstly, the real experimental data [6,7] revealed that the critical state line from very loose to loose states falls within a spindle instead of a simple line. The simplification to a straight line leads to inaccurate calculations of sand instability behaviour. Secondly, the unification of the critical state of the same sand specimen under different confining pressures is still under discussion. Attempts have been made—Li et al. [8] found that these lines are almost straight in the $\mathrm{e}-{\left({\mathrm{p}}^{\prime }/{\mathrm{p}}_{\mathrm{a}}\right)}^{\mathsf{\alpha }}$ plan, instead of the conventional $\mathrm{e}-\log {\mathrm{p}}^{\prime }$ plan. More recently, Luo and Zhang [9] proposed a new exponential form of the critical state by using the following equation with a positive and a negative constant $\mathrm{a},\mathrm{b}$, respectively:

$$\begin{array}{c}{\mathrm{e}}_{\mathrm{c}}=\mathrm{a}·{\mathrm{e}}^{\mathrm{b}·{\mathrm{p}}^{\prime }/{\mathrm{p}}_{\mathrm{a}}}\end{array}$$

(1)

The comparison presented in the literature between numerical simulations and experimental results in triaxial space indicates that the above equation is adequate for describing the dilatancy in the high density and effective confining pressure range. Nevertheless, the influence of the two mentioned parameters on liquefaction behaviour is still unclear. Therefore, the emphasis of this paper is to evaluate the influence of model constants on exponential-form critical line. In addition, special attention is paid to investigating the effect of the non-unified critical state on the determination of model parameters. The description of the model [9] used in this paper is briefly given in the next section, followed by a detailed synthesis of the effect of the models parameters on critical state.

2. Model Description

In this section, all formulas are described in triaxial conditions, where ${\mathrm{p}}^{\prime }=\left({\mathsf{\sigma }}_1^{\prime }+2{\mathsf{\sigma }}_3^{\prime }\right)/3$ and $\mathrm{q}={\mathsf{\sigma }}_1^{\prime }-{\mathsf{\sigma }}_3^{\prime }$. Generally speaking, sand deformation could be caused by grain sliding and sand crushing. However, within the scope of the present work, it is supposed that the rupture of sand specimens is solely caused by the relative displacement between soil grains. Under this assumption, the yield criterion surface appears as an open wedge with ${\mathrm{p}}^{\prime }$-axis in the

${\mathrm{p}}^{\prime }-\mathrm{q}$ space with its apex at the origin [4,9]:

$$\begin{array}{c}\mathrm{f}=\mathrm{q}-\mathsf{\eta }{\mathrm{p}}^{\prime }=\frac{\mathrm{q}}{{\mathrm{p}}^{\prime }}-\mathsf{\eta }=0\end{array}$$

(2)

In the framework of elasto-plasticity theory and the range of small deformations and rotations, the mechanical behaviour of sand could be divided into elastic and plastic parts, namely:

$$\begin{array}{c}\dot{{\mathsf{\epsilon }}_{\mathrm{q}}}={\dot{\mathsf{\epsilon }}}_{\mathrm{q}}^{\mathrm{e}}+{\dot{\mathsf{\epsilon }}}_{\mathrm{q}}^{\mathrm{p}}\end{array}$$

(3)

$$\begin{array}{c}\dot{{\mathsf{\epsilon }}_{\mathrm{v}}}={\dot{\mathsf{\epsilon }}}_{\mathrm{v}}^{\mathrm{e}}+{\dot{\mathsf{\epsilon }}}_{\mathrm{v}}^{\mathrm{p}}\end{array}$$

(4)

In the elastic part, the shear modulus $\mathrm{G}$ is expressed by the following empirical equation [4,10]:

$$\begin{array}{c}\mathrm{G}={\mathrm{G}}_0\frac{{\left(2.97 - \mathrm{e}\right)}^2}{\left(1+\mathrm{e}\right)}\sqrt{{{\mathrm{p}}^{\prime }\mathrm{p}}_{\mathrm{a}}}\end{array}$$

(5)

where $\mathrm{e}$ is the current void ratio, and ${\mathrm{p}}_{\mathrm{a}}$ is the atmospheric pressure for normalization ($\approx 100\mathrm{kPa}$). Note that the unit of the elastic shear modulus $\mathrm{G}$ is given by the unit of mean effective stress ${\mathrm{p}}^{\prime }$, which causes ${\mathrm{G}}_0$ to be a dimensionless model constant. On the basis of elasticity theory, the bulk modulus $\mathrm{k}$ is equal to

$$\begin{array}{c}\mathrm{k}=\mathrm{G}\frac{2\left(1+\mathsf{\nu }\right)}{3\left(1-2\mathsf{\nu }\right)}\end{array}$$

(6)

where $\mathsf{\nu }$ is the Poisson’s ratio, which is considered constant and independent of effective confining pressure and the density index.

By applying the loading index $\mathrm{L}$ [11], the plastic deviatoric strain increment can be expressed as

$$\begin{array}{c}{\dot{\mathsf{\epsilon }}}_{\mathrm{q}}^{\mathrm{p}}=\mathrm{L}·\frac{\partial \mathrm{f}}{\partial \mathrm{q}}=\mathrm{L}·1=\mathrm{L}=\frac{\frac{\partial \mathrm{f}}{\partial {\mathrm{p}}^{\prime }}·{\dot{\mathrm{p}}}^{\prime }+\frac{\partial \mathrm{f}}{\partial \mathrm{q}}·\dot{\mathrm{q}}}{{\mathrm{K}}_{\mathrm{p}}}=\frac{\dot{\mathrm{q}}-\mathsf{\eta }·{\dot{\mathrm{p}}}^{\prime }}{{\mathrm{K}}_{\mathrm{p}}}\\ =\frac{{\mathrm{p}}^{\prime }·\left(-\frac{\mathrm{q}}{{\left({\mathrm{p}}^{\prime }\right)}^2}·{\dot{\mathrm{p}}}^{\prime }+\frac{1}{{\mathrm{p}}^{\prime }}·\dot{\mathrm{q}}\right)}{{\mathrm{K}}_{\mathrm{p}}}=\frac{{\mathrm{p}}^{\prime }\dot{\mathsf{\eta }}}{{\mathrm{K}}_{\mathrm{p}}}\end{array}$$

(7)

where $\dot{\mathsf{\eta }}=\frac{\partial \mathsf{\eta }}{\partial {\mathrm{p}}^{\prime }}·{\dot{\mathrm{p}}}^{\prime }+\frac{\partial \mathsf{\eta }}{\partial \mathrm{q}}·\dot{\mathrm{q}}=-\frac{\mathrm{q}}{{\left({\mathrm{p}}^{\prime }\right)}^2}·{\dot{\mathrm{p}}}^{\prime }+\frac{1}{{\mathrm{p}}^{\prime }}·\dot{\mathrm{q}}$ and ${\mathrm{K}}_{\mathrm{p}}$ is the plastic hardening modulus. As mentioned earlier, the key issue of describing the deformation of sand is the implantation of the state parameter into the dilatancy law $\mathrm{d}={\dot{\mathsf{\epsilon }}}_{\mathrm{v}}^{\mathrm{p}}/{\dot{\mathsf{\epsilon }}}_{\mathrm{q}}^{\mathrm{p}}$ [9]:

$$\begin{array}{c}d={\dot{\mathsf{\epsilon }}}_{\mathrm{v}}^{\mathrm{p}}/{\dot{\mathsf{\epsilon }}}_{\mathrm{q}}^{\mathrm{p}}={\mathsf{\psi }}^{\mathrm{n}}-\frac{\mathsf{\eta }}{\mathrm{M}}\end{array}$$

(8)

where $\mathrm{n}$ is a model dimensionless constant, and $\mathrm{M}$ is the stress ratio at critical state. Then the deviatoric and volumetric strain increment can be composed as follows:

$$\begin{array}{c}\dot{{\mathsf{\epsilon }}_{\mathrm{q}}}={\dot{\mathsf{\epsilon }}}_{\mathrm{q}}^{\mathrm{e}}+{\dot{\mathsf{\epsilon }}}_{\mathrm{q}}^{\mathrm{p}}=\frac{\mathrm{q}\dot{}}{3\mathrm{G}}+\frac{{\mathrm{p}}^{\prime }·\dot{\mathsf{\eta }}}{{\mathrm{K}}_{\mathrm{p}}}\end{array}$$

(9)

$$\begin{array}{c}\dot{{\mathsf{\epsilon }}_{\mathrm{v}}}={\dot{\mathsf{\epsilon }}}_{\mathrm{v}}^{\mathrm{e}}+{\dot{\mathsf{\epsilon }}}_{\mathrm{v}}^{\mathrm{p}}=\frac{\dot{{\mathrm{p}}^{\prime }}}{\mathrm{K}}+\mathrm{d}\frac{{\mathrm{p}}^{\prime }·\dot{\mathsf{\eta }}}{{\mathrm{K}}_{\mathrm{p}}}\end{array}$$

(10)

Here, it is of great importance to recall that the flow rule used in this study is non-associated, as the flux of plastic strain is directly defined by the dilatancy equation. Meanwhile, the flux of plastic strain is controlled by model constant parameter n. State parameter is implanted into plastic hardening modulus ${\mathrm{K}}_{\mathrm{p}}$ [9,12]:

$$\begin{array}{c}{\mathrm{K}}_{\mathrm{p}}=\mathrm{C}\mathrm{G}\left(\frac{\mathrm{M}}{\mathsf{\eta }}-\mathsf{\psi }\right)\end{array}$$

(11)

where $\mathrm{G}$ is the elastic shear modulus at the current state. Note that the initial state parameter values ${\mathsf{\psi }}_{\mathrm{ini}}={\mathrm{e}}_0/{\mathrm{e}}_{\mathrm{c}}$, under different density index and confining pressures are possible to them to be the same (see Figure 1). Thus, a dimensionless model parameter $\mathrm{C}\left({\mathrm{e}}_0\right)$ standing for the effect of initial void ratio is needed in Equation (11) in order to compensate this error, and it reads:

$$\begin{array}{c}\mathrm{C}=\mathsf{\alpha }{\mathrm{e}}_0^{\mathsf{\beta }}\end{array}$$

(12)

where $\mathsf{\alpha },\mathsf{\beta }$ are model dimensionless parameters and ${\mathrm{e}}_0$ is the initial void ratio.

It is worth noting that the plastic hardening modulus ${\mathrm{K}}_{\mathrm{p}}$ has the following special values:

• ${\mathrm{K}}_{\mathrm{p}}\left(\mathsf{\eta }=0\right)=\underset{\mathsf{\eta }=0}{\lim }\left(\mathrm{CG}\left(\frac{\mathrm{M}}{\mathsf{\eta }}-\mathsf{\psi }\right)\right)=\infty$, at the beginning of the test;
• ${\mathrm{K}}_{\mathrm{p}}\left(\mathsf{\eta }=\mathrm{M},\mathsf{\psi }=1\right)=\underset{\mathsf{\eta }=\mathrm{M};\mathsf{\psi }=1}{\lim }\left(\mathrm{CG}\left(\frac{\mathrm{M}}{\mathsf{\eta }}-\mathsf{\psi }\right)\right)$, in achievement of the critical state.

From Equations (11) and (12), it can be seen that the hardening or softening behaviour is determined by the sign of $\left(\frac{\mathrm{M}}{\mathsf{\eta }}-\mathsf{\psi }\right)$. Sand is initially in a dense state when $\mathsf{\psi }<1$, which results in $\frac{\mathrm{M}}{\mathsf{\eta }}-\mathsf{\psi }>0$ and sand response is hardening. Sand is initially in a loose state when $\mathsf{\psi }<1$. Its response is hardening with $\left(\frac{\mathrm{M}}{\mathsf{\eta }}-\mathsf{\psi }\right)>0$ and softening with $\left(\frac{\mathrm{M}}{\mathsf{\eta }}-\mathsf{\psi }\right)<0$. Thus, the introduction of state parameter in the plastic modulus succeeds in distinguishing the hardening and softening response of sand. Finally, the relationship between the stress and strain increments can be established in a matrix form as:

$$\begin{array}{c}\left(\begin{array}{c}{\dot{\mathsf{\epsilon }}}_{\mathrm{q}}\\ {\dot{\mathsf{\epsilon }}}_{\mathrm{v}}\end{array}\right)=\left[\mathrm{S}\right]\times \left(\begin{array}{c}\dot{\mathrm{q}}\\ {\dot{\mathrm{p}}}^{\prime }\end{array}\right)=\left[\begin{array}{cc}\frac{1}{3\mathrm{G}}+\frac{1}{{\mathrm{K}}_{\mathrm{p}}}& -\frac{\mathsf{\eta }}{{\mathrm{K}}_{\mathrm{p}}}\\ \frac{\mathrm{d}}{{\mathrm{K}}_{\mathrm{p}}}& \frac{1}{\mathrm{K}}-\frac{\mathrm{d}\mathsf{\eta }}{{\mathrm{K}}_{\mathrm{p}}}\end{array}\right]\times \left(\begin{array}{c}\dot{\mathrm{q}}\\ {\dot{\mathrm{p}}}^{\prime }\end{array}\right)\end{array}$$

(13)

$$\begin{array}{c}\left(\begin{array}{c}\dot{\mathrm{q}}\\ {\dot{\mathrm{p}}}^{\prime }\end{array}\right)={\left[\mathrm{S}\right]}^{-1}\times \left(\begin{array}{c}{\dot{\mathsf{\epsilon }}}_{\mathrm{q}}\\ {\dot{\mathsf{\epsilon }}}_{\mathrm{v}}\end{array}\right)=\left[\begin{array}{cc}3\mathrm{G}+9{\mathrm{AG}}^2& -3\mathrm{AKG}\mathsf{\eta }\\ 3\mathrm{AKGd}& \mathrm{K}-{\mathrm{AK}}^2\mathsf{\eta }\mathrm{d}\end{array}\right]\times \left(\begin{array}{c}{\dot{\mathsf{\epsilon }}}_{\mathrm{q}}\\ {\dot{\mathsf{\epsilon }}}_{\mathrm{v}}\end{array}\right)\end{array}$$

(14)

where $\mathrm{A}$ is equal to $1/\left(\mathrm{K}\mathsf{\eta }\mathrm{d}-3\mathrm{G}-{\mathrm{K}}_{\mathrm{p}}\right)$.

3. Effect of Model Constant on Critical State

3.1. Liquefaction Behaviour in Loose State

The recommended identification method of model constants $\mathrm{a}$ and $\mathrm{b}$ is shown in Figure 2 [9,13]. Several critical state points need to be plotted and fitted in the $\mathrm{e}-{\mathrm{p}}^{\prime }/{\mathrm{p}}_{\mathrm{a}}$ space to a straight line where $\mathrm{a}$ is the interception value and $\mathrm{b}\left(\mathrm{b}<0\right)$ is the slope of this line.

The critical state of Hostun RF sand from very loose (${\mathrm{I}}_{\mathrm{Di}}=-0.10$) to loose (${\mathrm{I}}_{\mathrm{Di}}=0.20$) state is shown in Figure 3. It can be observed that the unification of final state could not be seriously achieved. The simplification of the linear fitting is only a rough hypothesis and the determination of the model parameter $\mathrm{a}$ is in general sophisticated.

Thus, a parametric study of $\mathrm{a}$ was carried out by simulating a specimen of density indices equal to 0.00 with different values of $\mathrm{a}=1.05,1.00,0.95$ under an undrained triaxial condition. The material constants according to their functions in separate columns are given in

Table 1. Model parameters and physical properties of Hostun RF sand.

Elastic Parameters	Critical State Parameters	Dilatancy Parameters	Hardening Parameters	Sand Property
${\mathrm{G}}_0$ = 176	M = 1.22	n = 0.42	α = 0.63	${\mathrm{e}}_{\max }$ = 1
ν = 0.176	${\mathrm{b}}_0$ = −0.030 a = 1.00	-	β = −2.5	${\mathrm{e}}_{\min }$ = 0.656

. One should also note that all parameters listed in Table 1 were dimensionless.

In Figure 4, it is clear that the drop of deviator stress $\mathrm{q}$ and the static liquefaction phenomenon became increasingly important with the decrease in $\mathrm{a}$. Dilatant behaviour was observed for $\mathrm{a}=1.05$, whereas very contractive behaviour was observed for $\mathrm{a}=0.95,1.00$, which meant that the result in loose state depended, to a large extent, on the value of $\mathrm{a}$. The critical state line in terms of effective confining pressure is plotted in Figure 5 and the initial state of the simulated case is represented by a black rhombus point. For $\mathrm{a}=1.05$, the initial state was clearly below the critical state, which caused state parameter $\mathsf{\psi }$ to be greater than 1 during all calculations. However, for the case of $\mathrm{a}=0.95, 1.00$, the initial states were both over the critical state, resulting in the state parameter $\mathsf{\psi }$ being lower than 1.

Therefore, a very contractive behaviour could be expected, which is consistent with the generation of excess pore water pressure, as shown in Figure 6. All the above findings suggest that the value of $\mathrm{a}$ is the fundamental point to successfully simulate the liquefaction behaviour of very loose sand.

According to Figure 5, the form of critical state defined by Equation (1) can be separated into two parts. Firstly, it was almost constant before 100 kPa. Then, this line decreased significantly with the increase of ${\mathrm{p}}^{\prime }$. The interception value can be obtained by the following equation:

$$\begin{array}{c}e\left({\mathrm{p}}^{\prime }=0\right)=\underset{{\mathrm{p}}^{\prime }\to 0}{\lim }{\mathrm{ae}}^{\mathrm{b}\left({\mathrm{p}}^{\prime }/{\mathrm{p}}_{\mathrm{a}}\right)}=\mathrm{a}\underset{{\mathrm{p}}^{\prime }\to 0}{\lim }{\mathrm{e}}^{\mathrm{b}\left(0/{\mathrm{p}}_{\mathrm{a}}\right)}=a\end{array}$$

(15)

Within the range of 100 kPa, ${\mathrm{e}}_{\mathrm{c}}$ increased slowly with the decrease in ${\mathrm{p}}^{\prime }$ and achieved its ultima value ${\mathrm{e}}_{\mathrm{c}}=\mathrm{a}$, while ${\mathrm{p}}^{\prime }=0$. Equation (15) gives the physical sense of the model constant $\mathrm{a}$, which represents the possible maximum void ratio of sand without any confining pressure. In a study associated with liquefaction, it is also reasonable to consider the maximum void ratio as a comparison basis for evaluating the contractive behaviour. Thus, the substitution of model constant $\mathrm{a}$ by a sand property ${\mathrm{e}}_{\max }$ is recommended by the authors.

The variation of critical state in terms of model constant $\mathrm{b}$ is plotted in Figure 7. In the range of effective confining pressure ${\mathrm{p}}^{\prime }$ greater than 100 kPa, the form of critical state was mainly controlled by model parameter $\mathrm{b}$ and the critical state line decreased sharply by reducing the value of $\mathrm{b}$. However, it seemed that the variation in $\mathrm{b}$ did not remarkably affect the state parameter for very low density index. Thus, the influence of model constant $\mathrm{b}$ on liquefaction behaviour was then studied. By altering $\mathrm{b}$ from −0.035 to −0.025, the effective stress paths were almost coincident (Figure 8), except for minor differences at high-strain level.

For ascertaining the predominant effect of the model constant $\mathrm{a}$, numerical simulations under confining pressures of 200 and 400 kPa were performed by considering $\mathrm{a}={\mathrm{e}}_{\max }=1.00$. It can be seen that there was a good agreement (Figure 9) between the experimental and numerical results. The model performance on liquefaction behaviour in very loose state was re-confirmed.

3.2. Dilatant Behaviour in Dense State

As mentioned earlier, the unification of critical state under various confining pressures is still under discussion. According to Figure 10, it appears that the critical mean effective stress ${\mathrm{p}}^{\prime }$ increased with the increase in initial confining pressure for all density indexes.

It was analysed in the last section that the performance of numerical results in loose and very loose states were insensitive to variation in model parameter $\mathrm{b}$ due to the fact that their initial density indexes were almost located in the first part of proposed critical state line, which was mainly determined by model constant $\mathrm{a}$. However, for dense specimens whose initial states fell totally within the second part of the proposed critical line, the determination of model constant $\mathrm{b}$ played a very important role instead of the model constant $\mathrm{a}$. A case of density index of 0.70 under confining pressure equal to 200 kPa was investigated by keeping model parameter $\mathrm{b}$ constant. The stress–strain curve is shown in Figure 11, and it is clear that the result is not ideal because the state parameter was over-estimated.

Therefore, it seems better to suppose that model parameter $\mathrm{b}=\mathrm{b}\left({\mathrm{p}}^{\prime },{\mathrm{b}}_0\right)$ is a function of initial confining pressure instead of a model constant, where ${\mathrm{b}}_0$ is determined by original model parameter identification method [9,13]. For a well-studied sand, the variable model parameter $\mathrm{b}=\mathrm{b}\left({\mathrm{p}}^{\prime },{\mathrm{b}}_0\right)$ can be directly determined by fitting the critical state line for different confining pressures. For a sand less-investigated, the authors recommend the following equation:

$$\begin{array}{c}b=\{\begin{array}{c}{\mathrm{b}}_0,\mathrm{while} {\mathrm{p}}_0^{\prime }<100 \mathrm{kPa}\\ -0.12{\mathrm{p}}_0^{\prime }{}^{-10{*\mathrm{b}}_0}, \mathrm{while} {\mathrm{p}}_0^{\prime }>100 \mathrm{kPa}\end{array}\end{array}$$

(16)

Figure 12 shows that with the correction of model parameter $\mathrm{b}$, the model was capable of correctly simulating a case for dense sand specimens. Some small differences can be attributed to the observation that deformation of dense sand was not homogenous and that the expansion of the sample’s section was significant.

4. Disussion

The presented model works well for cohesionless soil because it is able to imitate the possible contractive and dilatant behaviour of sand by incorporating the state parameter into its framework. However, the effect of isotropic consolidation on soil response is neglected. Moreover, the yield criterion is proposed under a ‘strong’ assumption that the sand particle is incompressible during shearing. However, particle breakage is inevitable, especially near the shear band. Thus, the necessity of a ${\mathrm{p}}^{\prime }$-controlling cap into model yield criterion is under discussion. On the basis of real measurements of Hostun RF sand, the authors proposed a correction formula of model parameter $\mathrm{b}$ in order to capture the influence of confining pressure on the critical state of tested sands. However, more investigations are still needed to verify whether it could be suitable for other sands.

5. Conclusions

The present work presents a parameter evaluation of two model constants based on the exponential form of critical state line. Undrained triaxial tests of clean Hostun RF sand were simulated. Special attention was paid to studying the model performance on sand liquefaction and dilatant behaviours. The principal findings were as follows:

• Model constant $\mathrm{a}$ had a predominant effect in loose state because it had a clear physical meaning. The substitution of $\mathrm{a}$ by ${\mathrm{e}}_{\max }$ is recommended, and the numerical results were consistent with the real experimental data. Moreover, in loose state, model constant$\mathrm{b}$ did not play an important role.
• In dense state, the influence of initial confining pressure on critical state cannot be neglected and $\mathrm{b}$ should be corrected with the initial confining pressure. An empirical equation is proposed on the basis of Hostun RF clean sand.