Discrete Element Study of Particle Size Distribution Shape Governing Critical State Behavior of Granular Material

Abstract

Granular soil is a porous medium composed of particles with different sizes and self-similar structures, exhibiting fractal characteristics. It is well established that variations in these fractal properties, such as particle size distribution (PSD), significantly influence the mechanical behavior of the soil. In this paper, a three-dimensional (3D) Discrete Element Method (DEM) is applied to study the mechanical and critical-state behavior of the idealized granular assemblages, in which various PSD shape parameters are considered, including the coefficient of uniformity (C_u), the coefficient of curvature (C_c), and the coefficient of size span (C_s). In addition, the same PSDs but with different mean particle sizes (D₅₀) are also employed in the numerical simulations to examine the particle size effect on the mechanical behavior of the granular media. Numerical triaxial tests are carried out by imposing axial compression under constant mean effective pressure conditions. A unique critical-state stress ratio in p′-q space is observed, indicating that the critical friction angle is independent of the shape of the PSD. However, in the e-p′ plane, the critical state line (CSL) shifts downward and rotates counterclockwise, as the grading becomes more widely distributed, i.e., the increasing coefficient of span (C_s). Additionally, a decrease in the coefficient of curvature (C_c) would also move the CSL downward but with negligible rotation. However, it is found that the variations in the mean particle size (D₅₀) and coefficient of uniformity (C_u) do not affect the position of the CSL in the e-p′ plane. The numerical findings may shed some light on the development of constitutive models of sand that undergo variations in the grading due to crushing and erosion, and address fractal problems related to micro-mechanics in soils.

1. Introduction

The particle size distribution (PSD) refers to different particle size fractions within granular material. It is well established that these fractions are fundamental in determining the texture of the material and are a key characteristic of granular soils. The distribution of these fractions significantly influences fundamental properties such as compaction, permeability, shear strength, and consolidation. Therefore, a comprehensive understanding of particle size distribution is essential for accurately predicting the behavior of granular soils in engineering applications. Recently, research on the effects of PSD has gained traction, particularly in the areas of particle breakage [1,2,3] and gap-graded assemblies [4,5].

Fractal theory has been widely used in describing and studying the particle size distribution (PSD) of granular material [6,7,8,9,10]. However, in geotechnical engineering, the wide particle size range of granular materials, such as sand and gravel, makes fractal representation challenging. Therefore, a more in-depth study of the influence of the PSD shape is necessary. Among the various PSD descriptors, the coefficient of uniformity (C_u) has been the most widely used to characterize the PSD shape [11,12,13,14,15]. However, since this parameter is relatively simple, it may not fully capture the complexities of PSD characteristics. As a result, other parameters, such as the coefficient of curvature (C_c) [16,17,18] and the coefficient of size span C_s [19], have started to receive more attention. Note, however, that these parameters are defined as a function of individual diameters (e.g., D₁₀, D₃₀, D₅₀, and D₆₀); hence, they may not always describe fully a given PSD. A more detailed investigation of these additional PSD shape descriptors is nevertheless warranted.

Critical state theory (CST), first formulated by Roscoe et al. [20], serves as the foundation for soil constitutive models. It describes an asymptotic state in which the ratio of the shear stress q to mean effective stress p′ and the volume reach a steady condition during further shearing. Numerous studies on particle breakage, e.g., [21,22,23], have shown that the critical state locus shifts as the PSD changes. Although the shear strength at the critical state appears to be unaffected by the coefficient of uniformity C_u [24,25], Li et al. [26] found that, for a broader range of C_u, deviations in the shear strength can occur. Furthermore, as the PSD widens, the critical state locus in the e-p′ plane initially shifts downward and may later move back upward [27,28,29]. This raises the following question: how does the PSD influence the critical state shear strength and the location of the critical state line (CSL) in the e-p′ plane?

While most studies have focused on the effect of C_u on the mechanical behavior of granular materials, other PSD shape descriptors have received limited attention. The objective of this study is to analyze the impact of four PSD descriptors—mean particle size D₅₀, coefficient of uniformity C_u, coefficient of curvature C_c, and coefficient of size span C_s—on the shear and critical state behaviors of granular materials_. These descriptors are defined in Equations (1)–(3):

$${\mathrm{C}}_{\mathrm{u}}={\displaystyle \frac{{\mathrm{D}}_{60}}{{\mathrm{D}}_{10}}}$$

(1)

$$C_c={\displaystyle \frac{{D_{30}}^2}{D_{10}{D}_{60}}}$$

(2)

$$C_s={\displaystyle \frac{D_{100}-D_0}{D_{100}+D_0}}$$

(3)

where D₁₀₀ is the maximum particle diameter, D₀ is the minimum particle diameter, D₁₀ is the 10% cumulative passing particle diameter, and D₆₀ is the 60% cumulative passing particle diameter.

To investigate these effects, the Discrete Element Method (DEM) [30], a powerful and efficient tool for studying granular materials, is employed. The DEM allows for isolating the effects of the PSD while minimizing the influence of other variables, such as the particle shape. Seven distinct PSDs are carefully constructed to isolate the influence of each shape descriptor. The results show that, although the evolution of deviatoric stress varies with different PSDs, the ultimate shear strength remains invariant across all PSD descriptors. Moreover, an increase in the size span C_s and a decrease in the coefficient of curvature C_c tend to shift the CSL downward in the e-logp′ plane. However, the position of the CSL in the e-logp′ plane is not significantly influenced by PSDs with particle size scaling alone or PSDs with varying C_u.

2. Critical State Framework and Simulation Details

2.1. Critical State Framework

The critical state theory (CST) proposed by Roscoe et al. [20] is a milestone in soil mechanics and has fundamentally shaped the field [31]. The critical state (CS) refers to the ultimate state for particulate materials, in which the material continues to deform in shear at a constant volume and under constant stresses. Mathematically, the critical state can be expressed as follows:

$$\dot{p}=0,\dot{q}=0,{\dot{\epsilon }}_v=0 \mathrm{but} {\dot{\epsilon }}_q\ne 0$$

(4)

where p is the mean effective normal stress, q is the deviatoric stress, ε_v is the volumetric strain, and ε_q is the deviatoric strain. The over-dot notation indicates the rate of change in each quantity.

Two sufficient and necessary conditions for the critical state were defined as follows:

$$\eta ={\eta }_c{(q/p)}_c=M, e=e_c={\widehat{e}}_c\left(p\right)$$

(5)

where η is the stress ratio, M is the critical state stress ratio, and e_c is the critical state void ratio, which is only related to p, resulting in a unique critical state line (CSL) in the e- space. The first part of Equation (5) defines a straight line in the p-q space, which is termed the critical state line in the p-q space. The slope M of this line is an intrinsic material constant, which is also a function of the stress Lode angle.

Although both conditions in Equation (5) must be simultaneously satisfied for critical state failure, they can be independently reached. Li and Dafalias [32] revisited this theory from the thermodynamics perspective and proposed an anisotropic critical state theory (ACST) to address the limitations of the CST.

In the context of a discrete system with N_f particles, the fabric of the system can be defined by the deviatoric tensor F, as follows [33,34]:

$$\mathbf{F}={\displaystyle \frac{1}{N_f}}\sum _{\alpha =1}^{N_f}{n}^{\alpha }\otimes n^{\alpha }-{\displaystyle \frac{1}{3}}\mathbf{I}$$

(6)

where I is the identity tensor, and $n^{\alpha }$ represents a characteristic orientation of the discrete system, which could include void vectors [35], particle orientations [36], or contact unit normals [37]. As further explained by Li and Dafalias [38], the fabric tensor must be thermodynamically consistent with the dissipation mechanism and thus should be a per-unit volume measure. To account for the relative orientation between the fabric anisotropy and loading directions, a fabric anisotropy variable (FVA), denoted as A, is introduced:

$$A=\mathbf{F}:\mathbf{n}=F{\mathbf{n}}_F:\mathit{n}=FN$$

(7)

where n is the deviatoric unit loading direction. The fabric tensor F has two non-trivial invariants: the norm F and the Lode angle θ_F, associated with the unit direction of F, denoted as n_F. The Lode angle θ_F is 0° for triaxial compression and 60° for triaxial extension. The scalar N represents the product of n_F and n, signifying the relative orientation between the loading and fabric directions. The additional condition of the fabric anisotropy at the critical state can be expressed as $A=A_c=1$. This condition implies that, during plastic deformation, the fabric evolves such that its norm gradually approaches a critical value and its direction aligns with the loading direction. The necessary and sufficient conditions for the critical state can therefore be written as

$$\eta ={\eta }_c{(q/p)}_c=M, e=e_c={\widehat{e}}_c\left(p\right),A=A_c=1$$

(8)

In this paper, the effect of four PSD shape descriptors is investigated under the anisotropic critical state framework. The contact-normal-based fabric tensor, which has been widely shown to efficiently reveal insights into sand behavior [39,40,41], is adopted to describe the internal fabric in this study.

2.2. Simulation Details

The open-source program YADE, developed by Kozicki and Donzé [42], is used to perform 3D Discrete Element Method (DEM) simulations in this study. The periodic boundary condition is applied to a cubic packing composed of unbreakable spheres. The cubic specimens, consisting of over 5000 particles, are considered as a representative volume element. All the particles within the cell move according to the specified strain increment Δε_ij, as shown in Equation (9) [43]:

$$\Delta x_i=\Delta {\epsilon }_{ij}{x}_j\Delta t$$

(9)

where Δx_i is the displacement of particle i, and Δt is the time step.

During the consolidation (isotropic compression) and shearing phases, an unbalanced force ratio below 0.001 was defined to ensure initial equilibrium after consolidation and quasi-static conditions during shear. Gravity is set to zero throughout all phases of the simulation. In the quasi-static state with zero gravity, increasing the particle density accelerates the simulation without affecting the accuracy of the results. To achieve this, all particles are assigned a solid density of ρ = 2650 × 10³ kg/m³, which is 1000 times the usual density. Additionally, a default numerical damping coefficient ξ = 0.7 is introduced to dissipate the kinetic energy of the particles, thereby accelerating the simulation and reducing the computational costs. A parallel sensitivity study indicated that using lower values of damping did not affect the simulation results.

A linear contact law is employed, and the particle’s normal stiffness is set to equal to the tangential stiffness, k = 10⁸ N/m. Initially, spheres with the desired number of particles are generated without contact. The cubic sample is then compressed isotropically to the target pressure of 100 kPa, 200 kPa, 500 kPa, 1000 kPa, and 1500 kPa. Before shearing, the coefficient of friction between particles is set to 0.5. Subsequently, the samples are sheared to larger axial strains, at which point, the critical state can be reached.

The servo mechanism for strain rate loading, developed on the basis of the method of Thornton and Zhang [44], follows the general form of

$$\dot{\epsilon }=G\left({\sigma }_t-{\sigma }_c\right)$$

(10)

where G is the gain parameter, determined through trial and error, σ_t is the target stress, and σ_c is the current stress. The mean effective stress is kept constant during the triaxial shearing process.

In this study, four PSD descriptors (D₅₀, C_u, C_c, and C_s,) are investigated. The seven curves are intentionally designed to isolate the effect of each descriptor, as shown in Figure 1. On the basis of the definitions of these descriptors, it is clear that they are related to particle sizes D₀, D₁₀, D₃₀, D₅₀, D₆₀, and D₁₀₀. Therefore, in the PSD design, changes in these specific particle sizes were considered. Additionally, D₉₀ is also included in the analysis. To begin with, four PSDs (PSD1 to PSD4) are designed, all having the same D₅₀ but differing in their maximum and minimum particle sizes. Next, PSD5 is selected through trial and error, using the above nine particle sizes corresponding to different cumulative percentages, ensuring that PSD5 differs from PSD2 only by C_s and from PSD3 only by C_u. For the effect of D₅₀, a mono-sized PSD is chosen as a control. To further investigate this effect, another mono-sized PSD with a larger D₅₀ is generated and named PSD6. A similar method is used to generate PSD7 based on PSD4 in order to study the role of C_c. This design approach allows for the separate evaluation of the influences of C_s, C_u, C_c, and D₅₀ on the behavior of the granular material. The details of the PSD descriptors are summarized in

Table 1. Descriptors of particle size distribution.

PSD Type	D50 (mm)	Cu	Cc	Cs
PSD1	0.26	1.000	1.000	0
PSD2	0.26	1.914	1.013	0.545
PSD3	0.26	2.564	1.002	0.765
PSD4	0.26	9.595	1.007	0.965
PSD5	0.26	1.914	1.013	0.765
PSD6	0.52	1.000	1.000	0
PSD7	0.26	9.595	0.761	0.965

.

3. Results

3.1. Typical Macroscopic Responses

The shear strength evolution of a group of samples with seven different PSDs (with p₀′ = 500 kPa) under a constant mean effective stress path is presented in Figure 2. It is evident that the pace of deviatoric stress evolution varies across the different PSDs during the early stages of shearing (ε_a ≤ 10%). However, with further shearing, all specimens eventually reach the same shear strength.

To better understand the impact of the various PSD descriptors on the deviatoric stress evolution, these curves are grouped into four categories in Figure 3. From these plots, it is clear that the coefficient of uniformity (C_u) has almost no effect on the deviatoric stress evolution, as the curves for C_u (represented by the yellow and blue lines in Figure 3b) almost overlap from the initial shearing state to the critical state.

On the other hand, the mean particle size (D₅₀) appears to have a slight effect on the evolution rate of the deviatoric stress, as shown in Figure 3a. However, the overall impact of these two descriptors (C_u and D₅₀) is relatively limited.

In contrast, the deviatoric stress evolution rate is significantly influenced by the other two descriptors: the coefficient of curvature (C_c) and the coefficient of size span (C_s). Specifically, samples with a larger C_c and smaller C_s exhibit higher stiffness and higher peak deviatoric stresses. This indicates that both C_c and C_s play a more pronounced role in controlling the shear strength behavior compared to C_u and D₅₀.

The volumetric strain responses are shown in Figure 4. The samples with PSD4 and PSD2 exhibit slightly dilative behavior before reaching the critical state, consistent with previous observations of strain-softening responses. In contrast, the other five samples display volumetric contraction.

To better understand the impact of different PSD descriptors on volume evolution, the results are grouped into four categories in Figure 5. For each group, despite having the same initial void ratios, the volumetric response influenced by the four PSD shape descriptors varies. Figure 5a demonstrates that all samples exhibit the same tendency of volumetric contraction until the critical state is reached, with negligible differences in volumetric strain. This suggests that the mean particle size (D₅₀) has little effect on both the volume evolution and the critical state void ratio.

Figure 5b shows that the coefficient of uniformity (C_u) has no influence on the critical state void ratio or the volumetric evolution throughout the shearing process. In contrast, the coefficient of curvature (C_c) and the coefficient of size span (C_s) significantly affect the volumetric response, as shown in Figure 5c,d. Specimens with PSD4 and PSD2 exhibit dilative behavior during shearing, while specimens with PSD7 and PSD5 experience volumetric contraction.

3.2. Typical Microscopic Responses

The coordination number (Z_m) serves as a scalar measure of the contact features in granular materials. Defined as the ratio of twice the number of contacts to the number of particles, it quantifies the average number of contacts per particle and provides critical insights into the microstructural evolution of the material. The mathematical expression for Z_m is as follows:

$$Z_m=2{\displaystyle \frac{N_c-N_p^1}{N_p-(N_p^1+N_p^0)}}$$

(11)

Here, N_p and N_c represent the total numbers of particles and contacts, respectively. $N_p^1$ and $N_p^0$ are the numbers of particles with only one contact and no contacts, respectively. This definition excludes particles with fewer than two contacts, as they do not contribute to stress transmission.

Figure 6 shows the evolution of Z_m with axial strain. Consistent with findings from previous DEM studies, Z_m decreases with increasing axial strain. The influence of PSD descriptors on the Z_m evolution is illustrated in Figure 7. Similarly to the macroscopic responses, D₅₀ and C_u exhibit a negligible influence on the coordination number at the critical state, as demonstrated in Figure 7a,b.

In contrast, Figure 7c,d reveal that the coordination number is significantly influenced by C_c and C_s. The samples with higher C_c values tend to have larger coordination numbers, indicating a more interconnected contact network. Conversely, the samples with lower C_s values show a similar trend, reinforcing the role of size span in influencing particle contacts. This finding suggests that a wider size span may increase the proportion of floating particles—those with insufficient contacts to contribute to stress transmission—thereby reducing the overall coordination number.

In addition to the coordination number, the internal structure, or fabric, of granular soils is also investigated. The fabric of granular materials refers to the collective microstructural features within the Representative Elemental Volume (REV), including the spatial arrangement of grains, the distribution of voids, and the interactions between particles [45].

In this study, a contact normal vector-based deviatoric second-order fabric tensor F is used to characterize the microstructure of the granular assemblies. This approach is derived from Kanatani’s equation [46] where f(n) is the density function that defines the distribution of contact normals. The fabric tensor G is the symmetric second-order tensor, and F represents its deviatoric part. The mathematical formulation is as follows:

$$G_{ij}=\underset{n}{\int }f\left(\mathit{n}\right)n_i{n}_{j}dn$$

(12)

$$f\left(n\right)={\displaystyle \frac{1}{4\pi }}\left(1+F_{ij}{n}_i{n}_j\right)$$

(13)

$$F_{ij}={\displaystyle \frac{15}{2}}\left(G_{ij}-{\displaystyle \frac{1}{3}}{\delta }_{ij}\right)$$

(14)

$$F=\sqrt{F_{ij}{:F}_{ij}}$$

(15)

Following Yang and Wu [47], this deviatoric fabric tensor is normalized by dividing by the volume, such that $F^{\prime }={\displaystyle \frac{F}{v}}={\displaystyle \frac{F}{1+e}}$, where e is the void ratio. As previously mentioned, F has two invariants: the norm F and unit direction n_F. The direction n_F is related to the Lode angle of the fabric, θ_F, through the expression $\cos 3{\theta }_F=\sqrt{6}\mathrm{tr}{n}_F^3$.

The evolution of the fabric tensor during shearing was examined for all the PSDs. Figure 8 shows that the fabric norm increases with the axial strain for all samples, regardless of the PSD. Although the rate of fabric evolution varies, all the samples eventually reach a similar fabric value at the critical state. Regarding the fabric direction, a significant reduction in the fabric Lode angles is observed at the very beginning of shearing (less than 2–3% axial strain), which occurs much more rapidly than the evolution of the fabric norm (which takes approximately 30–50% of the axial strain). Figure 9 demonstrates that the fabric Lode angles eventually drop to zero, accompanied by slight fluctuations, indicating that the fabric direction quickly aligns with the loading direction.

3.3. Critical State Responses

As mentioned before, three planes are usually studied under anisotropic critical state theory, which are the p′-q plane, e-p′ plane, and fabric space. In this study, all the critical stress states obtained for the different PSDs are shown in the p′-q plane in Figure 10. It can be concluded that all samples exhibit a unique critical stress ratio M_c ≈ 0.77, which is consistent with previous Discrete Element Method (DEM) results [48,49,50]. The data from all samples can be fitted by a single straight line, suggesting that none of the PSD shape descriptors are correlated with the critical state strength.

Figure 11a illustrates the relationship between the critical state void ratio and mean effective stress in a logarithmic coordinate system. It can be observed that the location of the critical state line (CSL) in the e-logp′ plane varies with the PSD shape. Specifically, both higher C_s and lower C_c lead to a lower CSL position. In general, it can be stated that a wider PSD results in a lower CSL location. Interestingly, although the PSDs of samples 3 and 5 differ due to variations in C_u, their CSLs coincide. Similarly, the CSLs of parallel PSDs (PSD1 and PSD6) appear to overlap.

To investigate this further, the critical state lines for non-cohesive soils are straightened in the e⁻ − (p′/p_a)^ξ plane by use of the approach by Li and Wang’s approach [51], as expressed in Equation (16):

$$e_c=e_{\Gamma }-{\lambda }_c{\left(p\prime /p_a\right)}^{\xi }$$

(16)

where p_a is the atmospheric pressure (101.325 kPa), $e_{\Gamma }$ is the intersection of the CSL on the e-axis, ${\lambda }_c$ is the slope of the CSL, and ξ is a constant parameter used for fine-tuning. As noted by [52], ξ is not a sensitive parameter and is set to 1.0 for simplicity.

Figure 11b shows the results of the seven PSDs grouped into five categories. In general, a wider PSD (with larger C_u or C_s) corresponds to a lower CSL position. The slopes of the CSLs also vary within a narrow range. Thus, the positional changes of the CSL due to PSD variations can be described by two main components: a downward translational shift and a counterclockwise rotational shift. These two shifting modes are consistent with previous experimental findings. It appears that C_c moves the CSL downward with negligible rotation, while C_s causes both a downward movement and a counterclockwise rotation.

Figure 12 further illustrates the effect of PSD descriptors on the critical state void ratio, e_c. It is evident that D₅₀ and C_u have no significant effect on the critical state void ratio. However, Figure 12c,d show that the CSL for samples with higher C_c and lower C_s is positioned higher in the e-p′ plane. Interestingly, the intersection of the CSL seems to be related to the coordination number. As shown in Figure 13, as the $e_{\Gamma }$ increases, the coordination number at the critical state also increases.

Figure 14 presents the critical values of F′_c from simulations conducted with seven different PSDs. The coefficient of variation (CoV) is used to quantify the degree of deviation of the fabric norm from its mean value. Overall, the critical fabric norm is found to be a unique value, F′_c = 0.293. This suggests that the critical state fabric norm is independent of both the critical void ratio e_c and the mean effective stress p′, which is consistent with the anisotropic critical state theory.

4. Discussion and Conclusions

Although the effect of the PSD shape on the critical state behavior of granular materials has been widely studied, particularly concerning the coefficient of uniformity (C_u), the relationship between the PSD shape and the location of the critical state line (CSL) remains unclear. This paper expands on previous studies by considering four PSD shape descriptors—the mean particle size (D₅₀), coefficient of uniformity (C_u), coefficient of uniformity (C_c), and coefficient of uniformity (C_s). A series of systematic 3D-DEM simulations was conducted to analyze the effects of these different PSD shape descriptors on the mechanical behavior of granular materials. Seven distinct PSDs were carefully designed so that the influence of each descriptor could be investigated separately. Particular emphasis was placed on critical state behavior, interpreted within the framework of anisotropic critical state theory (ACST). A series of constant-mean effective stress drained triaxial tests were performed, and the critical state behavior was examined on both the macroscopic and microscopic scales in terms of the critical stress ratio, critical state void ratio, and critical state fabric norm. According to the simulation results, the following conclusions can be drawn:

• Isolation of the effect of PSD shape descriptors:

The DEM numerical strategy employed in this study allows for the isolation of the effects of PSD shape descriptors on the mechanical behavior of a granular assembly. It was found that not all descriptors significantly impact the mechanical responses. Specifically, D₅₀ and C_u had little to no effect on most of the mechanical properties, including the deviatoric stress, volumetric strain, and coordination number. Note that the simulations shown here are performed by use of periodic boundaries (free from boundary effects); hence, the effect of D₅₀ is not expected; however, it is clear from the results presented that the increased particle masses resulting from scaling up the PSD have a slight effect. However, both C_c, and C_s had a significant impact on the macroscopic mechanical properties, microscopic responses, and critical state properties. This also emphasizes the need to accurately describe the PSD. In some cases, it may be necessary to use alternative PSD descriptors (such as when gap-graded soils or soils with a significant amount of fines are considered).

2.

Critical state strength and void ratio:

The critical state strength appears to be insensitive to the PSD shape descriptors considered in this study. In contrast, the critical state void ratio was significantly affected by the PSD shape. Generally, a wider PSD (with higher C_u and C_s values) results in a lower CSL position in the e-p′ plane. While both C_u and C_s represent the range of particle sizes, C_s significantly affects the CSL position, while C_u has a minimal effect.

3.

CSL Shifting Modes:

Two distinct shifting modes of the CSL in the e-p′ plane were observed: downward translational shifts and counterclockwise rotational shifts. An increase in C_s caused both a downward movement and counterclockwise rotation of the CSL. In contrast, a decrease in C_c resulted in a downward shift with negligible rotation.

4.

Coordination Number and CSL Intersection:

The intersection of CSL appears to be related to the coordination number. As the e_Г increases, the coordination number at the critical state also increases. Moreover, the coordination number stabilizes at relatively small strain levels (<10%), while the fabric continues to evolve towards the critical state.

5.

Microstructure Evolution:

The evolution of the microstructure is quantified using the deviatoric fabric tensor, which captures both the direction and norm of the fabric. The fabric direction evolves towards the loading direction at the early stages of shearing. The fabric norm evolves more slowly and appears to reach a unique value for all seven PSDs. The critical state fabric norm is independent of both the critical void ratio and the mean effective stress. Furthermore, while the critical fabric norm is generally consistent across materials, samples with different PSDs tend to converge on the same value for the critical state fabric norm.

This study demonstrates that the PSD shape parameters have distinct effects in the two critical state spaces. On the one hand, the stress ratio at the critical state is shown to be independent of the PSD. While previous research has suggested similar findings, it has primarily considered C_u alone. This study extends the validity of these results by showing that the critical state shear strength is independent not only of C_u, but also of D₅₀, C_u, and C_c. Thus, we can confidently state that there is no relationship between the PSD and critical state shear strength.

On the other hand, the simulation results highlight that the critical state void ratio is significantly affected by the PSD shape. Notably, C_u has no impact on the dilatancy or critical state void ratio, which contrasts with some previous studies. This discrepancy may arise from those studies considering only C_u and not controlling for other PSD shape descriptors, such as C_c and C_s, which were found to be the main factors influencing both the macroscopic and microscopic responses in this study. It appears that the pattern of these responses cannot be captured by any single descriptor. Therefore, a more comprehensive description of the PSD needs to be developed.

These numerical findings may offer valuable insights for developing constitutive models of sand that undergo grading variations, such as those caused by crushing and erosion, and for addressing fractal problems related to micro-mechanics in soils.