The Constitutive Model of Rockfill Based on Property-Dependent Plastic Potential Theory for Geomaterials

Abstract

To better control the strength and deformation of the roadbed, a constitutive model of rockfill was established based on property-dependent plastic potential theory for geomaterials. The effect of the particle gradation on the anisotropy was described in the model. According to the effect of the particle grading and crushing on the fractal dimension, the fractal theory and fabric tensor were introduced to establish the yield and failure criteria of the rockfill. By combining the property-dependent concepts of the materials and the results of the rockfill strength test, a critical state line considering the microstructure, fractal dimension, particle breakage, and stress state of the rockfill was established. The dilatancy equation was derived based on the novel potential theory and the hardening criterion affected by the critical state was established. A constitutive model of the rockfill in the general stress space was established under the framework of the novel potential theory. The 3D strength and its intensity change in the π plane were simulated through the drainage strength test results, which verified the description of the critical state under various stress paths. By simulating the stress-strain relationship, the validity and rationality of the model were verified.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Rockfill is a common engineering material that is largely used in railway, transportation, and highway roadbeds construction projects due to its high filling density, strong permeability, good compacting performance, and high strength. Studying more about the stress, strain, and deformation characteristics of rockfill materials under various loads in detail and establishing a reasonable and effective constitutive model were of great significance to the theoretical analysis of rockfill engineering application and the numerical calculation and analysis of transportation geotechnics [1–3].

The constitutive relationship is not only the core of modern soil mechanics, but also an important entry point for studying the mechanical properties of rockfill materials. Compared with sand, the particle size and the void of rockfill material are larger. Under complex stress conditions, the spatial fabric and particle size distribution of rockfill are usually more prone to slipping and breaking. The complex macroscopic mechanical properties often show a tendency for nonlinear forms of strength and volume change, which affect the stability of a rockfill fabric [4]. This has brought some difficulties in studying the constitutive relationship of rockfill. Some researchers have done a large amount of work and achieved positive research results such as the strain-softening and dilatancy of rockfill [5], material constants and particle characteristics [6], load and creep [7, 8], particle breakage and the relationship between intermediate principal stress coefficient b value and particle breakage [9], loading path [10], wetting deformation [11], the coupling effect of mean effective stress p and deviatoric stress $q$ on deformation [12], and so on.

In the context of generalized plasticity theory, Wang et al. [13] and Liu et al. [14] established a constitutive model of rockfill to describe the particle breakage considering critical state and dilatancy. Xiao and Liu [15] proposed a critical state line for particle crushing called the breakage critical state line (BCSL) and established an elastoplastic model for rockfill materials considering the state dependence and particle crushing based on critical state soil mechanics. Liu et al. [16] and Fang et al. [17] proposed an elastoplastic constitutive model and a state-dependent 3D multimechanism boundary surface model, respectively, by introducing state-related parameters, but the models still had certain limitations in the adaptation and expansion of complex stress paths. Liu and Chen [18] established an exponential-parabolic nonlinear elastic constitutive model reflecting the porosity and density of rockfill material based on the improved Nanshui double yield surface model. Besides, Brito et al. [19] established a new model specifically describing soil-rockfill mixtures (SRM).

In summary, many researchers have established a constitutive model that could describe the hardness, dilatancy, loading, and unloading of rockfill materials suitably under complex stress paths based on elastic-plastic theory, subloading surface theory, generalized plasticity theory, boundary surface theory, and others. These constitutive models have achieved positive and effective results in describing the basic strength and deformation characteristics of materials, such as friction and dilatancy. However, the existing models still have the shortcomings of narrow application range, many model parameters, and complicated forms. Meanwhile, these models cannot accurately describe the isotropy and anisotropy of rockfill materials. Furthermore, some models also cannot reflect the plastic deformation that is caused by the rotation of the principal stress. Therefore, determining how to correctly describe the deformation characteristics of rockfill materials, quoting the small amount of necessary mechanical parameters, and unifying the anisotropy and isotropic characteristics of the material characteristics have become the research difficulties. This study tries to adopt the potential theory to solve the above problems. Compared with the traditional constitutive models, the mathematical principle is clear and is not based on a plastic postulate. It connected physical assumptions and mathematical foundations established by the constitutive model. Thus, it forms a complete set for a theoretical system. The main principle of potential theory involves taking the principal stress and principal strain of the main space or its increment as a mathematical vector and using the idea of vector fitting to fit the known vector [20, 21].

Under the framework of the theory of continuum mechanics, Li et al. [22] linked the strain distribution rule of geotechnical materials with its fabric properties and proposed the property-dependent plastic potential theory. The proposed theory was verified and applied to sand constitutive models with a good description result. Compared to the noncoaxial plastic theory proposed by Gutierrez et al., the theory of Li et al. has a better description of the noncoaxial characteristics. However, it is necessary to consider the differences in the fabrics and the mechanical properties between sand and rockfill, so it is not easy to directly apply the potential theory of sand to rockfill. The particle fractal and fabric of rockfill material can be described as the microscopic characteristics of the material, and the size effect builds a bridge between the sand research and the rockfill material. The key idea of the plasticity theory related to material properties in this research was that, according to the material’s meso-structural properties, it was assumed that the material properties could be described by a fabric tensor. Then, the strain distribution of the material could be affected by the material properties. Based on this, in this research, the property-dependent plastic potential theory will be extended to describe anisotropy characteristics of rockfill, and a constitutive model is proposed for rockfill with general stress space. Based on the fractal dimension, the novel state variables will be introduced. Additionally, in this research, the yield failure criterion of rockfill will be proposed, the dilatancy equation of rockfill will be derived by proposed potential theory, hardening criterion considering the influence of critical state will be established, and the rationality of its 3D strength description will be verified through large-scale triaxial test results and the determination of effectiveness.

2. Fractal Theory

Rockfill has the particular properties of coarse particles with large sizes and its easily broken characteristics. There is a significant relationship among the changing gradation along with the stress, dilatancy, and strength-deformation characteristics in rockfill materials [23, 24]. Therefore, it is important to study the mechanical properties of rockfill materials by quantitatively describing its gradation of rockfill materials which approaches the true properties and status of raw materials. Zhang and Zhang [25] found that the particle yield and compression state were closely related to the particle-size distribution and fractal curve characteristics of particle breakage through a one-dimensional compression test. Besides, the study found that the fractal dimension $D$ could be used to describe the rockfill particle-size distribution. This was especially true for the two indicators $C_u$ and $C_c$ used to characterize the particle-size distribution of the large particle rockfill material of the earth-rock dam and highway roadbed. The indicators $C_u$ and $C_c$ could accurately determine the size distribution and the optimal content of fine particles. This means that the fractal theory had a strong theoretical and practical significance for the study of constitutive relationships [26].

The fractal relationship for the particle size of the rockfill and the particle mass and particle volume was described by the fractal model. It could be expressed as$$\begin{array}{}\text{(1)}& \lg \frac{Mr<d_{\max }}{M_T}=3-D\lg \frac{d_{\max }}{{\lambda }_V},\end{array}$$where $M$ is the mass of the certain particles, M_T is the total mass of the particles, and d_max and ${\lambda }_V$ stand for the maximum particle size and the characteristic particle size, respectively. $D$ is the fractal dimension of the rockfill particle size, which was defined as follows:$$\begin{array}{}\text{(2)}& D=3-\frac{\ln Mr<d_i/M_T}{\ln d_i/d_{\max }}.\end{array}$$

This shows that as the nonuniformity coefficient of the rockfill gradation becomes larger, the number of particles $Nr$ that are less than r from equation (2) increases. Thus, the dimension $D$ was larger. Formula (2) better reflected the grading of rockfill, so it could be used to represent the grading of rockfill.

Studies have shown that the crushing of rockfill particles becomes increasingly severe with an increase in confining pressure [27]. The reason for this is that the particles inside the rockfill are easily broken with the multiple effects occurring under higher confining pressures. This usually caused the quality of particles with larger particle sizes to decrease. However, the relative density of the rockfill material had little effect on the particle crushing. The specific behavior was as follows. The finer the initial rockfill gradation was (i.e., the larger the initial fractal dimension was), the smaller the crushing degree was. The experiment confirmed that the fractal dimension $D$ had a good linear relationship with ${\sigma }_3/P_a$, which could be obtained by binary linear regression analysis:$$\begin{array}{}\text{(3)}& D=l+\kappa \frac{{\sigma }_3}{P_a}+\beta D_0,\end{array}$$where $l$, $\kappa$, and $\beta$ are model parameters that could be calculated from the experimental data.

3. Constitutive Model Construction

3.1. Elastic Description

The elastic shear modulus of rockfill that describes its pressure hardness characteristics based on the empirical equation (4) can be found in the work by Richart et al.:$$\begin{array}{}\text{(4)}& G=G_0\frac{{2.97-e}^2}{1+e}{p\cdot p_a}^{0.5},\end{array}$$where e is the current void ratio, P_a is the atmospheric pressure (P_a = 101 kPa), and $G_0$ is the material constant. Therefore, the bulk modulus of elasticity can be expressed as$$\begin{array}{}\text{(5)}& K=G_0\frac{21+\nu }{31-2\nu }\cdot p_a\frac{{2.97-e}^2}{1+e}{\frac{p}{p_a}}^{0.5},\end{array}$$where $v$ is Poisson’s ratio. Equations (4) and (5) indicate that the strength of a material increases with the enhancement of the stress level. This conforms to the pressure hardness characteristic for which the strength of the rockfill increases with the confining pressure.

3.2. Yield Surface

Based on the work of the failure criterion for anisotropic sand and rockfill by Li et al. [28, 29], the novel anisotropic state variable was introduced, and the yield criterion in this research could be written as$$\begin{array}{}\text{(6)}& f=q-1+\varsigma AMg{\theta }_{\sigma }p=0,\end{array}$$where $M$ is the internal variable of hardening and $\varsigma$ is the weight coefficient, which was usually equal to 0.2. $p$ and $q$ are the mean stress and the generalized deviatoric stress, respectively. $p={\sigma }_{kk}/3$, $q=\sqrt{3J_2}$, $J_2=s_{ij}{s}_{ij}/2$, and $s_{ij}={\sigma }_{ij}-{\delta }_{ij}p$. ${\delta }_{ij}$ is the Kronecker tensor (when $i\ne j$, ${\delta }_{ij}$ = 0 and when $i=j$, ${\delta }_{ij}$ = 1). ${\theta }_{\sigma }={\sin }^{-1}3\sqrt{3}{J}_3/2J_2^3/3$, $J_3=s_{ij}{s}_{jk}{s}_{ki}/3$, and ${\theta }_{\sigma }$ is the stress Lode angle. $A$ is the anisotropic state variable. $g{\theta }_{\sigma }$ is a function of ${\theta }_{\sigma }$ that reflects the influence of the intermediate principal stress. The elliptic function of $g{\theta }_{\sigma }$ in equation (6) was proposed by William and Warnke [30]:$$\begin{array}{}\text{(7)}& g{\theta }_{\sigma }=\frac{21-c^2\cos \pi /6+{\theta }_{\sigma }+2c-1\sqrt{41-c^2{\cos }^2\pi /6+{\theta }_{\sigma }+c5c-4}}{41-c^2{\cos }^2\pi /6+{\theta }_{\sigma }+{2c-1}^2},\end{array}$$where $c$ is the ratio of the triaxial tension strength to the triaxial compression strength. To ensure the convexity of the yield surface, the law of compression and elongation was used. The yield surface had to satisfy the following conditions: when ${\theta }_{\sigma }$ was equal to π/6 or −π/6, $\text{d}g{\theta }_{\sigma }/\text{d}{\theta }_{\sigma }$ was equal to 0. When $g-\pi /6$ was equal to $c$, $g\pi /6$ was equal to 1. $M$_f is the peak failure stress ratio, which stands for equation (8) with the peak internal friction angle $\phi$_f:$$\begin{array}{}\text{(8)}& M_f=\frac{6\sin {\phi }_f}{3-\sin {\phi }_f}.\end{array}$$

It was found that the strength envelope of the rockfill materials exhibited a nonlinear trend. With the increase of the confining pressure, the slope of the envelope gradually decreased and the longitudinal intercept increased. Therefore, the Mohr–Coulomb yield surface had some limitations in describing the strength characteristics of the rockfill. The test showed that the relationship between the ratio of the deviatoric stress and the atmospheric pressure $\tau /P_a$ and the ratio of normal stress and the atmospheric pressure $\sigma /P_a$ followed the power function in$$\begin{array}{}\text{(9)}& \frac{\tau }{p_a}=A_c{\frac{\sigma }{P_a}}^{B_c}.\end{array}$$

The peak internal friction angle decreased with the increase of the confining pressure, which could be expressed as$$\begin{array}{}\text{(10)}& {\phi }_f=\arctan A_c{\frac{\sigma }{P_a}}^{B_c-1},\end{array}$$where $A$_c and $B$_c are the strength parameters. The amplitude parameters of the x₁−x₂ and x₂−x₃ planes could be obtained similarly as follows:$$\begin{array}{}\text{(11)}& \begin{array}{c}{A}_c=f{D}_0,e,\\ B_c=fe.\end{array}\end{array}$$

$A$_c was proportional to the fractal dimension $D$₀ and inversely proportional to the void ratio e, while the fractal dimension $D$₀ had little effect on $B$_c.

3.3. Orthotropic Fabric Tensor

The anisotropic state variable based on the research of the quantitative description of mesoscopic fabric by Li et al. [31] in equation (12) was expressed as$$\begin{array}{}\text{(12)}& \mathbf{A}=\frac{{\overline{\sigma }}_{ij}}{{\sigma }_m}\overline{F}{}_{ij}-{\frac{{\overline{\sigma }}_{ij}}{{\sigma }_m}\overline{F}{}_{ij}}_0=\frac{{\overline{\sigma }}_{ij}}{{\sigma }_m}\overline{F}{}_{ij}-A_0,\end{array}$$where ${\overline{\sigma }}_{ij}$ is stress tensor. In particular, ${\overline{\sigma }}_{ij}$ is principal stress tensor and ${\sigma }_m$ is the partial tensor of the stress tensor when the fabric tensor is expressed in terms of the main fabric space. In equation (12), the equation ${\overline{\sigma }}_{ij}/{\sigma }_m{\overline{F}}_{ij}=\overline{A}$ represents the secondary anisotropic state variable induced by an external load. Thus, $A$₀ is the native anisotropic state variable. In conclusion, according to the work by Li et al. [31], the mesoscopic fabric tensor was given by$$\begin{array}{}\text{(13)}& {\overline{F}}_{ij}=\begin{array}{ccc}\frac{1+a_1+a_2+a_1{a}_2}{3+a_1+a_2-a_1{a}_2}& & \\ & \frac{1+a_1-a_2-a_1{a}_2}{3+a_1+a_2-a_1{a}_2}& \\ & & \frac{1-a_1+a_2-a_1{a}_2}{3+a_1+a_2-a_1{a}_2}\end{array},\end{array}$$where a₁ and a₂ are orthotropic amplitude parameters that could be expressed as$$\begin{array}{}\text{(14)}& \begin{array}{c}{a}_1=\frac{1}{2N}\sqrt{{{\displaystyle \sum _{k=1}^{2N}{\cos }^2{\theta }_1^k-{\sin }^2{\theta }_1^k{\cos }^2{\alpha }^k}}^2+{{\displaystyle \sum _{k=1}^{2N}\sin 2{\theta }_1^k\cos {\alpha }^k}}^2},\\ a_2=\frac{1}{2N}\sqrt{{{\displaystyle \sum _{k=1}^{2N}{\cos }^2{\theta }_1^k-{\sin }^2{\theta }_1^k{\sin }^2{\alpha }^k}}^2+{{\displaystyle \sum _{k=1}^{2N}\sin 2{\theta }_1^k\sin {\alpha }^k}}^2},\\ a_3=\frac{1}{2N}\sqrt{{{\displaystyle \sum _{k=1}^{2N}{\sin }^2{\theta }_1^k\cos 2{\alpha }^k}}^2+{{\displaystyle \sum _{k=1}^{2N}{\sin }^2{\theta }_1^k\sin 2{\alpha }^k}}^2},\end{array}\end{array}$$where a₁, a₂, and a₃ are the F_x1−F_x3, F_x1−F_x2, and F_x2−F_x3 orthotropic amplitude parameters on three orthogonal surfaces. According to the microscopic test, all three parameters could be determined. The range of the values allowed the amplitude parameters to be in the range of 0 to 1. It reflected the anisotropy degree of the materials on the three deposition surfaces. θ^k is the angle between the long axis orientation of the kth particle and the corresponding coordinate axis. N is the total number of measured rockfill particles and ${\alpha }^k$ is the kth angle between the projection of the long axis of the particle on the horizontal plane and the x_k axis.

The definition of the 3D fabric equation (13) required that F_ij be equal to 1. There were only two independent variations for the determined fabric amplitude parameters on the a₁, a₂, and a₃ surfaces. The expressions of the three orthogonal fabric tensors (F₁, F₂, and F₃) obtained by (14) were functions of a₁ and a₂, a₁ and a₃, and a₂ and a₃, respectively. Therefore, only the information for the microscopic arrangement of the rockfill on any two of the three orthogonal surfaces of the particles was obtained through experiments, so the orthotropic fabric tensors of the rockfill material could be determined completely.

4. Model Framework

4.1. Property-Dependent Plastic Potential Theory

To reflect the anisotropy of the rockfill fabric and the uniqueness of the rockfill critical state, the consideration of rockfill mesoscopic fabric material properties related to the plastic potential theory was put forward based on material status-related dilatancy theory. It was assumed that the material properties could be represented by the fabric tensor and that the rockfill material properties would influence the strain distribution. For an incremental property of the geotechnical materials, the general expression of the property-dependent plastic potential related to rockfill materials can be written as$$\begin{array}{}\text{(15)}& \text{d}{\epsilon }_{ij}^p=\text{d}\lambda \frac{\partial g}{\partial {\sigma }_{il}}{\overline{\mathbf{F}}}_{lj},\end{array}$$where $g$ is the plastic potential function, dλ is the increment of the plastic factor, and ${\overline{F}}_{ij}$ is the tensor of the fabric tensor F_ij after the transformation of equation (15). According to the transformation characteristics of equation tensors, equation (15) was further simplified to$$\begin{array}{}\text{(16)}& {\mathbf{F}}_{ij}=p_F{\delta }_{ij}+{\tilde{s}}_{ij}=p_F{\overline{F}}_{ij},\end{array}$$where $p_F$ and ${\tilde{s}}_{ij}$ are the spherical tensor and the partial tensor of the fabric tensor, respectively.

Equations (15) and (16) could be rearranged as another expression of the material property potential theory:$$\begin{array}{}\text{(17)}& \text{d}{\epsilon }_{ij}^p=\text{d}\lambda \frac{\partial g}{\partial {\sigma }_{ij}}+\text{d}\lambda \frac{\partial g}{\partial {\sigma }_{il}}{\tilde{s}}_{lj},\end{array}$$where $\text{d}\lambda \partial g/\partial {\sigma }_{ij}$ is traditional plastic potential theory. When the material was isotropic, $\text{d}\lambda \partial g/\partial {\sigma }_{ij}{\tilde{s}}_{lj}$ was equal to 0, and the plastic factor dλ was a nonnegative proportional coefficient in equation (17).

4.2. Hardening Law

The hardening law of internal variables was adopted in the incremental hyperbolic form proposed by Li and Dafalias [32]. It was given by $$\begin{array}{}\text{(18)}& \text{d}{\epsilon }_s^p=\frac{M}{h_{s}G{M}_p-M}\text{d}M,\end{array}$$where $G$ is the shear modulus. $M$_p is the peak stress ratio in a conventional triaxial compression test, and h_s is the material constant.

It must be noted that the peak stress ratio $M$_p had to not be a constant, but rather a quantity related to the state parameters ψ. In this research, the equation for rockfill was adopted, $M_p=M_{cs}\exp -k_p\psi$, as suggested by Li and Dafalias, and k_p was the constant of the model. The state parameter ψ was expressed as$$\begin{array}{}\text{(19)}& \psi =e-e_c,\end{array}$$where e_c is the critical void ratio which could be determined from the following equation:$$\begin{array}{}\text{(20)}& e_c=e_{\Gamma }-{\lambda }_c{\frac{p}{p_a}}^{\xi },\end{array}$$

The results showed that there was a good linear relationship between the critical void ratio intercept e_Г and the initial void ratio e₀, which satisfied the relationship$$\begin{array}{}\text{(21)}& e_{\Gamma }=c^{\prime }+b^{\prime }{e}_0.\end{array}$$

Moreover, there was a linear strip relationship between the critical void ratio intercept e_Г and the fractal dimension $D$, so the critical void ratio intercept e_Г was expressed as$$\begin{array}{}\text{(22)}& e_{\Gamma }=c^{\prime }-a^{\prime }D+b^{\prime }{e}_0,\end{array}$$where a′, b′, and c′ are the model parameters that could be calculated from the experimental data.

For the value of h_s, Li and Dafalias [32] considered that it was related to the void ratio, taken as h_s = h₁−h₂e, where h₁ and h₂ were the constants of the model. In addition, the model parameters were also affected by anisotropy, and its expression could be rearranged as$$\begin{array}{}\text{(23)}& h_s=h_1-h_{2}e1+k_{h}A,\end{array}$$where k_h is the model parameter and $A$ is the defined anisotropic state variable in equation (12), when the material is isotropic, $A\equiv 0$.

To consider the effect of the anisotropy, the expression of the peak stress ratio $M$_p could be changed into the function of an anisotropic state variable:$$\begin{array}{}\text{(24)}& M_p=M_{cs}g{\theta }_{\sigma }{e}^{-k_p\psi A},\end{array}$$where $\psi A$ is the state parameter related to the initial anisotropy described in equation (24) and k_p is the model parameter. The expression of the peak stress ratio showed that the peak stress ratio was not a constant but rather an obvious function of the anisotropic state variable.

4.3. Dilatancy Equation

The dilatancy equation could be deduced from property-dependent plastic potential theory. Based on the basic theory of rockfill mechanics in a critical state, when the rockfill reached the critical state, the volume change was equal to zero and the stress ratio reached the critical state pressure ratio ηc. Lü et al. found that stress dilatation is a function related to the state of geomaterials [33]. At this point, the deformation of the rockfill only had shear strain, so its energy relationship satisfied$$\begin{array}{}\text{(25)}& \text{d}{W}^p=p\text{d}{\epsilon }_v^p+q\text{d}{\gamma }^p\cos {\theta }_{\sigma }-{\theta }_{\text{d}\epsilon }={\eta }_{c}p\text{d}{\gamma }^p,\end{array}$$where η_c is the stress ratio under the critical state. According to the energy relationship of the critical state given in equation (25), the dilatancy equation based on the potential theory related to material properties could be obtained as follows:$$\begin{array}{}\text{(26)}& \frac{\text{d}{\epsilon }_v^p}{\text{d}{\gamma }^p}={\eta }_c-\cos {\theta }_{\sigma }-{\theta }_{\text{d}\epsilon }\frac{q}{p}.\end{array}$$

The dilatancy equation is usually determined by conventional triaxial tests, and it can be expressed as follows:$$\begin{array}{}\text{(27)}& d=d_0{\eta }_c-\cos {\theta }_{\sigma }-{\theta }_{\text{d}\epsilon }\frac{q}{p},\end{array}$$where d₀ is the dilatation coefficient, which was determined by the experiment, θ_σ is the stress Lode angle, and θ_dε is the strain increment Lode angle in principal stress space. When the stress ratio was greater than the critical state (d < 0) for the dilatancy equation, rockfill dilatancy occurred. When the stress ratio was less than the critical stress (d > 0), the rockfill contracted. When the stress ratio was equal to the critical state (d = 0), the rockfill was in the state of phase, as shown in Figure 1.

[figure omitted; refer to PDF]

In the principal stress space, the strain increment θ_dε of the Lode angle was defined as$$\begin{array}{}\text{(28)}& {\theta }_{\text{d}\epsilon }=\arctan \frac{1}{\sqrt{3}}\frac{2\text{d}{\epsilon }_2-\text{d}{\epsilon }_1-\text{d}{\epsilon }_3}{\text{d}{\epsilon }_1-\text{d}{\epsilon }_3},\end{array}$$where the main strain was determined by property-dependent plastic potential theory in equation (25). Assuming $q/p=\eta$, equation (29) could be rearranged as$$\begin{array}{}\text{(29)}& \text{d}q=\eta \text{d}p+p\text{d}\eta ,\end{array}$$

By substituting equation (29) into equation (27), we obtained the following:$$\begin{array}{}\text{(30)}& \frac{\text{d}p}{p}=-\frac{\text{d}\eta }{d+\eta },\end{array}$$

The plastic potential function was obtained by integrating equation (30):$$\begin{array}{}\text{(31)}& g=q-p\frac{d_0{\eta }_c}{d_0{c}_{\psi }-1}1-{\frac{p}{p_0}}^{d_0{c}_{\psi }-1}=0,\end{array}$$where $c_{\psi }$ is the state parameter, the value could be obtained from the equation:$c_{\psi }=\cos {\theta }_{\sigma }-{\theta }_{\text{d}\epsilon }$ under the critical state of rockfill.

Considering the influence of anisotropy on the critical state stress ratio η_c, the expression of η_c in 3D space was defined as follows:$$\begin{array}{}\text{(32)}& {\eta }_c=1+\varsigma A{M}_{cs}g{\theta }_{\sigma }{e}^{k_d\psi A},\end{array}$$where k_d is the model parameter.

4.4. Incremental Stress-Strain Relationship

Based on the expression method of the strain increment, the total strain could be expressed as$$\begin{array}{}\text{(33)}& \text{d}{\epsilon }_{ij}=\text{d}{\epsilon }_{ij}^e+\text{d}{\epsilon }_{ij}^p,\end{array}$$where $\text{d}{\epsilon }_{ij}$ is the total strain increment and $\text{d}{\epsilon }_{ij}^e$ is the elastic strain increment, the expression for which is shown in equation (34). According to the potential theory related to the material properties, $d{\epsilon }_{ij}^p$ is the plastic strain increment, which is expressed in equation (17), as follows:$$\begin{array}{}\text{(34)}& \text{d}{\epsilon }_{ij}^e=C_{ijkl}^e\text{d}{\sigma }_{kl}.\end{array}$$

According to equations (17), (33), and (34), the total strain increment was derived as $\text{d}{\sigma }_{ij}=D_{ijkl}^e\text{d}{\epsilon }_{kl}-\text{d}\lambda \partial g/\partial {\sigma }_{is}{F}_{sj}$:$$\begin{array}{}\text{(35)}& \text{d}{\epsilon }_{ij}=C_{ijkl}^e\text{d}{\sigma }_{kl}+\text{d}\lambda \frac{\partial g}{\partial {\sigma }_{ik}}{F}_{kj}.\end{array}$$Equation (35) was simplified as follows:$$\begin{array}{}\text{(36)}& \text{d}{\sigma }_{ij}=D_{ijkl}^e\text{d}{\epsilon }_{kl}-\text{d}\lambda \frac{\partial g}{\partial {\sigma }_{is}}{F}_{sj}.\end{array}$$

Equation (40) could be rearranged as$$\begin{array}{}\text{(41)}& \text{d}\lambda =\frac{\partial f/\partial {\sigma }_{ij}{D}_{ijkl}^e\text{d}{\epsilon }_{kl}}{A_p+\partial f/\partial {\sigma }_{mn}{D}_{mnpq}^e\partial g/\partial {\sigma }_{su}{F}_{ut}}.\end{array}$$where $D_{ijkl}^e$ is the elastic stiffness matrix. dλ is also related to the hardening parameter $A$_p, so it could be rearranged as$$\begin{array}{}\text{(37)}& \text{d}\lambda =\frac{1}{A_p}\frac{\partial f}{\partial {\sigma }_{ij}}\text{d}{\sigma }_{ij}.\end{array}$$

According to the consistency equation df equal to 0 on the yield surface, the equation became$$\begin{array}{}\text{(38)}& \text{d}f=\frac{\partial f}{\partial {\sigma }_{ij}}\text{d}{\sigma }_{ij}+\frac{\partial f}{\partial M}\text{d}M=0,\end{array}$$where f is the yield function and $M$ is the hardening parameter.

Based on Equations (37) and (38), this led to the expression below:$$\begin{array}{}\text{(39)}& \text{d}\lambda A_p=-\frac{\partial f}{\partial M}\frac{\partial M}{\partial {\epsilon }_{ij}^p}\text{d}{\epsilon }_{ij}^p.\end{array}$$

After combining equations (36)–(39), the expression was given by$$\begin{array}{}\text{(40)}& \frac{\partial f}{\partial {\sigma }_{ij}}{D}_{ijkl}^e\text{d}{\epsilon }_{kl}-\frac{\partial f}{\partial {\sigma }_{mn}}{D}_{mnpq}^e\text{d}\lambda \frac{\partial g}{\partial {\sigma }_{su}}{p}_F{\overline{F}}_{ut}-A_p\text{d}\lambda =0.\end{array}$$

Equations (37) and (38) were combined as follows:$$\begin{array}{}\text{(42)}& A_p=-\frac{\partial f}{\partial M}\frac{\text{d}M}{\text{d}\lambda }=-\frac{\partial f}{\partial M}\frac{\partial M}{\partial {\epsilon }_{ij}^p}\frac{\partial g}{\partial {\sigma }_{ik}}{F}_{kj}.\end{array}$$

Because the hardening law took the plastic shear strain as the hardening parameter in the form of an incremental hyperbola, it could be obtained according to the flow law:$$\begin{array}{}\text{(43)}& A_p=-\frac{\partial f}{\partial M}\frac{\text{d}M}{\text{d}\lambda }=-\frac{\partial f}{\partial M}\frac{\partial M}{\partial {\epsilon }_q^p}\sqrt{{\frac{\partial g}{\partial q}{F}_{kj}}^2+{\frac{1}{q}\frac{\partial g}{\partial {\theta }_{\sigma }}{F}_{kj}}^2}.\end{array}$$

The hardening law equation (18) was substituted into equation (43) to obtain the hardening function as$$\begin{array}{}\text{(44)}& A_p=g{\theta }_{\sigma }p\frac{h_{s}G{M_p-M}^2}{p{M}_p{}^2}\sqrt{{\frac{\partial g}{\partial q}{F}_{kj}}^2+{\frac{1}{q}\frac{\partial g}{\partial {\theta }_{\sigma }}{F}_{kj}}^2}.\end{array}$$

4.5. Model Parameter Determination

There were four groups and seventeen model parameters in this research. The elastic parameters and the fractal parameters could be measured using the basic test data of the rockfill materials. The critical state parameters were tested with a triaxial tensile test. The state-related parameters were determined with the drained triaxial test. All model parameters are shown in Table 1. The constitutive model parameters were divided into four groups: elastic parameter, fractal parameter, critical parameter, and state parameter. The specific methods for determining these parameters were as follows.

Table 1

Parameters of the state-dependent constitutive model for rockfill.

Elastic parameter	Fractal parameter	Critical parameter	State parameter
$G_0=190$	l = 0.744	MCS = 1.727	d0 = 2.267
$v$ = 0.3	k = 0.008	$a^{\prime }$  = 0.032	cψ = 0.98
	β = 12.8	$b^{\prime }$ = 0.527	h1 = 0.46
		$c^{\prime }$ = 0.313	h2 = 0.78
		kh = 0.12	kp = 4.95
		λc = 0.14	kd = 0.458

4.5.1. Elastic Parameter

$G_0$ was the material constant. It can be obtained from the three-dimensional stress-strain curve of the rockfill material. $v$ was Poisson’s ratio reflecting the transverse deformation of the material and can be obtained by the triaxial compression test.

4.5.2. Fractal Parameter

$l$, $\kappa$, and $\beta$ were the fractal parameters. The initial fractal dimension $D$₀ and fractal dimension $D$ could be obtained from equation (2). Then according to equation (3), the fractal parameters could be obtained from fitting the fractal dimension $D$ and ${\sigma }_3/P_a$ through the binary linear regression analysis.

4.5.3. Critical Parameter

According to the measured initial void ratio $e_0$, the values of $a^{\prime }$, $b^{\prime }$, and $c^{\prime }$ were obtained from equations (21) and (22). For rockfill materials, the value of the parameter ${\lambda }_c$ was obtained by normalizing mean effective stress after obtaining the critical state line through triaxial test data. Regarding the other critical parameters $k_h$ and $M_{cs}$, for details, please refer to relevant research for determination [32, 34].

4.5.4. State Parameter

$d_0$ was the expansion coefficient which can be obtained from the dilatancy equation and conventional triaxial test of rockfill material. $k_p$ and $k_d$ were model parameters that can be obtained from equations (24) and (30) in the test, respectively. $c_{\psi }$ was defined as $c_{\psi }=\cos {\theta }_{\sigma }-{\theta }_{\text{d}\epsilon }$. It can be obtained from the stress Lode angle in the critical state and the principal stress space strain increment Lode angle. $h_1$ and $h_2$ were model parameters, please refer to the relevant research for details [34].

5. Test Verification

5.1. Introduction to the Experiment

The experiment adopted the consolidated drained triaxial test results of Cai [23] by using a large triaxial instrument. The test specimen was dolomitic limestone. The density of the particle was 2.77 g/cm³ and the particle size was less than 800 mm. Furthermore, the grading was good. The nonuniformity coefficients C_u and C_c were 35.48 and 1.35, respectively. There were four types of gradation for the rockfill samples and the relative densities were 0.65, 0.75, 0.90, and 1.0. To characterize the strength characteristics of the rockfill under different relative compactness and confining pressures, the triaxial shear test was conducted under four confining pressures: 300 kPa, 600 kPa, 1000 kPa, and 1500 kPa.

5.2. Experimental Characteristics and 3D Description of the Fractal Dimension

Tapias et al. found that grain crushing largely resulted in the deformation of rockfill particles under different pressures [35]. The experiment showed that when the contact force was small, contact-related crushing dominated, although equivalent volume division dominated beyond given yield stress. The experimental characteristics of the fractal dimension reproduced these observations and trends. Figure 2 shows the relationship curves between the grain gradation and the fractal dimension under the confining pressure of 1000 kPa for $D$_r equal to 0.90. The initially graded 1 to 4 samples showed obvious particle fragmentation after the test. The particle fragmentation was mainly concentrated in the range of 80 mm to 40 mm. The results of particle crushing were about 5%–7% higher than the average particle mass percentage before crushing. The fractal dimension of crushing was proportional to the confining pressure. The dependence of stress path was also one of the important characteristics of rockfill material properties [36, 37]. The large-scale triaxial compression test showed that the stress path had a hierarchical evolution effect on the critical state and the particle fragmentation of the rockfill materials.

[figures omitted; refer to PDF]

Besides, Guo et al. found that the particle breakage and particle size distribution had an important effect on the critical state line of the rockfill under stress [38–40]. The critical state line of rockfill was very similar to sand. The effective mean normal stress $p^{\prime }$, deviatoric stress $q$, and volume changes of the rockfill were no longer changing after reaching the critical state. Then, the continuous increment of the shear deformation would cause different degrees of damage. However, because of the different initial material properties (particle size, initial gradation, etc.) of the rockfill, the critical state was difficult to determine accurately. Therefore, we would approximate achieve the critical state of the rockfill when the deviatoric stress and volume deformation were stable in the shear test. So, it is generally believed that the stable state was when the shear strain was greater than 20%. Figure 3 shows that the critical state point of the rockfill material is linear and unique in the $p^{\prime }$∼$q$ space. This is consistent with the research results on critical state line of rockfill material [40–42].

[figure omitted; refer to PDF]

Figure 4 shows the $e\sim {p^{\prime }/p_a}^{\xi }$ curve of the rockfill. The slope of the critical void ratio of the four-gradation rockfill was essentially the same, but the y-axis intercepts were not the same, which reflected the compaction of the rockfill with different levels and different relative densities. The experimental results showed that the critical void ratio intercept was closely related to the initial void ratio and the fractal dimension. Figure 5 shows that the critical void ratio intercept had a good linear relationship with the initial void ratio and the fractal dimension, which further reflected the necessity of the fractal dimension when describing the void of the rockfill particles.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Figure 6(a) reflects the failure surfaces of the rockfill materials at different levels. The figure indicates that the gradation, fractal dimension, and stone-piling strength were correlated. The greater the inhomogeneity coefficient of the rockfill gradation was, the greater the fractal dimension was and the greater the yield strength was. Figures 6(b) and 6(c) show the variation rules of the strength parameters and the strength failure surface on the plane, respectively. The figures show that the strength failure lines on the π plane were similar shapes in equilateral curved triangles and they were symmetrical for the condition in which the orthotropic amplitude parameters a₁ and a₂ were equal. With an increase of the $A$_c or $B$_c, the strength failure line on the π plane increased along the coordinate axis. The increase of the strength parameter $A$_c essentially reflected the compactness of the arrangement of particles with the decrease of the void ratio. The aggregation of the particle action showed the characteristic of the stronger yield strength macroscopically. However, the increase of the strength parameter $B$_c made the yield strength decrease because the parameter $B$_c was directly proportional to the void ratio. When the amplitude a₁ and the amplitude a₂ changed, the yield strength in the vicinity of σ₃ and σ₂ changed in response.

[figures omitted; refer to PDF]

Figure 7 shows the relationship between the void ratio and the intermediate principal stress coefficient. It is shown in Figure 7(a) that the degree of anisotropy of the rockfill materials decreased with the increase of the intermediate principal stress coefficient from the microscopic characteristics of the rockfill particles. However, on the macro level, the critical void ratio decreased and the yield strength increased. Figure 7(b) shows the denser arrangement of rockfill particles on the microscopic level when the gradation nonuniformity coefficient and the initial fractal dimension were larger under the same pressure. The test results showed that the critical void ratio increased with the yield strength.

[figures omitted; refer to PDF]

5.3. Deformation Simulation of the Triaxial Drainage Test

The test results and the model simulation results of three kinds of rockfill gradation were compared at $D$_r = 0.75, as shown in Figure 8. It can be seen from the stress-strain curve that the constitutive equation established in this research could better predict the dilatancy characteristics of the rockfill. The gradation unevenness coefficient and the initial fractal dimension had a positive correlation to the impact of the rockfill swelling. Moreover, as deviatoric stress σ₁−σ₃ increased at the peak state, the softening of the rockfill became more obvious. The test results showed that the confining pressure σ₃ had an important influence on the process of the softening to the hardening of the rockfill. Specifically, the volume from contractive to dilative was increasingly obvious with the increase of the peak value of the deviatoric stress. The finer the gradation was (the initial fractal dimension $D$₀ was larger), the greater the peak value of the deviatoric stress σ₁−σ₃ was. At that time, it changed from the hardening type to the softening type, and the volume changed from contractive to dilative.

[figures omitted; refer to PDF]

6. Conclusions

Because of the research status of large particle sizes, compression crushing characteristics, and lack of relatively large-scale triaxial test results and the difficult description on the constitutive equation, this research established a constitutive model based on the test results of rockfill. Additionally, the research method for combining the macro and the micro effects of particle gradation with its anisotropy was described, along with the plastic state theory of the material state. The main conclusions were as follows:

(1) A constitutive model of rockfill was established based on property-dependent plastic potential theory for geomaterials. The analysis of the model could be used to characterize the particle gradation of the rockfill and the law of particle breakage under loading conditions. At the same time, the microscopic arrangement and macroscopic yield characteristics of the rockfill were connected by the fabric tensor. The hardening law could better describe the anisotropic mechanical response of the mesoscopic fabric, fractal dimension, particle breakage, and stress state to the rockfill. Modeling based on property-dependent plastic potential theory for geomaterials made the model give a better description of the mechanical characteristics of the principal stress axis rotation caused by complex stress conditions.

(2) A conventional rockfill triaxial drainage test was used to verify the model. The new model established the relationship between the mesoscopic particles and the stress state of the rockfill and used anisotropic state variables to reflect the anisotropy caused by the particle arrangement, such as initial anisotropy and stress-induced anisotropy. The model also considered the influences of the fractal dimension, particle fragmentation, gradation, and other factors on the critical state line of the rockfill. The performance of the model was verified by the test results in describing the critical state under various stress paths. In this research, the intensity change in three dimensions and its intensity change rule on the π plane were simulated, and the validity and the rationality of the model were verified with the simulation.

(3) The triaxial tests and simulations confirmed that the fractal dimension could reflect the crushing of the rockfill particles and the scale of the rockfill during the loading of the stress space internal force, as well as the gradation change of the rockfill and its critical void ratio and yield strength deformation. In this research, based on the fractal dimension theory, the two material properties of sand and rockfill were organically linked, which could not only quantitatively describe the particle grading, but also clearly express the void ratio strength parameters $A$_c and $B$_c.

(4) The experimental results showed that the gradation and the initial fractal dimension determined the material properties of the rockfill, the peak value of the deviatoric stress had an important effect on the softening of the rockfill, and the influence of the confining pressure was mainly in the process from contractive to dilative. The constitutive model accurately simulated this result. However, based on the limited experimental data and the need for practical application, the constitutive equation needs to be further improved.

Acknowledgments

This research was financially supported by the Key R & D Program of Ningxia Hui Autonomous Region Projects of International Cooperation and Exchanges (no. 2018DWHZ0084), Ningxia Science & Technology Innovation Leading Talent Project (no. KJT2019001), the National Nature Science Foundation of China (no. 51669027), and the National Key R & D Program of China (no. 2017YFC0504404), These supports are gratefully acknowledged.

References

[1] L. Zeng, X. F. Yao, J. H. Zhang, Q. F. Gao, J. C. Chen, Y. T. Gui, "Ponded infiltration and spatial-temporal prediction of the water content of silty mudstone," Bulletin of Engineering Geology and the Environment, vol. 79, pp. 5371-5383, DOI: 10.1007/s10064-020-01880-1, 2020.

[2] Z. M. He, Y. X. Liu, H. L. Tang, Y. H. Xing, H. B. Bian, "Experimental study on cumulative plastic deformation of coarse-grained soil high-grade roadbed under long-term vehicle load," Advances in Civil Engineering, vol. 2018,DOI: 10.1155/2018/8167205, 2018.

[3] J. H. Peng, J. H. Zhang, J. Li, Y. S. Yao, A. S. Zhang, "Modeling humidity and stress-dependent subgrade soils in flexible pavements," Computers and Geotechnics, vol. 120, 2020.

[4] X. X. Zhang, "Research on cyclic elasto-plastic constitutive model for rockfill materials and its application," 2015. Ph.D. thesis

[5] Z. Y. Cai, S. Y. Ding, Q. T. Bi, "Numerical simulation of strength and deformation characters of rockfill," Chinese Journal of Rock Mechanics and Engineering, vol. 28 no. 7, pp. 1327-1334, 2009.

[6] A. Varadarajan, K. G. Sharma, S. M. Abbas, A. K. Dhawan, "Constitutive model for rockfill materials and determination of material constants," International Journal of Geomechanics, vol. 6 no. 4, pp. 226-237, DOI: 10.1061/(asce)1532-3641(2006)6:4(226), 2006.

[7] Z. Z. Fu, S. S. Chen, Q. M. Zhong, Y. J. Zhang, "Modeling interaction between loading-induced and creep strains of rockfill materials using a hardening elastoplastic constitutive model," Canadian Geotechnical Journal, vol. 56 no. 10, pp. 1380-1394, DOI: 10.1139/cgj-2018-0435, 2019.

[8] J. H. Zhang, J. H. Peng, A. S Zhang, "Prediction of permanent deformation for subgrade soils under traffic loading in southern China," International Journal of Pavement Engineering,DOI: 10.1080/10298436.2020.1765244, 2020.

[9] W. L. Guo, L. Chen, "A stress-dilatancy relationship for rockfill incorporating particle breakage and intermediate principal-stress ratio," KSCE Journal of Civil Engineering, vol. 23 no. 7, pp. 2847-2851, DOI: 10.1007/s12205-019-0279-8, 2019.

[10] M. C. Liu, Y. F. Gao, X. M. Huang, "Study on elasto-plastic constitutive model of rockfills with nonlinear strength characteristics," Chinese Journal of Geotechnical Engineering, vol. 27 no. 3, pp. 294-298, 2005.

[11] E. Bauer, "Constitutive modelling of wetting deformation of rockfill materials," International Journal of Civil Engineering, vol. 17 no. 4, pp. 481-486, DOI: 10.1007/s40999-018-0327-7, 2019.

[12] S. H. Liu, Y. Sun, C. M. Shen, Z. Y. Yin, "Practical nonlinear constitutive model for rockfill materials with application to rockfill dam," Computers and Geotechnics, vol. 119,DOI: 10.1016/j.compgeo.2019.103383, 2020.

[13] Z. J. Wang, S. S. Chen, Z. Z. Fu, "Viscoelastic-plastic constitutive model for creep deformation of rockfill materials," Chinese Journal of Geotechnical Engineering, vol. 36 no. 12, pp. 2188-2194, 2014.

[14] E. L. Liu, S. S. Chen, G. Y. Li, Q. M. Zhong, "Critical state of rockfill materials and a constitutive model considering grain crushing," Rock and Soil Mechanics, vol. 32 no. 2, pp. 148-154, 2011.

[15] Y. Xiao, H. L. Liu, "Elastoplastic constitutive model for rockfill materials considering particle breakage," International Journal of Geomechanics, vol. 17 no. 1,DOI: 10.1061/(asce)gm.1943-5622.0000681, 2017.

[16] S. H. Liu, C. M. Shen, H. Y. Mao, Y. Sun, "State-dependent elastoplastic constitutive model for rockfill materials," Rock and Soil Mechanics, vol. 40 no. 8, pp. 2891-2898, 2019.

[17] H. L. Fang, Y. H. Cai, W. J. Wang, "State-dependent 3D multi-mechanism bounding surface model for rockfills," Chinese Journal of Geotechnical Engineering, vol. 40 no. 12, pp. 2164-2171, 2018.

[18] D. H. Liu, H. Chen, "Relationship between porosity and the constitutive model parameters of rockfill materials," Journal of Materials in Civil Engineering, vol. 31 no. 2,DOI: 10.1061/(asce)mt.1943-5533.0002598, 2019.

[19] A. Brito, J. R. Maranha, L. M. M. S. Caldeira, "A constitutive model for soil-rockfill mixtures," Computers and Geotechnics, vol. 95, pp. 46-56, DOI: 10.1016/j.compgeo.2017.09.009, 2018.

[20] G. H. Yang, G. X. Li, "Constitutive theory of soils based on the generalized potential theory," Chinese Journal of Geotechnical Engineering, vol. 29 no. 4, pp. 594-597, 2007.

[21] G. H. Yang, G. X. Li, "Mathematical foundation of constitutive models of geotechnical material and generalized potential theory," Rock and Soil Mechanics, vol. 23 no. 5, pp. 531-535, 2002.

[22] X. F. Li, L. Kong, M. S. Huang, "Property-dependent plastic potential theory for geomaterials," Chinese Journal of Geotechnical Engineering, vol. 35 no. 9, pp. 1722-1729, 2013.

[23] Z. Y. Cai, X. M. Li, L. Han, Y. F. Guan, "Critical state of rockfill materials considering particle gradation and breakage," Chinese Journal of Geotechnical Engineering, vol. 38 no. 8, pp. 1357-1364, 2016.

[24] Y. Xiao, M. Q. Meng, A. Daouadji, Q. S. Chen, Z. J. Wu, X. Jiang, "Effects of particle size on crushing and deformation behaviors of rockfill materials," Geoscience Frontiers, vol. 11 no. 2, pp. 375-388, DOI: 10.1016/j.gsf.2018.10.010, 2020.

[25] J. R. Zhang, B. W. Zhang, "Fractal pattern of particle crushing of granular geomaterials during one-dimensional compression," Advances in Civil Engineering, vol. 2018, 2018.

[26] J. J. Ji, Q. Yao, F. M. Wu, H. T. Li, "Evaluation of particle size distribution of granular blasting materials based on the fractal theory," Advances in Mechanical Engineering, vol. 11 no. 10,DOI: 10.1177/1687814019881562, 2019.

[27] F. X. Zhao, S. C. Chi, X. F. Mi, "Gradation evolution model based on particle breakage characteristics for rockfill materials," Chinese Journal of Geotechnical and Engineering, vol. 41 no. 9, pp. 1707-1714, 2019.

[28] X. F. Li, M. S. Huang, J. G. Qian, "Failure criterion of anisotropic sand with method of macro-meso incorporation," Chinese Journal of Rock Mechanics and Engineering, vol. 29 no. 9, pp. 1885-1892, 2010.

[29] X. F. Li, R. J. Li, J. H. Zhang, Q. Wang, "Three-dimensional strength characteristic of rockfill," China Journal of Highway and Transport, vol. 33 no. 9, pp. 54-62, 2020.

[30] K. J. William, E. P. Warnke, "Constitutive model for the triaxial behavior of concrete," Proceedings of International Association for Bridge and Structure Engineering, pp. 117-131, .

[31] X. F. Li, Q. Wang, J. F. Liu, W. Wu, F. C. Meng, "Quantitative description of microscopic fabric based on sand particle shapes," China Journal of Highway and Transport, vol. 29 no. 10, pp. 29-36, 2016.

[32] X. S. Li, Y. F. Dafalias, "Dilatancy for cohesionless soils," Géotechnique, vol. 50 no. 4, pp. 449-460, DOI: 10.1680/geot.2000.50.4.449, 2000.

[33] X. L. Lü, M. S. Huang, J. E. Andrade, "Predicting the initiation of static liquefaction of cross-anisotropic sands under multiaxial stress conditions," International Journal for Numerical and Analytical Methods in Geomechanics, vol. 41 no. 17, pp. 1724-1740, DOI: 10.1002/nag.2697, 2017.

[34] X. S. Li, "A sand model with state-dapendent dilatancy," Géotechnique, vol. 52 no. 3, pp. 173-186, DOI: 10.1680/geot.2002.52.3.173, 2002.

[35] M. Tapias, E. E. Alonso, J. Gili, "A particle model for rockfill behaviour," Géotechnique, vol. 65 no. 12, pp. 975-994, DOI: 10.1680/jgeot.14.p.170, 2015.

[36] Y. F. Jia, B. Xu, S. C. Chi, B. Xiang, D. Xiao, Y. Zhou, "Particle breakage of rockfill material during triaxial tests under complex stress paths," International Journal of Geomechanics, vol. 19 no. 12,DOI: 10.1061/(asce)gm.1943-5622.0001517, 2019.

[37] F. W. Ning, J. M. Liu, X. J. Kong, D. G. Zou, "Critical state and grading evolution of rockfill material under different triaxial compression tests," International Journal of Geomechanics, vol. 20 no. 2,DOI: 10.1061/(asce)gm.1943-5622.0001550, 2020.

[38] C. S. Chang, Y. B. Deng, "Modeling for critical state line of granular soil with evolution of grain size distribution due to particle breakage," Geoscience Frontiers, vol. 11 no. 2, pp. 473-486, DOI: 10.1016/j.gsf.2019.06.008, 2020.

[39] W. L. Guo, Z. Y. Cai, Y. L. Wu, Z. Z. Geng, "Estimations of three characteristic stress ratios for rockfill material considering particle breakagefill material considering particle breakage," Acta Mechanica Solida Sinica, vol. 32 no. 2, pp. 215-229, DOI: 10.1007/s10338-019-00074-x, 2019.

[40] G. Yang, Y. Jiang, S. Nimbalkar, Y. F. Sun, N. H. Li, "Influence of particle size distribution on the critical state of rockfill," Advances in Civil Engineering, vol. 2019, 2019.

[41] Y. Xiao, H. L. Liu, X. M. Ding, Y. M. Chen, J. S. Jiang, W. G. Zhang, "Influence of particle breakage on critical state line of rockfill material," International Journal of Geomechanics, vol. 16 no. 1,DOI: 10.1061/(asce)gm.1943-5622.0000538, 2016.

[42] S. H. Bian, B. Wu, Y. Z. Ma, "Modeling static behavior of rockfill materials based on generalized plasticity model," Advances in Civil Engineering, vol. 2019,DOI: 10.1155/2019/2371709, 2019.