# Fractional Plasticity Model for Coarse Aggregates Incorporating Particle Breakage under Triaxial Loading

Authors: Jiangong Yang; Yifei Sun
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

Fractional plasticity has been originally proposed for capturing the state-dependent nonassociated behaviour of sand, with no use of state variables and additional plastic potentials. However, for coarse aggregates, significant particle breakage would occur, which could influence the stress-strain behaviour of the material, by shifting the critical state line in thee−lnp′ plane. In this note, an extension of the fractional plasticity for constitutive modelling of coarse aggregates is made by incorporating particle breakage under triaxial loading. The developed model is validated by simulating a series of laboratory test results of different coarse aggregates under triaxial loading, where a good agreement between the model simulations and test results is observed.


---

## Body

## 1. Introduction

Coarse granular aggregates, e.g., rockfill and ballast, have been widely used in constructing public facilities, due to its high mechanical qualities [1–4]. However, it has been long found that when subjected to external loads, particle breakage would occur in coarse aggregates, which significantly influenced the strength and deformation behaviour of the material [5–9]. For example, reduced volumetric dilatancy usually took place with increasing particle breakage [7]. The Mohr–Coulomb strength decreased as the stress level increased because a significant extent of particle breakage occurred in the material. Most importantly, particle breakage would change the critical state behaviour of coarse aggregates. For example, Yu [10] reported that particle breakage would shift down the critical state line (CSL) of rockfill in the e−lnp′ plane, where e and p′ are the void ratio and mean effective principal stress, respectively. Similar observations have been also reported elsewhere [11]. However, traditional constitutive approaches for soils were mainly based on linear critical state lines, which may not be suitable for coarse granular soils that exhibit high extent of particle breakage. To solve this problem, a series of elastoplastic models incorporating the concept of particle breakage have been proposed [12, 13]. However, an additional plastic potential was usually required for unified modelling of the nonassosiated stress-strain behaviour of coarse soils. Recently, a novel approach, the fractional plasticity [14–17], has been developed for modelling the state-dependent nonassociated behaviour for geomaterials, where no use of the additional plastic potentials was required. However, this approach did not consider the effect of particle breakage on the constitutive behaviour of highly crushable soils.Therefore, this study attempts to extend the original fractional plasticity for crushable coarse aggregates, by incorporating the concept of particle breakage. For the sake of simplicity, only homogenous and isotropic coarse aggregates under triaxial loads are considered here. For more complex loading conditions, one can refer to Sun et al. [15].


## 2. Fractional Plasticity considering Particle Breakage

### 2.1. Plastic Loading Direction

Following previous studies [14, 15], the loading surface (f) can be formulated as(1)f=qM2+p′2−p′p0′=0,where p0′ denotes the intercept of f with the p′-axis, M is the critical state stress ratio, [p′=σ3′+2σ3′/3] and q [=σ1′−σ3′] are the mean effective and deviator stresses, respectively, and σ1′ and σ3′ are the first and third effective principal stresses, respectively. The loading direction (n) can be obtained from the first-order derivatives of f, such that(2)n=nv,nsT=M2−η24η2+M4−2M2η2+η4,2η4η2+M4−2M2η2+η4T.


### 2.2. Plastic Flow Direction

Instead of using an additional plastic potential, the plastic flow direction (m) in fractional plasticity [14, 15] is directly obtained by conducting the fractional derivatives of f, such that(3)m=mv,msT=dg1+dg2,11+dg2T,where(4)dg=∂αf/∂p′α∂αf/∂qα,=M1+αp′−pc′+2−αpc′−p0′/2q−qc+2−αqc,in which α=expk1−p′/p0′, is used, considering the initial isotropic stress-state, where q = 0 and dg⟶+∞. k is a material constant. Note that the definition of the fractional derivatives can be found in appendix. A dependence of the plastic flow direction on the material state is observed in equation (4), where not only p′ and q but alsop′−pc′ and q–qc determine the value of dg. qc can be calculated from the relative position between the current stress point and CSL in the p′−q plane [14], such that(5)qc−qp′−pc′=q−qcpc′−p′=M,while pc′ can be obtained from the CSL in the e−lnp′ plane. As shown in many studies [7, 10], particle breakage would shift down the CSL. Therefore, a CSL with evolving intercept should be used for constitutive modelling. Considering the fact that both strain variation [18] and stress variation [2] could induce significant particle breakage, an energy-based CSL suggested by Liu and Zou [11] is used in this study for determining pc′:(6)pc′=prexpeΓ0−eΓ0wpawp+b−e/λ,where pr = 1 kPa is the unit pressure, λ is the gradient of the CSL in the e−lnp′ plane, e is the current void ratio, wp=∫σij′dεij is the plastic energy density, σij′ and εij are the stress and strain tensors, respectively, and eΓ0, a, and b are material constants. Through equation (6), the effect of particle breakage on the stress-strain behaviour of coarse aggregates can be thus considered.


### 2.3. Hardening Modulus

The plastic modulus,H, can be obtained by using the consistency condition together with the evolving law, i.e., Δh=Lh¯, for the hardening parameter (h). Here, Δ indicates increment. Hence,(7)H=−∂f∂hh¯,where L=mTΔσ′/H and h¯ can be defined by the law for deviatoric hardening. During loading, f would evolve until reaching the critical state, where the evolution should be bounded by the peak stress surface. Therefore, one can write(8)h¯=χGMp−ηq,where G=G02.97−e2/1+ep′pa is the shear modulus [14]. G0 is a material constant. pa=100kPa is the atmospheric pressure. χ is a model parameter and can be a function of the void ratio, e.g., χ=χ11−e, for better model performance. χ1 is a model constant. Mp=M1+βpc/p′−1 is the peak stress ratio [19], in which β is a material constant. Further substituting equations (1) and (8) into equation (7), H can be obtained as(9)H=χGMpη−1.In addition to the plastic strain, coarse aggregates may also experience elastic strain during loading, which can be described via Hooke’s elastic law. Accordingly, the elastoplastic relation for coarse aggregates can be provided as(10)Δεv=Δp′K+nvHΔp′+nsHΔqmv,Δεs=Δq3G+nvHΔp′+nsHΔqms,where the bulk modulus K=2+2νG/3/1−2ν and ν is Poisson’s ratio.


## 2.1. Plastic Loading Direction

Following previous studies [14, 15], the loading surface (f) can be formulated as(1)f=qM2+p′2−p′p0′=0,where p0′ denotes the intercept of f with the p′-axis, M is the critical state stress ratio, [p′=σ3′+2σ3′/3] and q [=σ1′−σ3′] are the mean effective and deviator stresses, respectively, and σ1′ and σ3′ are the first and third effective principal stresses, respectively. The loading direction (n) can be obtained from the first-order derivatives of f, such that(2)n=nv,nsT=M2−η24η2+M4−2M2η2+η4,2η4η2+M4−2M2η2+η4T.


## 2.2. Plastic Flow Direction

Instead of using an additional plastic potential, the plastic flow direction (m) in fractional plasticity [14, 15] is directly obtained by conducting the fractional derivatives of f, such that(3)m=mv,msT=dg1+dg2,11+dg2T,where(4)dg=∂αf/∂p′α∂αf/∂qα,=M1+αp′−pc′+2−αpc′−p0′/2q−qc+2−αqc,in which α=expk1−p′/p0′, is used, considering the initial isotropic stress-state, where q = 0 and dg⟶+∞. k is a material constant. Note that the definition of the fractional derivatives can be found in appendix. A dependence of the plastic flow direction on the material state is observed in equation (4), where not only p′ and q but alsop′−pc′ and q–qc determine the value of dg. qc can be calculated from the relative position between the current stress point and CSL in the p′−q plane [14], such that(5)qc−qp′−pc′=q−qcpc′−p′=M,while pc′ can be obtained from the CSL in the e−lnp′ plane. As shown in many studies [7, 10], particle breakage would shift down the CSL. Therefore, a CSL with evolving intercept should be used for constitutive modelling. Considering the fact that both strain variation [18] and stress variation [2] could induce significant particle breakage, an energy-based CSL suggested by Liu and Zou [11] is used in this study for determining pc′:(6)pc′=prexpeΓ0−eΓ0wpawp+b−e/λ,where pr = 1 kPa is the unit pressure, λ is the gradient of the CSL in the e−lnp′ plane, e is the current void ratio, wp=∫σij′dεij is the plastic energy density, σij′ and εij are the stress and strain tensors, respectively, and eΓ0, a, and b are material constants. Through equation (6), the effect of particle breakage on the stress-strain behaviour of coarse aggregates can be thus considered.


## 2.3. Hardening Modulus

The plastic modulus,H, can be obtained by using the consistency condition together with the evolving law, i.e., Δh=Lh¯, for the hardening parameter (h). Here, Δ indicates increment. Hence,(7)H=−∂f∂hh¯,where L=mTΔσ′/H and h¯ can be defined by the law for deviatoric hardening. During loading, f would evolve until reaching the critical state, where the evolution should be bounded by the peak stress surface. Therefore, one can write(8)h¯=χGMp−ηq,where G=G02.97−e2/1+ep′pa is the shear modulus [14]. G0 is a material constant. pa=100kPa is the atmospheric pressure. χ is a model parameter and can be a function of the void ratio, e.g., χ=χ11−e, for better model performance. χ1 is a model constant. Mp=M1+βpc/p′−1 is the peak stress ratio [19], in which β is a material constant. Further substituting equations (1) and (8) into equation (7), H can be obtained as(9)H=χGMpη−1.In addition to the plastic strain, coarse aggregates may also experience elastic strain during loading, which can be described via Hooke’s elastic law. Accordingly, the elastoplastic relation for coarse aggregates can be provided as(10)Δεv=Δp′K+nvHΔp′+nsHΔqmv,Δεs=Δq3G+nvHΔp′+nsHΔqms,where the bulk modulus K=2+2νG/3/1−2ν and ν is Poisson’s ratio.


## 3. Parameter Identification and Model Validation

### 3.1. Parameter Identification

There are ten parameters (G0, ν, M, λ, eΓ0, a, b, k, β, and χ1) in the extended fractional plasticity model, which can be determined via triaxial tests. G0 and ν can be calculated from the εs−q curve at the initial loading stage, where G0≈1+eG/2.97−e2/p′pa0.5 and ν≈3εs−2εv/6εs+2εv. M is determined by measuring the gradient of the CSL in the p′−q plane. Another four critical state parameters, λ, eΓ0, a, and b, can be determined by fitting the critical state points in the e−lnp′ plane, as comprehensively discussed in Liu and Zou [11]. k can be determined from the phase transformation state of the material (dg = 0), where(11)k=p0t′p0t′−pt′ln2M2pct′−2ηt2pt′2M2pct′−M2pt′−ηt2pt′,where et, pt′, pct′, p0t′, and ηt are the void ratio, mean effective stress, critical mean stress, intercept of f with p′-axis, and stress ratio, respectively, at the phase transformation state. Moreover, β can be determined from the peak stress point of coarse aggregates under drained loading (H = 0), such that(12)β=ηp/M−1pcp′/pp′−1,where ep, pp′, pcp′, and ηp are the void ratio, mean effective stress, critical mean stress, and stress ratio, respectively, at the peak state. In addition, χ1 can be obtained by fitting the ε1−q relationship of samples with different e0 [14]. The values of each model parameter are listed in Table 1.Table 1
Model parameters.
MaterialG0νMλeΓ0abβkχ1Shah Nehar rockfill [2]200.351.40.040.56601200.050.086Purulia dam rockfill [2]400.351.460.030.5312750.20.051.5Railroad ballast [5]800.201.90.232.11060000.120.21.2Limestone ballast [20]900.251.750.0781.181504500.230.052Crushed basalt [21]730.252.040.191.831050000.450.11.8


### 3.2. Model Validation

To validate the extend fractional plasticity approach, a series of triaxial test results of different coarse aggregates are simulated and discussed. Details of the material properties and test setup can be found in the relevant literatures [2, 5, 20, 21] and thus are not introduced here for the sake of simplicity. The model simulations of the stress-strain behaviour of coarse aggregates, e.g., rockfill [2] and ballast [5, 20, 21], are represented in Figures 1–6.Figure 1
Model simulations of the test results of Shah Nehar rockfill [2].
(a)(b)Figure 2
Model simulations of the test results of Purulia dam rockfill [2].
(a)(b)Figure 3
Model simulations of the test results of dense ballast [5].
(a)(b)Figure 4
Model simulations of the test results of medium dense ballast [5].
(a)(b)Figure 5
Model simulations of the test results of limestone ballast [20].
(a)(b)Figure 6
Model simulations of the test results of crushed basalt [21].
(a)(b)Figure1 shows the model predictions of the test results of Shah Nehar rockfill [2]. e0 of 0.27 with σ3′ of 200–800 kPa was used for model implementation. It can be found that the extended fractional plasticity model can well describe the typical stress-strain response, including volumetric contraction/dilation and strain hardening/softening, of Shah Nehar rockfill. Figure 2 presents the model simulations of the stress-strain behaviour of Purulia dam rockfill [2]. It was reported that the material mainly consisted of angular/subangular metamorphic rock particles. e0 of 0.252 with σ3′ of 600–1200 kPa was used for model implementation. It can be also found from Figure 2 that the model can reasonably capture the stress-strain behaviour of Purulia dam rockfill under different loading pressures.Figures3–4 simulate the test results of ballast with two different relative densities [5]. e0 = 0.620–0.686 (dense ballast) and e0 = 0.711–0.735 (medium dense ballast) with σ3′ of 100–400 kPa were used for model simulation. It is found that the extended model can provide reasonable predictions of the stress-strain behaviour of ballast with different relative densities. However, a slight underestimation of the peak stress under low σ3′ is observed. This can be attributed to the low dilatancy predicted by the developed model. A higher extent of dilatancy would result in better model simulations. Figure 5 shows the model simulations of the triaxial behaviour of limestone ballast [20]. e0 = 0.721–0.683 with σ3′ of 10–60 kPa were used for model simulation. An overestimation of the volumetric dilation of limestone ballast can be found in Figure 5. However, the overall predictions of the stress-strain relations are satisfactory. Figure 6 shows the model predictions of stress-strain behaviour of crushed basalt [21]. e0 = 0.627–0.696 with σ3′ of 50–300 kPa was used for model predictions. It can be found that the proposed model simulates well the relationship between volumetric strain and shear strain for crushed basalt. An underestimation of the peak stress under low σ3′ is also observed, but the overall evolution trend of the deviator stress with the shear strain is captured.


## 3.1. Parameter Identification

There are ten parameters (G0, ν, M, λ, eΓ0, a, b, k, β, and χ1) in the extended fractional plasticity model, which can be determined via triaxial tests. G0 and ν can be calculated from the εs−q curve at the initial loading stage, where G0≈1+eG/2.97−e2/p′pa0.5 and ν≈3εs−2εv/6εs+2εv. M is determined by measuring the gradient of the CSL in the p′−q plane. Another four critical state parameters, λ, eΓ0, a, and b, can be determined by fitting the critical state points in the e−lnp′ plane, as comprehensively discussed in Liu and Zou [11]. k can be determined from the phase transformation state of the material (dg = 0), where(11)k=p0t′p0t′−pt′ln2M2pct′−2ηt2pt′2M2pct′−M2pt′−ηt2pt′,where et, pt′, pct′, p0t′, and ηt are the void ratio, mean effective stress, critical mean stress, intercept of f with p′-axis, and stress ratio, respectively, at the phase transformation state. Moreover, β can be determined from the peak stress point of coarse aggregates under drained loading (H = 0), such that(12)β=ηp/M−1pcp′/pp′−1,where ep, pp′, pcp′, and ηp are the void ratio, mean effective stress, critical mean stress, and stress ratio, respectively, at the peak state. In addition, χ1 can be obtained by fitting the ε1−q relationship of samples with different e0 [14]. The values of each model parameter are listed in Table 1.Table 1
Model parameters.
MaterialG0νMλeΓ0abβkχ1Shah Nehar rockfill [2]200.351.40.040.56601200.050.086Purulia dam rockfill [2]400.351.460.030.5312750.20.051.5Railroad ballast [5]800.201.90.232.11060000.120.21.2Limestone ballast [20]900.251.750.0781.181504500.230.052Crushed basalt [21]730.252.040.191.831050000.450.11.8


## 3.2. Model Validation

To validate the extend fractional plasticity approach, a series of triaxial test results of different coarse aggregates are simulated and discussed. Details of the material properties and test setup can be found in the relevant literatures [2, 5, 20, 21] and thus are not introduced here for the sake of simplicity. The model simulations of the stress-strain behaviour of coarse aggregates, e.g., rockfill [2] and ballast [5, 20, 21], are represented in Figures 1–6.Figure 1
Model simulations of the test results of Shah Nehar rockfill [2].
(a)(b)Figure 2
Model simulations of the test results of Purulia dam rockfill [2].
(a)(b)Figure 3
Model simulations of the test results of dense ballast [5].
(a)(b)Figure 4
Model simulations of the test results of medium dense ballast [5].
(a)(b)Figure 5
Model simulations of the test results of limestone ballast [20].
(a)(b)Figure 6
Model simulations of the test results of crushed basalt [21].
(a)(b)Figure1 shows the model predictions of the test results of Shah Nehar rockfill [2]. e0 of 0.27 with σ3′ of 200–800 kPa was used for model implementation. It can be found that the extended fractional plasticity model can well describe the typical stress-strain response, including volumetric contraction/dilation and strain hardening/softening, of Shah Nehar rockfill. Figure 2 presents the model simulations of the stress-strain behaviour of Purulia dam rockfill [2]. It was reported that the material mainly consisted of angular/subangular metamorphic rock particles. e0 of 0.252 with σ3′ of 600–1200 kPa was used for model implementation. It can be also found from Figure 2 that the model can reasonably capture the stress-strain behaviour of Purulia dam rockfill under different loading pressures.Figures3–4 simulate the test results of ballast with two different relative densities [5]. e0 = 0.620–0.686 (dense ballast) and e0 = 0.711–0.735 (medium dense ballast) with σ3′ of 100–400 kPa were used for model simulation. It is found that the extended model can provide reasonable predictions of the stress-strain behaviour of ballast with different relative densities. However, a slight underestimation of the peak stress under low σ3′ is observed. This can be attributed to the low dilatancy predicted by the developed model. A higher extent of dilatancy would result in better model simulations. Figure 5 shows the model simulations of the triaxial behaviour of limestone ballast [20]. e0 = 0.721–0.683 with σ3′ of 10–60 kPa were used for model simulation. An overestimation of the volumetric dilation of limestone ballast can be found in Figure 5. However, the overall predictions of the stress-strain relations are satisfactory. Figure 6 shows the model predictions of stress-strain behaviour of crushed basalt [21]. e0 = 0.627–0.696 with σ3′ of 50–300 kPa was used for model predictions. It can be found that the proposed model simulates well the relationship between volumetric strain and shear strain for crushed basalt. An underestimation of the peak stress under low σ3′ is also observed, but the overall evolution trend of the deviator stress with the shear strain is captured.


## 4. Conclusions

This study enriched the original fractional plasticity by incorporating particle breakage. Several main findings are summarised as follows:(1)
A state-dependent fractional plasticity was developed by considering the effect of particle breakage on the critical state behaviour of coarse aggregates. The loading rule and plastic flow rule were obtained by conducting first-order derivative and fractional-order derivative at the yielding surface, respectively, while the state-dependent hardening modulus was developed from consistency condition.(2)
The proposed model had ten parameters which can be determined from laboratory tests. Through simulating a series of test results of different coarse aggregates, it was found that the model can reasonably capture the stress-strain behaviour of coarse aggregates with particle breakage.



---

*Source: 6642461-2021-03-22.xml*