A state-dependent hyperelastic-plastic constitutive model considering shear-induced particle breakage in granular soils

An elastic–plastic constitutive model considering particle breakage for simulation of crushable granular soils behavior is proposed. In the model, elastic strain rates are derived from a modified Helmholtz free energy function, and the influence of plastic shear work-induced particle breakage on the elastic properties of sand is taken into account as an elastic–plastic coupling mechanism. A stress ratio-driven mechanism is employed for calculation of the plastic strain rates. The proposed model is capable of tracking the evolution of the grain size distribution (GSD) due to shear-induced particle breakage. The evolving breakage index of Einav (2007) (J Mech Phys Solids 55(6):1274–1297, 2007) is interrelated to the plastic shear work to avoid overestimation of shear-induced particle breakage in loose sands. A direct comparison between the model simulations and laboratory data has been carried out for five series of drained/undrained monotonic and cyclic triaxial tests covering a wide range of initial states. For the sake of comparison, predicted behaviors from a hypoplastic constitutive model specially developed for crushable granular soils are also included. It is shown that the proposed constitutive model can provide reasonable predictions using a single set of parameters for each series of the laboratory data.

Particle breakage results in permeability reduction around perforations [26,85], increase in the settlement of rockfill dams [1], creep in piles embedded in sands [68], reduction in the drained peak friction angle [6,7,46], rapid long-runout motion of landslides [83] and plays a fundamental role in stress-strain behavior of brittle granular soils [87,104,108,114]. In the past, a wide spectrum of experimental studies has been conducted to investigate the effect of particle breakage on the location of the critical state line (CSL) in crushable granular soils. According to the findings, in the presence of particle breakage, the critical state is reached at lower void ratios [3,30,48,57,78,101,120,126].
There exist two foremost traditions to consider particle breakage in constitutive models for granular soils: (i) modification of yield surface, dilatancy, plastic hardening and enabling the possibility of non-unique CSL through introduction of elegant breakage-dependent free energy and energy dissipation functions e.g., [14,22,23,81,101], and (ii) phenomenological methods using breakage-driven relocation of CSL in the e vs. p 0 plane in conjunction with critical state compatible statedependent constitutive equations e.g., [19,39,47,69,109,115,118,121]. The former clever route guarantees thermodynamics consistency, however, its rich conceptual theoretical mechanics/mathematics corroborates some researchers prefer adoption of the second, but commonly simpler method. For instance, Einav [22,23] introduced the theory of breakage mechanics with the purpose of modeling the behavior of crushable granular material from micromechanical considerations within the framework of hyperplasticity, guaranteeing thermodynamical consistency. Some studies extended the theory of breakage mechanics by connecting energetics with the micromechanics and investigating the variation of permeability in cataclasite zones [81], with inclusion of porosity as a state variable [90] and combined it with the Cosserat continuum through an elastic upscaling incorporating Cosserat state variables [14]. Tengattini et al. [101] introduced a porositydependent extension to the theory of [22,23] wherein the critical state is predicted in an unforced natural way, rather than imposed a priori. Tengattini et al. [101] succeeded to derive yield function and dilatancy from the rate of dissipation, but accurate simulation of stress-softening behavior is still lacking [12]. The large tear drop-shaped yield function, the lack of kinematic hardening, and the absence of the state parameter of Been and Jefferies [4] as a simple, but effective means for distinction of the dense from loose states, and the lack of phase transformation e.g., Ishihara et al. [50] for modeling of initial contraction prior of the subsequent dilation in medium-dense and dense sand subjected to undrained shear prompt the requirement of coupling between hyperelasticity and the bounding surface plasticity for crushable soils.
A constitutive model in a hyperelastic-plastic frame for crushable granular materials is suggested here. Owing to simplicity of the basic constitutive equations, versatility in simulation of various aspects of the mechanical behavior of granular soils, clear physical meaning of the parameters in one hand and ease of calibration using both data of monotonic and cyclic tests on the other, the state-dependent bounding surface model proposed by Dafalias and Manzari [16] has been adopted as benchmark to establish a hyperelastic-plastic constitutive model for crushable granular soils. A generalized Jiang and Liu [52]-type Helmholtz free energy function has been employed to derive elastic strain rates as well as the pressure-dependent soil elastic moduli. Hyperplasic theories for the granular soils with pressure-dependent moduli can take into account stress-induced anisotropy of the elastic response in a natural unforced way and improve model predictions under cyclic loading e.g., [24,31,38]. A breakage index depending on plastic shear work (say function of shear stress and plastic shear strain rate) is proposed to quantify particle breakage during shear. The breakage effect on the elastic stiffness is not verified, and the formulation follows constitutive equations proposed by Einav [22]. Owing to the irreversible nature of particle breakage, the breakage acts as an elastic-plastic coupling variable. Comparisons between the proposed model predictions and experimental data indicate reasonable performance of the model. For the sake of comparison, predictions obtained from the hypoplastic model of Engin et al. [25] specially developed for simulation of crushable sands are also included.

Fundamentals of constitutive equations for crushable sands
A definition of an index quantifying the extent of particle breakage is crucial in the mechanics of crushable media. Various types of relative breakage indexes have been proposed to quantify breakage in crushable granular soils in terms of the breakage-induced evolution of the GSD curve [59,65,67]. Einav [22] introduced a breakage index (i.e., B) varying from zero at initial state (i.e., nil breakage) to 1 at the ultimate grading which the crushable soil finally will approach. B is mathematically defined as the ratio of the area between the current and initial GSDs (i.e., B t ) to the area between the ultimate and initial GSDs (i.e., B p ) [see Fig. 1]: in Eq. (1), D stands for the grain size. D m and D M indicate, respectively, the minimum and maximum particle sizes, and F 0 (D), F(D), and F u (D) represent the initial, current,  Fig. 1 Definition of the particle breakage index of Einav [22] and ultimate GSDs, respectively. Evolution of particle breakage would terminate asymptotically once a GSD approaches its ultimate fractal distribution e.g., [22]. For this reason, F(D) evolves from F 0 (D) to F u (D) through [22]: Stiffness and mobilization of shear strength in sand-fine mixtures depend on fines content, and contribution of the coarse and fine fractions in the load transfer mechanism is a function of coarse to fine size ratio e.g., [17,61,66,75,88,102]. Applying the discrete element method (DEM), Einav [22] suggested that the energy stored in soil grains depends on their size, and thus, larger particles have greater stored energy compared to the fines. When a crushable granular soil is subjected to shear, the breakage-induced fines surround larger particles, but do not participate in the load-bearing microstructure actively [22,127], consequently, they (i.e., fines) do not store energy [22]. Therefore, the elastic energy functions in addition to the effective stress or strain invariants should also incorporate some measures of GSD. According to the breakage mechanics theory of Einav [23], the Helmholtz free energy function can be expressed in terms of the elastic strains and breakage index by applying the gradation curve as the average weighted function on microscopic variables, through a statistical homogenization process as: where e e v and e e q are elastic volumetric and shear strains, respectively. H 1 ðe e v ; e e q Þ is the Helmholtz free energy function for nil particle breakage, and t, a proximity index indicating the distance between the ultimate and the initial GSDs, is determined as: where in D h i 2 0 and D h i 2 u are, respectively, the second-order moments of the initial and ultimate GSDs.
The effective stress variables in triaxial space can be calculated from the modified Helmholtz free energy function as: wherein p 0 ¼ r 0 a þ 2r 0 r À Á =3 and q¼ r 0 a À r 0 r are, respectively, the mean principal effective stress and shear stress (of note, r 0 a and r 0 r are the axial and radial principal stresses acting on a triaxial specimen). Further differentiation of Eq. (5) yields the rates of the effective stress invariants:  (7) necessitates that the total elastic strain rates become functions of _ e e v , _ e e q and _ B. The particle breakage is an irreversible phenomenon [see Guo and Zhu [33] for discussion] which itself is capable of affecting the elastic response of crushable sands through the participation of _ B, H; vB and H; qB in Eq. (7). Therefore, the total elastic volumetric and shear strain rates in Eq. (7) are decomposed into two terms with completely different physical interpretations: (i) the term (I) which is reversible upon effective stress reversal, and (ii) the breakage-dependent term (II) which is irreversible. Accordingly, the total elastic strain rates are irreversible with effective stress rates unless under nil particle breakage (say _ B ¼ 0) which occurs on the condition that sand behaves purely elastic. For this reason, B as an irreversible variable can be interconnected to a proper hardening variable to bear the role of an elastic-plastic coupling factor. In this sense, B increases progressively with plastic strains in the elasticplastic regime of the behavior and remains unchanged ( _ B ¼ 0) given that soil behaves purely elastic. In the modern continuum mechanics, free energy functions incorporate a supplemental internal variable depending on the hardening parameter(s) to account for the effect of elastic-plastic coupling e.g., [13,40,41]. Of note, the coupling between elasticity and plasticity has been addressed in the pertinent literature of geomechanics e.g., [13,28,29,31,40,41,64]. Using the concept of alternative decomposition of strains for coupled materials by Collins and Houlsby [13], the total volumetric and shear strain rates can be expressed as the elastic reversible and irreversible rates: In Eq. (8), _ e er x , _ e ei x and _ e i x with x 2 v; q f g are, respectively, the elastic reversible, elastic irreversible, and total irreversible strain rates. Using Eqs. (7) and (8) Rearranging the terms in the first part of Eq. (9) yields: Equation (10) resembles the theory of Graham and Houlsby [32] for anisotropic soils in the triaxial space: where G, K and J are the elastic shear, bulk and shear-volumetric moduli, respectively. J denotes the stress-induced anisotropy impact on the elastic behavior e.g., [5,23,31,37,38]. The consideration of the influence of anisotropy in the elastic domain plays a significant role in predicting shear strength of sands subjected to undrained cyclic loading e.g., [28,29,60,82] as well as of fine-grained soils e.g., [95,96,97,98,99,111]. Now, comparing Eqs. (10) and (11) results in the following elastic moduli: Using resonance column tests, Iwasaki and Tatsuoka [51] reported that at a fixed void ratio, the shear modulus is influenced by variations in the GSD curve. Wichtmann and Triantafyllidis [110] revealed that G decreases with increasing uniformity coefficient in granular soils. Given that the GSD varies with respect to the progression of particle breakage [3,35,59,65] then particle breakage can affect the elastic moduli of soils. As the Helmholtz free energy is a function of particle breakage according to Eq. (3), also the elastic moduli will be influenced by the change in GSD.
The hyper-elastic part has been conjugated to an elastoplastic theory which is the Simple Anisotropic Sand Model (SANISAND) proposed by Dafalias and Manzari [16]. Many studies extended the SANISAND version proposed by Dafalias and Manzari [16] by adding a second yield criteria e.g., [101], by taking into account the anisotropic behavior of soils e.g., [60], by omitting the elastic range of the behavior e.g., [63], or by adding a memory surface [70]. But, the original version is not only simple but also robust enough to predict the behavior of granular materials over a wide range of initial confining pressures and dry densities. Besides this, the 2004 version is still most widely used in geotechnical practice and also in the calculation of boundary value problems in research [56,71].
3 A hyperelastic-plastic model accounting for particle breakage and elastic-plastic coupling Jiang and Liu [52] suggested two Helmholtz free energy functions for nonlinear hyperelastic response of soils dependent on p 01=3 (i.e., Hertz-type elasticity) and p 01=2 (i.e., Goddard-type elasticity). Ashkar and Lashkari [2] suggested an extended Jiang and Liu [52] Helmholtz free energy function to obtain hyperelastic response dependent on p 0v with v 2 ½0; 1Þ. A revised form of the latter Helmholtz free energy function with the purpose of incorporating the particle breakage is adopted here (see Eq. 13): wherein p 0 0 is confining pressure at nil elastic strain and p ref ¼ 100 [kPa] is a reference normalizing stress. B in Eq. (13) remains unchanged when the soil behaves purely elastic. However, B evolves progressively in the elastic-plastic regime of the behavior to reach its asymptotic value at the critical state. Equation (12) in combination with Eq. (13) results in the following pressure-dependent hyperelastic moduli: where g ¼ q=p 0 is stress ratio, and G ¼ G 0 FðeÞ and K ¼ K 0 FðeÞ are non-dimensional parameters dependent on the void ratio e with FðeÞ ¼ ð2:97 À eÞ 2 =ð1 þ eÞ proposed by Hardin and Richart [36]. G 0 and K 0 are material parameters. In Eq. (14), h is defined as: Hence, Eq. (14) describes the hyperelastic moduli as pycnotropy (i.e., void ratio-dependent) as well as barotropy (i.e., pressure-dependent) functions. Of note, Eq. (14) signifies that Eq. (13) is capable of reproducing p 0v -dependent hyperelastic moduli through the entire range of v 2 ½ 0; 1Þ; however, the basic Helmholtz free energy functions of Jiang and Liu [52] are only compelling for v ¼ 1=3 and 1=2. Moreover, the evolving GSD as an elastic-plastic coupling phenomenon is captured through the explicit participation of B in the hyperelastic and Helmholtz breakage coupling moduli [see the second part of Eq. (9)]: Appendix 1 illustrates the evaluation of the hyperelastic formulation of the proposed model.

Yield function
A narrow open wedge-type yield function is adopted from Dafalias and Manzari [16]: wherein m defines opening of the yield function (per default: m = 0.001). Back-stress ratio, i.e., a, is a kinematic hardening variable indicating rotation of the bisector of the yield function with respect to the e vs. p 0 -axis see [16]. Soil behaves purely elastic (i.e., elastic strains become completely reversible) as long as the effective stress state is located inside the yield function [say f(g, a)\0], whereas plastic strains are generated once the effective stress state reaches the yield function [i.e., f(g, a) = 0] and attempts to change the back-stress ratio (say kinematic hardening variable).

CSL and state parameter
To investigate the critical state properties of crushable granular soils, results of a series of drained triaxial tests on Tacheng rockfill from [118,120] are discussed here. The tests were conducted at four initial void ratios (i.e., 0.189, 0.244, 0.285, and 0.317), and three initial confining pressures (i.e., 0.4, 0.8, and 1.6 MPa). Figure 2a demonstrates the CSLs of Tacheng rockfill in the e vs: ðp 0 =p ref Þ n plane with n ¼ 0:7. The CSL incorporating the effect of grain crushing in the e-p 0 -B space becomes [39]:  whereby n is a model constant, and k is the slope of the CSL in the e vs:ðp 0 =p ref Þ n plane. ''a'' is a model parameter controlling the pace of relocation of CSL in the e vs:ðp 0 =p ref Þ n plane with particle breakage. e i and e u are, respectively, soil parameters related to the initial (corresponding to B ¼ 0) and ultimate (as B ! 1) positions of the CSL in the e vs:ðp 0 =p ref Þ n plane. Increase in fines owing to breakage of coarser parent particles in crushable granular soils causes gradual changes in the intercept (i.e., e C ) and slope (i.e., k) of the CSL in the e vs:ðp 0 =p ref Þ n plane e.g., [3,127]. However, change in k is practically negligible in comparison with the change in e C e.g., [10,30,62,78]. Breakage CSLs can be established considering critical states with the same breakage indices, as illustrated in Fig. 2a and are used herein solely for the determination of e i and e u Recently, Lashkari et al. [62] suggested that change in k becomes noticeable only once the fines fraction becomes so large that the mechanical behavior of the mixture transfers from the ''fines-incoarse'' regime to the ''coarse-in-fine'' one. In Fig. 2b, the predicted e C vs: B curves using Eq. (18) [4] which is the Euclidean distance between the current soil state from the breakage-dependent CSL in the e vs: p 0 plane can be incorporated into the constitutive model: wherein e is the current void ratio and e cs is the critical state void ratio calculated from Eq. (18) corresponding to the current values of p 0 and B.

Plastic strain rates
A stress ratio-driven mechanism in triaxial space suggested by Dafalias and Manzari [16] is employed here for the calculation of the plastic shear strain rate: whereby L is the loading index representing the magnitude of the plastic shear strain rate. The Macaulay brackets In Eq. (20), s is determined according to the location of the stress ratio with respect to a as the bisector of the yield function, s=1 if g À a ¼ m and s = -1 when ag = m. The plastic hardening modulus K p is defined by e.g., [31,64]: wherein h and c h are model constants. a in is the initial value of a at the beginning of the most recent shear loading. The dilatancy, d, defined as the ratio between the rate of irreversible volumetric strain and the rate of irreversible shear strain is expressed as e.g., [31,64]: In Eqs. (23), M is the slope of the CSL in the q-p 0 plane. M d and M b are dilatancy and peak stress ratios, respectively. n d and n b are model constants, determining the impact of the state parameter on the dilatancy and peak stress ratios, respectively. Now, using Eqs. (12), 9(b), (20), and (22), the volumetric plastic strain rate is calculated from:

Rates of the effective stress invariants and loading index
The rates of the effective stress invariants are calculated using Eqs. (6), (8), and (12) as: Recalling that _ g ¼ ½ _ q À g _ p 0 =p 0 , substitution of Eqs. (24) and (25) (20)] renders the following relation for the loading index:

Breakage index
Magnitude and spatial distribution of particle level contact forces in consequence of imposed stress path and strain level are key controlling factors for an individual particle to crush e.g., [20]. It is believed that the element level plastic work can be implemented as an appropriate factor taking into account the combined influence of the effective stress path and plastic strain magnitude e.g., [53,69,84]. Void ratio also governs the stress-strain-strength response of the granular soils and affects the evolution of microstructure and breakage of particles under both drained and undrained conditions. Particles in loose assemblies have greater space for sliding, rolling, and soil skeleton adjustment than dense assemblies at the same effective stress level. This phenomenon results in larger volumetric strains, thus dissipating greater plastic work in loose granular soils. Table 1 shows the breakage index B, plastic work (i.e., W p ) and plastic shear work (i.e., W p q ) for loose and dense Tacheng rockfill specimens under drained condition. According to Table 1, B for dense specimens (e 0 = 0.189) is greater than that of the loose specimens (e 0 = 0.317) at the same initial confining pressure. On the other hand, W p for the dense specimen is less than that of the loose ones. Hence, the use of W p would overestimate B of  the loose specimens compared to the dense ones. On contrary, W p q data are consistent with the trending of B. Therefore, we propose to interrelate B to W p q , which reads: Figure 3 illustrates the variations of B with W p q and e 0 for Tacheng rockfill. A hyperbolic relationship between B and W p q normalized by the initial void ratio, i.e., W p q =e 0 , in the following form can fit reasonably the laboratory data: wherein x is a dimensionless model constant and changes pace of B with W p q =e 0 . To demonstrate the function of x, the predicted g vs: e 1 and e v vs: e 1 curves for x = 1000, 2000, 5000, and 10,000 [-] obtained from the constitutive model for a Tacheng rockfill specimen with e 0 = 0.244 sheared under drained condition from p 0 0 ¼ 1:6 ½kPa are plotted in Fig. 4. At a fixed axial strain, B increases with decrease in x. Accordingly, the peak stress ratio vanishes and the volume change response tends to contract with increase in B. Since B associated with elevated pressures is beyond the scope of this research, W p q , rather than W p , was adopted here to circumvent overestimation of the shearinduced breakage in loose sands. Implementation of a narrow but closed yield function [instead of the wedge-type yield function of Eq. (17)] together with a second plastic hardening mechanism that is activated under loading with limited or even nil change in g can further improve the constitutive model predictive capacity. Under the latter circumstance, e p v or work due to the plastic volumetric strain must be hired for the second yield mechanism accounting for particle breakage occurring under isotropic or K 0 compression. Using Eq. (28), rate of the breakage index (i.e., _ B) becomes: Equation (29) signifies that B remains unchanged as long as soil behaves purely elastic (of note, _ W p q ¼ 0 when soil behaves purely elastic) and increases progressively with W p q until the ultimate GSD is attained at B ¼ 1.  4 Calibration of the model The proposed model parameters can be classified into five different categories: (i) parameters related to the elastic response, (ii) CSL parameters, (iii) dilatancy parameters, (iv) plastic hardening modulus parameters, and (v) parameters related to breakage of particles (see Table 2). G 0 can be determined from experimental data of bender element and resonant column tests or by using the initial slope of the q vs: e q curves, where small strains prevail. G and v can be found using the data of the elastic shear modulus in the small strain levels such as results of bender element and resonant column tests. At extremely low shear strain levels, Eq. (14) yields: a b d c f Fig. 6 The proposed constitutive model predictions versus experimental data of Toyoura sand for: a and b four undrained tests on medium-dense specimens (e 0 = 0.735); c and d seven undrained tests on medium-loose (e 0 = 0.833) and loose (e 0 = 0.907) specimens; e and f six drained tests on loose to very loose specimens (data from Verdugo and Ishihara [106]) By taking logarithm from both sides of Eq. (30), it can be re-written as: Hence lnðG p ref Þ and v can be determined, respectively, as the intercept and slope of the best straight line fitted to the experimental data in the ln G vs: lnðp 0 =p ref Þ plane. K 0 % ð0:7 to 0:9Þ G 0 is a reasonable assumption for various sands.

Mean Effective Stress, p' [MPa]
Proposed Model Experimental data Proposed Model Experimental data a b According to Eq. (4), t is obtained by using initial and ultimate GSDs. In the lack of ultimate grading curve, is a satisfactory assumption [e.g., Einav [22]]. The proposed range ofis 2.5 to 2.8 [e.g., Sammis et al. [93]]. Otherwise, t can be derived by the calculation of the second-order moments [ D h i 2 of the initial and ultimate GSDs see Fig. 5a]. b 1 and b 2 can be determined through a curve fitting process by using the initial and ultimate GSDs along with Eq. (2). Figure 5b Fig. 2). Once e i and e u were determined, e U can be obtained using Eq. (18). Figure 5d represents the calculation of e U for the Kurnell sand. n d, can be determined using the results of monotonic undrained triaxial tests on medium-dense and dense specimens exhibiting phase transformation (i.e., certain effective stress state at which the initial contraction turns into dilation). According to Eq. (23), n d is the slope of the best line fitted to the phase transformation data in the lnðM d =MÞ vs: w d plane wherein M d and w d are, respectively, the values of g and w at the phase transformation. Likewise, the plastic hardening modulus becomes zero at the peak shear strength in the conventional drained tests on medium-dense and dense specimens. According to Eq. (23), n b is the slope of the best line fitted to the peak stress ratio data in the lnðM=M b Þ vs: w b plane wherein M b and w b are the peak stress ratio and the corresponding state parameter, respectively. Parts ''e'' and ''f'' of Fig. 5 show the determination of n b and n d for the Kurnell sand.
Using Eqs. (22) and (9) a for the conventional drained triaxial compression tests (s = ? 1), the following equation can be obtained e.g., [31]: a b d c Fig. 9 The proposed and hypoplastic models predictions for mobilization of stress ratio and volume change behavior of six Tacheng rockfill specimens sheared under drained condition from p 0 0 ¼ 0:4; 0:8 and 1:6 ½MPa: a and b three specimens with e 0 = 0.244; c and d three specimens with e 0 = 0.189 (data from [117,119] d 0 is determined using Eq. (32) for drained triaxial tests on medium-dense and dense specimens following phase transformation and prior reaching the critical state. Figure 5g presents determination of d 0 for the Kurnell sand. h and c h can be calculated by fitting the predicted to the q vs: e q curves of the undrained and drained triaxial tests covering a wide range of initial void ratios. x manages the pace of particle breakage owing to shear loading and can be interconnected to W p q =e 0 through Eq. (28). Figure 5h illustrates determination of x for the Kurnell sand. 5 The constitutive model predictive capacity

Drained and undrained monotonic and cyclic triaxial tests on Toyoura sand
Verdugo and Ishihara [106] carried out an extensive series of undrained and drained triaxial compression tests covering a wide range of initial states (i.e., p 0 0 ¼ 0:1 to 3:0 ½MPa and e 0 ¼ 0:735 to 0:996 ½À) on Toyoura sand specimens prepared using moist tamping method. A majority of the Toyoura sand particles were made of quartz, and accordingly, no tangible sign of particle breakage within the effective stress range investigated was traced by Verdugo and Ishihara [106]. In the testing program, each specimen was isotropically consolidated to the target p 0 0 value, sheared under undrained/drained condition until e 1 % 25 ½% was reached and thereafter, unloaded to the isotropic effective stress state. Using the model parameter listed in Table 2 Fig. 6. Without changing the model parameters for Toyoura sands, the undrained behaviors of two specimens subjected to cyclic shear from p 0 0 % 0:30 ½kPa are simulated and plotted against the corresponding laboratory data in Fig. 7. As can be seen, the proposed model is capable of capturing Toyoura sand behavior subjected to both monotonic and cyclic loadings with different initial confining pressures and void ratios to a satisfactory extent. Of note, particle Table 3 The hypoplastic model parameters for simulation of drained triaxial compression tests on Tacheng rockfill breakage for the experiments performed on Toyoura sand by Verdugo and Ishihara [106] was negligible and thus, e i and e u were set identical.

Drained triaxial tests on Tacheng rockfill
Tacheng rockfill collected from Shangrila County of China is a well-graded gravel with fines content of 1.8% and uniformity coefficient of 5.54 see [117,119]. Large-scale triaxial specimens were prepared in five consecutive layers. To form each layer, Tacheng rockfill was carefully dropped into the split mold and then compacted using a vibrator and later saturated through vacuum saturation. Specimens were consolidated isotropically under p 0 0 ¼ 0:4; 0:8 and 1:6 ½MPa and then sheared under conventional drained condition. In Figs. 8 and 9, the proposed model simulations obtained from the use of the model parameters listed in Table 2 for Tacheng sand are illustrated against the experimental data of initially loose to dense Tacheng rockfill reported by [117,119]. Engin et al. [25] revised the hypoplastic model proposed by von Wolffersdorff [107] in order to improve its ability to simulate the crushable granular soils behavior (see Appendix 2). The nonlinear incremental formulation of hypoplasticity adopts neither a yield function nor decomposition of strain into elastic and plastic parts. Engin et al. [25] incorporated a history dependency into the hypoplasticity formulation to capture the effect of particle breakage. For the sake of comparison, the predicted behaviors from the hypoplastic model are included in the figures. The model was implemented and calibrated with the parameters listed in Table 3. The Tacheng rockfill specimens exhibit more dilative behavior when crushing is negligible under low confining stress (e.g., p 0 0 ¼ 0:4 ½MPa). Increase in particle breakage in Figs. 8 and 9 caused the volume change behaviors become progressively contractive. The overall behavior is well reproduced by the proposed model. On the other side, the hypoplastic model tends to overestimate the peak stress ratio and failed to predict realistically both the dilative and contractive responses. For the drained tests carried out under p 0 0 ¼ 1:6 ½MPa, the evolved GSDs measured at the end of the experiments (say e 1 ¼ 15 ½%) are plotted against the model predictions in Fig. 10. Particle breakage led to a broadening of GSD and increased the fines content. The model in conjunction with Eq. (2) was successful to replicate the change in GSD from its initial one to those reached at the end of the tests.

Drained triaxial tests on Cambria sand
Cambria sand is a uniformly graded sandy soil with subrounded particles covering sizes from 0.83 to 2.00 mm. Specimens for the conventional drained triaxial compression tests were prepared using dry pluviation method. For six Cambria sand specimens with e 0 =0.52 sheared from p 0 0 ¼ 2:1; 4:0; 5:8; 8:0; 11:5 and 15:0 ½MPa, the proposed model predictions for the mobilization of shear stress with axial strain, volumetric strain with axial strain and GSDs at the end of experiments agree reasonably with the corresponding laboratory data in Fig. 11. Of note, the model parameters used in simulations illustrated in Fig. 11 are listed in Table 2.

Drained and undrained triaxial tests on Kurnell sand
According to Russell and Khalili [91], Kurnell sand (with uniformity coefficient of 1.83 and mean particle size of 0.31 mm) is a predominantly quartz sand collected from sand dunes at Kernell, Sydney, Australia. Using the model parameters reported in Table 2 for Kurnell sand, the mechanical behavior of six Kurnell sand specimens subjected to conventional drained shear under p 0 0 ¼ 0:76; 1:41; 2:39; 3:0; 5:7 and 7:8 ½MPa with e 0 = 0.615, 0.634, 0.680, 0.641, 0.654, and 0.661 is plotted against their corresponding data in Fig. 12. To predict undrained behavior of six Kurnell sand specimens in Fig. 13, all model parameters except h ¼ 90 ½À are considered the same as drained tests which are illustrated in Table 2. An explanation for this assumption is the fact that the strain rate and the size of the samples were different between drained and undrained tests. Moreover, enlarged lubricated end plates were used for conducting drained tests on Kurnell sand, while undrained tests were performed by using rough end plates. In Figs. 12 and 13, the proposed model shows its predictive capacity in simulation of 12 drained and undrained triaxial tests performed on sand specimens prepared through dissimilar methods within a wide range of initial void ratios and confining stresses.

Undrained cyclic triaxial tests on Aio sand
Aio sand (with uniformity coefficient of 2.74 and specific gravity of 2.64) collected from Yamaguchi prefecture in the south-west of Honshu Island located in Japan. Aio sand contains sub-angular to angular grains with the maximum and minimum void ratios equal to 0.958 and 0.582, respectively see [43,113]. Figures 14, 15 and 16 present the undrained cyclic triaxial tests performed by Wu et al. [113] and Hyodo et al. [43] on Aio sand. Continuous porewater pressure buildup caused reduction in the mean effective stress. Furthermore, axial strain was accumulated mainly toward the extension side. The mechanical behavior of both samples was well reproduced by the proposed constitutive model. The breakage-induced progressive downward relocation of the CSL in the e vs: p 0 plane causes an increase in state parameter and thereafter, a tangible decrease in plastic hardening modulus (thorough decrease in M b ). The presence of _ B as the third term in the nominator of Eq. (23) a and the breakage-induced decrease in K p together magnify generation of irreversible strains.
Consequently, particle breakage can affect the evolution of hysteresis loops under cyclic loads. Figure 17 demonstrates the evolution of GSD after experiments performed by Hyodo et al. [43] versus model predictions. Figure 18 presents the variation of breakage index with W p q =e 0 for cyclic tests performed by Hyodo et al. [43] on Aio sand. Of note, e i and e u for Aio sand have been determined using monotonic tests performed by Hyodo et al. [44] and cyclic tests reported by Hyodo et al. [43] and Wu et al. [113]. Regarding Figs. 14, 15 and 16, the proposed model simulations are in a satisfactory agreement with the experimental results.

Conclusion
In this paper, a new constitutive model able to simulate the mechanical behavior of granular materials with and without shear-induced particle breakage was proposed. The main contents are summarized as follows: a b d c • The hyperelastic strain rates were driven based on a generalized Jiang and Liu-type Helmholtz free energy function to guarantee the reversibility of the elastic strains. In order to account for the effect of particle breakage and coupling between plastic strains and hyperelastic properties, an additional hyperelastic-plastic coupling variable was incorporated into the Helmholtz free energy function. • A breakage index evolving with the plastic shear work was employed to quantify shear-induced particle breakage and to play the role of elastic-plastic coupling variable. It was realized that the plastic work overestimates the shear-induced particle breakage in loose sands. However, it was shown that a reasonable correspondence between the shear-induced plastic shear work normalized by the initial void ratio and shearinduced breakage exists. By relating the breakage index to a kinematic hardening parameter, coupling emerged into the formulation. This makes the pressure-dependent moduli evolve as a function of the hardening variable. • The numerical analyses are performed for a series of triaxial tests on both materials with permanent grains as well as with grain crushing. The comparisons indicate that the established model is applicable to capture stress-strain behavior of granular materials over a wide range of initial void ratios and confining pressures with the same set of model constants. Both the monotonic and cyclic loading have been considered. • A hypoplastic model is used for comparison purposes.
Comparisons revealed a discrepancy between predictions and experimental data considering the peak stress ratio and the dilatancy with the same set of material parameters. a b d c

. Evaluation of the hyperelastic formulation
One of the main advantages of employing a hyperelastic formulation is the energy conservative framework for any arbitrary closed-loop effective stress path. To bring an example in this regard, the closed stress cycle presented in Fig. 19a has been simulated. By considering the behavior as purely elastic, the elastic strain components are calculated using Eqs.7 and 14. G 0 = 50, K 0 = 0.75G 0 , v = 0.95, and t = 0.6 have been assumed for the numerical simulations. In Fig. 19b-   both the strains and the total work, which is positive for the assumed stress path. It should be noted that changing the sense of the rotation (say E?D?C?B?A instead of A?B?C?D?E) does not affect the magnitude of elastic strain components. In order to further evaluate the elastic part of the model, the contours of constant elastic strain proposed by Houlsby et al. [37] have been used. The suggested graphic depicts isolines (e P = ffiffi ffi 3 p =3 e v = constant; e Q = 3/2 e q = constant) in the isomorphic strain space. Figure 20 demonstrates such isolines for the proposed model. Unlike isotropic elasticity, the deviatoric portion of stress affects the isotropic portion of strain and vice versa. Such an anisotropic response makes the contours of constant e P and e Q to be curved. Apart the well-known advantages of the hyperelasticity, this framework imposes mathematical requirements which should be fulfilled. In the development of hyperelastic frameworks, the elastic stiffness should fulfill the following equations [42]: These requirements render the following relation, which should be fulfilled by the herein proposed model: The above equation should be considered for the calibration of the elastic parameters.
Appendix 2: A hypoplastic model for crushable granular soils [25] The model has been developed as an improved version of the grain crushing hypoplastic model of Rohe [89], which is based on the hypoplastic model by von Wolffersdorff [107]. The basic hypoplastic model can be described by a nonlinear tonsorial equation as: wherein r and e are effective stress and strain tensors, respectively. The fourth-order tensor L and the secondorder tensor N are, respectively, associated with the linear and nonlinear parts of the behavior. Both L and N are functions of effective stress and strain rate tensors. In general, in the hypoplastic framework, neither the definition of a yield surface nor a decomposition of strain rates is necessary. Furthermore, the nonlinear irreversible strains are defined by means of the degree of nonlinearity and the flow rule. Here, only the enhancement of Engin et al. [25] will be summarized. The attentive reader is referred to the cited references for further details about hypoplasticity.
Given the fact that grain crushing is irreversible, Engin et al. [25] the reference void ratios at the densest state, critical state, and loosest states, respectively, through: with: The subscript x denotes either minimum or maximum void ratio modification. In Eqs. (36) and (37), a p , a q , b p , b q , C u0 , c min , and c max are model parameters.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.