A Generalised Plastic Model for Gravelly Soils Considering Evolution of Void Ratio and Particle Breakage

Gravelly soils exhibit complicated mechanical behaviours closely related to particle breakage and relative density state. To better capture the mechanical responses of gravelly soils, a generalised plastic model considering evolution of void ratio and particle breakage was developed within the framework of critical state soil mechanics. In the model, particle breakage effect was described by incorporating breakage index to deviate the critical state line off the ideal position. A differential equation relating increment of void ratio to variation of volumetric strain was used to depict the evolution of current void ratio. It indirectly reflected the relative density state of gravelly soils. The model was applied to conducting numerical simulations for a series of triaxial tests on four types of gravelly soils. Comparisons between the test data and the modelling results indicated that considerations of void ratio evolution and particle breakage could better simulate the stress-dependent dilatation/contraction behaviours of gravelly soils.


INTRODUCTION
Gravelly soils (including coarse-grained aggregates, rockfill materials and railway basalts) are widely used as construction materials in the various geotechnical projects (e.g. roadbeds and dams). Over the past several decades, researchers have made efforts to investigate the mechanical properties of gravelly soils [1]. Experimental studies have indicated that particle breakage greatly affects the mechanical behaviours of gravelly soils [2,3]. In a microstructure view, particle breakage changed the grain distributions. Then, the mechanical responses of gravelly soils have changed remarkably with the reconfiguration of particle grading [3].
Based on the fundamental observations, particle breakage has been generally taken into account to build the constitutive models for gravelly soils [3][4][5][6]. In the constitutive descriptions, effect of particle breakage was described using a variable of breakage index. Particle grading evolution was often used to define the breakage index [2,3]. Investigations have indicated that particle breakage index is correlated to plastic-dissipation energy and formulated as a hyperbolic function of plastic work [5].
Besides the constitutive frameworks of traditional elastoplasticity [7], generalised plasticity [8] and critical state soil mechanics [9] have been often combined to construct the constitutive models of gravelly soils. In this group of models, as an important state parameter, particle breakage index has been incorporated into the equation of critical state line [3,5,6]. Such mathematical treatments can describe the translation effects of particle breakage on the critical state line of gravelly soils [5]. On the other hand, generalised plasticity is a type of simple constitutive framework, due to without consideration on the complex failure surfaces of geomaterials [8]. Thus, a number of generalised plastic models have been developed for gravelly soils based on the critical state concept [5,[10][11][12]. These models considering particle breakage could well capture the main mechanical behaviours of gravelly soils, including stress-dependent volumetric dilatancy and nonlinear strength and deformation.
However, most of these models have not made clear expositions on the variations of void ratio during volume changing process. As we all know, current void ratio of a soil mass closely relates to and evolves with the volume deformation. By identifying the current, maximum and minimum void ratio values, the relative density state of soil mass could be estimated. As for gravelly soils, relative density state (dense or loose sates relative to critical state) decides whether they exhibit volume dilatation or contraction behaviours [12]. Considering the relationship between void ratio and relative density state, the evolution of void ratio should be included in the constitutive descriptions. In this way, the changing process of relative density state synchronous with the variation of volume deformation can be captured precisely. Therefore, besides particle breakage, it is reasonable to take the evolution of void ratio into account to build the constitutive models for gravelly soils. In this respect, there are some research attempts of considering evolution of void ratio in the constitutive models of granular soils [6,13].
In this study, an evolution equation of void ratio was incorporated into a generalised plastic model that was established within the framework of critical state soil mechanics. In the model, current void ratio and particle breakage index were regard as two state-dependent variables updating with the variation of volumetric strain and shifting the critical state line, respectively. The model was used to simulate the mechanical responses of four types of gravelly soils during triaxial tests. By analysing the simulation results, constitutive modelling effectiveness of considering evolution of void ratio and particle breakage was evaluated.

Generalised Plasticity Framework
In this study, generalised plasticity framework firstly proposed by Pastor and Zienkiewicz [8] was used to build the constitutive model for gravelly soils. In a generalised plasticity model, incremental stress-strain relationship is expressed as follows: where dσ and dε represent the increments of stress σ and strain ε, respectively. D ep is the elastoplastic stiffness matrix, explicitly expressed as follows: where D e is the elastic constitutive stiffness matrix. g and f are the normalised plastic flow and plastic loading direction, respectively. H denotes the plastic modulus.

Elastic Behaviour
Nakai [14] provided the following Eqs 3, 4 to describe the elastic and plastic volumetric strains of sands under isotropic consolidation.
where p (σ 11 + σ 22 + σ 33 )/3 is the mean stress. ε v ε 11 + ε 22 + ε 33 denotes the volumetric strain. c t and c e are the compression and resilience indexes. m is a material constant. p 0 represents the initial value of mean stress. p a 0.098MPa denotes the engineering atmosphere pressure. Using Eq. 3 and Poisson's ratio ], the elastic bulk modulus K and shear modulus G can be defined as follows:

Critical State Considering Particle Breakage
As for gravelly soils, particle breakage remarkably changes the critical void state. This phenomenon was observed in triaxial tests [15,16]. In this study, a sigmoid function suggested by Liu et al. [3] was adopted to define the critical state line (CSL) influenced by particle breakage. That is: where e c means the critical-state void ratio, e max and e min are the maximum and minimum critical state void ratios of a gravelly soil mass, respectively. λ and ξ are the material constants. B r is the particle breakage index. In this study, a function proposed by Lade et al. [17] is adopted to link particle breakage index B r to plastic dissipated energy w p . That is: where a and b are both the material constants. w p denotes the plastic dissipated energy that is calculated using:

Evolution of Void Ratio
In soil mechanics, it is always assumed that the volume of solid phase V s is constant during loading and volume change is only related to the void phase V v . Assuming V s 1, void phase volume V v is equal to void ratio e. Thus, the change of void ratio causes an increment of volumetric strain. This relationship is formulated with: Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 767821 The second formula in Eq. 9 is a differential equation. It has a solution expressed as follows: In Eq. 10, the parameter A can be determined using the initial condition. That is ε v 0, e e 0 0A 1 + e 0 . The current void ratio e, however, is not always just onto the CSL, but usually above or below the CSL. To identify the relative density state of soils, Been and Jefferies [18] proposed the following state parameter ψ: where the state variable ψ describes the evolution of dilatation curve from contraction line towards the critical state line. ψ < 0 indicates that the relative density state of a gravelly soil is denser than the critical state. The gravelly soil may undergo a volumedilatation behaviour. ψ > 0 means the gravelly soil is in the looser state than the critical state. The gravelly soil may exhibit a volume-contraction behaviour.

Nonlinear Strength and Dilatancy Behaviours
Gravelly soils generally exhibit significant nonlinear strength behaviours. A power function is widely used to describe the nonlinear strength behaviours of geomaterials [3,19]: where q 3(σ ij − pδ ij )(σ ij − pδ ij )/2 is the deviatoric stress. Correspondingly, σ ij is the second-order stress tensor. δ ij is the Kronecker tensor, with i j, δ ij 1; i ≠ j, δ ij 0. g(θ) is a shape function of the failure criterion on the deviatoric plane. For simplification, g(θ) 1 is used in this study.
In Eq. 12, M c0 and n f are both material constants. Using Eq. 12, the critical stress state ratio M c is defined as the tangent slope of the strength envelope, that is: The dilatation stress ratio M d q d /p d is used to describe the stress state of transition phase point, which locates the stress threshold for the onset of volume-dilatation behaviour. Following the studies of Dafalias and Manzari [20], the dilatation stress ratio M d in this work is defined as a function of the state variable ψ: where κ is a material constant for describing dilatation behaviour.
The dilatation coefficient d g is related to the shear stress ratio η q/p and the dilatation stress ratio M d using the following equation [11]: where c p is the plastic component of the equivalent shear strain

Plastic Flow and Loading Direction
In the p − q stress plane, the plastic flow g is given by Pastor et al. [8]: As the nonassociated flow rule is assumed, the plastic loading direction f is different with the plastic flow g. For simplification, f can be expressed in a formula similar to g, as follows: where d f is defined as the loading direction factor, which is a state variable related to critical state. Similar to d g , d f is defined as:

Plastic Modulus
Under the isotropic compression condition, there is q 0 and dq 0. The incremental relationship between plastic volumetric strain and mean stress can be simply expressed as follows: Additionally, the plastic volumetric strain increment can be determined by differentiating Eq. 4, that is: Combination Eqs 19, 20 leads to the formula of plastic modulus H: Eq. 21 is only suitable for the isotropic compression condition. Inspired by Chen et al. [10], a formulation of plastic modulus H, suitable for more general cases, was developed. That is expressed as follows:

MODEL CALIBRATION AND SIMULATIONS
This constitutive model has 15 material constants and can be classified into six groups. They can be determined using triaxial compression tests.

Compression and Elastic Constants
Material constants c t , c e and m can be determined using loading/ unloading ε v ∼ p curves under various isotropic compression stress states. In particular, by analysing the unloading curves, ε e v and ε p v can be accurately distinguished from ε v . Then, using Eqs 3, 4 to fit the test data pairs (p, ε e v ) and (p, ε p v ), respectively, the values of c t , c e and m can be determined.
Additionally, using the elastic constants: c t , c e and m, the value of K can be determined. Then, using K and Poisson's ratio ], the value of G can be calculated.

Initial void Ratio
Initial void ratio e 0 can be determined using the following formula: where G s is the specific weight of soil particles. ω and ρ 0 are the water content and density of soil mass, respectively. ρ w denotes the water density.

Critical state Constants
The critical state constants e max , e min , λ and ξ are determined from the triaxial tests under small confining pressures, which no obvious particle breakage occurred. Eq. 6 (with B r ≈ 0) is used to fit a number of data pairs (e c , p). Then, the values of e max , e min , λ and ξ can be determined.

Strength constants
The strength constants: M c0 and n f can be determined by fitting the strength data (p, q) at the failure states using the strength criterion in Eq. 12. Because of nonlinearity of Eq. 12, at least three triaxial compression tests under different confining pressures are needed.

Dilatation constants
Using the test data of plastic volumetric and shear strains, d g can be calculated by d g dε p v /d c p . Certainly, it is difficult to depart the plastic strain increment dε p from the total strain increment dε. For simplification, it is assumed that d g ≈ dε v /d c. By relating d g and η q/p, then using Eq. 15, α can be determined.
Using the test data at the dilatation point, the dilatation stress ratio M d can be calculated. With the value of initial void ratio e 0 , Eq. 10 is used to compute the current void ratio e. Then, using Eqs 6, 11, the state variable ψ can be determined. In this way, the dilatation constant κ can be determined using the formula ln(M d /M c ) κψ to fit the test data pairs [ln(M d /M c ), ψ].

Particle breakage Parameters
Using Eq. 8, the plastic dissipated energy w p can be calculated. By comparing the particle size gradation before and after the tests, particle breakage index B r can be estimated [2,3]. Then, using Eq. 7, the particle breakage parameters a and b can be determined.
If B r cannot be explicitly estimated, an available method is combining Eqs 6, 7 to fit a number of test data pairs (e c , p, w p ) on the CSL to indirectly determine the values of a and b.

Performance of the Proposed Model
The proposed model was used to simulate the mechanical behaviours of four gravelly soils: Hekouchun rockfill [21], Wudongde rock-soil aggregate [22] and crushed latite basalt [23]. The material constants of the gravelly soils are given in Table 1. The test data and simulation results are shown in Figures  1-3. In these figures, the scatter symbols and curves represent the test data and modelling results, respectively.
Two factors: particle breakage and evolution of void ratio were taken into account to formulate the constitutive model. To observe their influences on the modelling effectiveness, three computing conditions were comparatively analysed. They are termed as: A "EVR + PB", which represent that the modelling predictions consider the compound effects of evolution of void ratio and particle breakage. B "PB_NoEVR", which represent that the modelling predictions only consider the effect of particle breakage. C "EVR_NoPB", which represent that the modelling predictions only consider the evolution of void ratio.
In the model, the shapes of numerical modelling q ∼ ε a curves are mainly controlled by the material constants: c t , c e , m, ], M c0 and n f . Thus, EVR + PB, PB_NoEVR and EVR_NoPB obtain almost the same modelling predictions of q ∼ ε a curves.

Hekouchun rockfill
Hekouchun rockfill material was a type of dolomitic limestone, which would be adopted to construct the Hekouchun reservoir dam. Cai et al. [21] conducted a series of large-scale consolidated and drained triaxial shear tests on Hekouchun rockfill material with an initial relative density D r 60%. The model parameters of the Hekouchun rockfill material are summarised in Table 1. The test data and modelling results are shown in Figure 1.
As seen from Figure 1 in both stress-strain and volume deformation responses, there is a good agreement between the test data and modelling predictions. Figure 1A shows that the model is capable of capturing the deformation and strength behaviours of Hekouchun rockfill material over a range of confining pressures from 0.3 to 1.5 MPa. Figures 1B,C show the testing volumetric strains and modelling results predicted by EVR + PB, PB_NoEVR and EVR_NoPB. It is found that, on the whole, EVR + PB achieves better simulations, but PB_NoEVR and EVR_NoPB overestimates and under-predicts the volumetric contraction strains, respectively.
The evolution curves of current and critical-state void ratios modelled by EVR + PB are shown in Figure 1D. The ideal CSL without considering particle breakage is plotted using a red dash curve. It is found that the realistic CSLs are influenced by particle breakage and gradually depart from the ideal CSL with the increase of mean stress. When the confining pressure increases, the volume-contraction behaviour of Hekouchun rockfill material becomes more prominent. The evolution rule of current void ratio is consistent with the volume changing rules.   Figure 2. Figure 2A shows that the modelling stress-strain relations exhibit good agreements in accordance with the test data. Figures  2B,C show that EVR + PB provides more accurate modelling predictions. PB_NoEVR overestimates volume contraction strains. EVR_NoPB tends to over-predict the volumedilatation behaviour. In Figure 2D, it is found that EVR + PB modelling curves of current and critical-state void ratios have an intersection point at p/p a 25. This indicates that, at this mean stress state, the volume deformation transits from dilatation to contraction. Therefore, it is reasonable to deduce that, under the confining pressure of 1.6 MPa, Wudongde rock-soil aggregate undergoes a phase change from relative dense to loose state.

Crushed latite Basalt
Salim and Indraratna [23] used a series of confining pressures: 0.05, 0.1, 0.2, and 0.3 MPa to conduct large-scale triaxial tests on a type of crushed latite basalt. Although the confining pressures were relatively lower, the volumetric strains had still undergone a transition from dilatation to contraction. Table 1 lists the material constants of the crushed latite basalt. Figure 3 shows the test data and modelling predictions. With an increase of confining pressure, PB_NoEVR predicts larger volumetric contraction strains. This phenomenon may be related to that PB_NoEVR overestimates the gap between the current void ratio and the CSL. Figure 3D shows the EVR + PB modelling curves of current and critical-state void ratios. Upward-bending evolutions of current void ratios indicates that volume-dilatation responses occur under the confining pressures of 0.05, 0.1, and 0.2 MPa. Under the confining pressure of 0.3 MPa, when the mean stress ratio p/p a > 7.9, volume deformation of the crushed latite basalt transits from dilatation to contraction.

CONCLUSIONS
In this study, a generalised plastic model considering evolution of void ratio and particle breakage was developed for gravelly soils. The model was used to simulate the mechanical behaviours of four gravelly soils. By comparing three computing conditions: EVR + PB, PB_NoEVR, FIGURE 2 | Model predations and experimental results of the Wudongde rock-soil aggregate (test data were from [22] and EVR_NoPB, the effects of void ratio evolution and particle breakage on the modelling predictions were analysed. Following conclusions can be summarised: 1. Because the evolution of void ratio and particle breakage are not considered in the constitutive descriptions of stiffness and strength behaviours, the three computing conditions: EVR + PB, PB_NoEVR, and EVR_NoPB predict almost the same relations between deviatoric stress and axial strain. 2. Under the PB_NoEVR computing condition, current void ratio is constant. In this condition, the proposed model overestimates the volume-contraction strains. The higher the confining pressure, PB_NoEVR will more greatly overestimate the volume-contraction strains. 3. Under the EVR_NoPB computing condition, the ideal CSL is used and current void ratio evolves with the variation of volumetric strain. Overall, EVR_NoPB underestimates the volume-contraction strains and overestimates volumedilatation strains. 4. In general, the model considering EVR + PB achieves relatively better simulation results for the volume dilatation/contraction responses of gravelly soils. It is reasonable for us to believe that the constitutive model considering EVR + PB could better reflect the volume deformation and failure mechanism of gravelly soils.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
Conception and design of study was mainly contributed by RW. Numerical simulations was mainly contributed by JZ, YL. Drafting the manuscript was mainly contributed by JZ.

FUNDING
The work presented in this paper were financially supported by Open  Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 767821