The influence of particle geometry and the intermediate stress ratio on the shear behavior of granular materials

The behavior of granular materials is very complex in nature and depends on particle shape, stress path, fabric, density, particle size distribution, amongst others. This paper presents a study of the effect of particle geometry (aspect ratio) on the mechanical behaviour of granular materials using the discrete element method (DEM). This study discusses 3D DEM simulations of conventional triaxial and true triaxial tests. The numerical experiments employ samples with different particle aspect ratios and a unique particle size distribution (PSD). Test results show that both particle aspect ratio (AR) and intermediate stress ratio (b=(σ2′-σ3′)/(σ1′-σ3′))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(b=({\upsigma }_{2}'-{\upsigma }_{3}')/({\upsigma }_{1}'-{\upsigma }_{3}'))$$\end{document} affect the macro- and micro-scale responses. At the macro-scale, the shear strength decreases with an increase in both aspect ratio and intermediate stress ratio b values. At the micro-scale level, the fabric evolution is also affected by both AR and b. The results from DEM analyses qualitatively agree with available experimental data. The critical state behaviour and failure states are also discussed. It is observed that the position of the critical state loci in the compression (e-p′)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(e-p')$$\end{document} space is only slightly affected by aspect ratio (AR) while the critical stress ratio is dependent on both AR and b. It is also demonstrated that the influence of the aspect ratio and the intermediate stress can be captured by micro-scale fabric evolutions that can be well understood within the framework of existing critical state theories. It is also found that for a given stress path, a unique critical state fabric norm is dependent on the particle shape but is independent of critical state void ratio.


Introduction
Granular materials are commonly encountered in nature and industry. For example, they are found in landslides and avalanches, as raw minerals for extraction and transport, cereal storage, powder mixing, etc. Cohesive frictional materials like concrete are also made of granulates [1]. The mechanical response of the granular materials is very complex and can be affected by confining pressure, density, stress path, particle shape, particle size distribution, amongst other factors [2][3][4]. It is widely recognized that particle characteristics affect the mechanical behavior of granular soils, which has been investigated extensively by experimental work. For example, Cho et al. [5] explored the effects of particle shape on packing density and on the small-tolarge strain mechanical properties of sandy soils and found that the increasing particle irregularity causes an increase in the critical state angle of shearing resistance. Cavarretta et al. [3] and Yang and Luo [6] employed spherical glass beads and crushed angular beads within mixtures including sand to explore the relationship between the critical state behavior and particle shape. The latter experiments show that critical state behavior and liquefaction potential are affected significantly by particle shape. Granular materials experience different stress paths in nature [7]. Previous physical experiments and DEM simulations show that the intermediate stress ratio (b = (σ 2 − σ 3 )/(σ 1 − σ 3 )) affects the soil response [8][9][10]. More recently, Xiao et al. [11] performed a series of true triaxial experiments on coarse granular soil and found that b greatly affects the critical state line in both e − p and p − q space. These findings are questionable as most existing investigations unanimously confirm that the uniqueness of the critical state line (CSL) is independent of stress paths. In addition, careful scrutiny of the critical state data in the e − p plane by Xiao et al. [11] indicates that the data points are indeed very sparse and the discrepancy between different b cases for CSLs is insignificant, with the maximum void ratio difference only being less than 0.02 for p = 200 kPa. More critical state data are therefore required with a larger range of stresses and void ratios in order to draw convincing conclusions.
The discrete element method (DEM) developed by Cundall and Strack [12] has been demonstrated to be a promising tool to improve our understanding of the fundamental behavior of granular materials. A key feature of DEM simulations is that the same numerical sample can be subjected to different stress paths, therefore the sample repeatability can be ensured in the simulation, whilst also providing additional information which is not easily measured from physical experimental tests. In DEM simulations soil particles are idealized entities with spherical shape being the most common type of particle geometry considered in 3D DEM simulations. Many DEM based numerical simulations have been performed to investigate the effect of particle shape and intermediate stress ratio on the mechanical behavior of granular media. Nouguier et al. [13] performed 2D numerical simulations with samples made of irregular hexagons with different particle elongation (i.e. ratio of maximum to minimum length) and found that the angle of shearing resistance decreases as the particle elongation increases for samples with the same initial density. Zhao et al. [14] used a composite (clumped) particle model to create non-spherical particles and demonstrated that particle shape increased resistance to compaction. Kumari and Sitharam [15] studied the effect of particle aspect ratio on the monotonic shear behavior, but their study was limited to granular assemblies with only two aspect ratios and the particle grading of these DEM specimens was not unique. More importantly, their study only focused on the anisotropic responses in shear strength and the corresponding microscopic behavior, while no analyses of the quantification of critical-state fabrics were performed. Barreto and O'Sullivan [10] considered the effects of friction μ and the intermediate stress ratio b under constant mean effective stress ( p ) and constant b fixed conditions to postulate the mechanism that explains the effect of b on material response. Zhao and Guo [16] found that a unique CSL in the e-log p space, which was then confirmed by Yang and Wu [17] using anisotropic assemblies consisting of clumped particles. Huang et al. [18] used DEM to understand the effect of b on the critical state and discussed that the critical state line may not be unique, but the discrepancy between different b cases for CSLs is not very significant and it might be argued that their results are not accurate [17]. Bear in mind that the unique CSL assertion is derived from sound theoretical proof [19]. Also note that the studies by Zhao and Guo [16] and Huang et al. [18] only considered spherical particles. There is no prior research that systematically considers both particle shape and generalized stress conditions such as those measured on true triaxial tests using DEM.
This study presents a comprehensive DEM analysis of particle aspect ratio effects on the mechanical behavior of granular materials under generalized stress conditions. Shear strength, dilatancy characteristics, critical state response, failure criterion, amongst other features of granular materials' response are considered in this study. A series of DEM simulations including conventional triaxial compression tests, as well as constant p and constant b true triaxial tests are conducted in this study, in which numerical specimens with varying aspect ratio but a unique particle size distribution are considered. The main objective of this research was to investigate the critical state behavior of granular assemblies composed by grains with varying ARs, when subjected to generalized stress conditions (i.e not restricted to the triaxial space). Conventional triaxial compression tests under constant volume (undrained) and constant cell pressure (drained) were performed to assess the effect of AR on peak resistance and critical state behavior. Similarly, DEM simulations of true (drained) triaxial tests under constant p and different (constant) values of b were performed to assess the effect of both AR and b. The effect of AR might refer to different granular materials, while the effect of b refers to the effect of loading under generalized stress conditions. To differentiate these two effects, the simulation sets in which AR is kept constant were performed under different constant b values.

Particle characteristics
Particle shape in granular materials is generally complex and often characterized by elongation, roundness, angular-  ity, texture, etc. [5]. In this study, the elongation of particle characterized by the aspect ratio is considered. As illustrated in Fig. 1, the aspect ratio (AR) is defined as the ratio a 2 /a 1 between the minor and major elongation axes. For a spherical particle, AR = 1, and the higher the particle elongation the smaller AR. For the numerical DEM samples used in this study, a power law size distribution is chosen to generate particles to replicate a real sand grading and also avoid numerical difficulties [20,21]. As shown in Fig. 2, the particle size ranges vary between 0.4-1.2 mm equivalent diameters. To assess the effect of the aspect ratio of the particles, clumped particles with varying particle aspect ratio were created. To guarantee a fixed particle size distribution-PSD (Fig. 2) for all numerical specimens the volume of individual particle clumps were set at the same value as for spherical particles (AR = 1). The uniformity coefficient D 60 /D 10 of the resulting PSD is 1.33. Three different ratios were used for this study as illustrated in Table 1, which also shows their equivalent radii.

Simulation procedure
The open source DEM code YADE [22] was used in this study. A standard linear contact model was chosen. Its input parameters are listed in Table 2. A cubical assembly of particles with periodic boundaries was adopted to assess granular material behavior free from boundary effects. Ini- tially an assembly of particles with no contacts between clumps was generated and then isotropically compressed to desired mean effective stresses using a servo-control algorithm. A two-staged compression was performed. In the first stage, to obtain varying densities of the samples, different inter-particle friction coefficients were assigned onto the particles and the confining pressure was gradually increased to 180 kPa. In the second stage, the inter-particle friction coefficient was restored to a default value of μ = 0.5 and the confining pressure was then increased to the target mean stress. Static equilibrium after compression was ensured by reaching an unbalanced force ratio lower than 0.001 as recommended by Kuhn et al. [23]. After this, conventional triaxial compression under both constant cell pressure and constant volume conditions were carried out. In addition, constant p and constant b true triaxial simulations were also performed. In DEM simulations, density scaling is often used to reduce computational cost [1,4,24,25]. Thornton and Antony [26] noted that the use of particle density scaling has no effect on the contact forces and displacements in case of quasi-static conditions. In this study, particle density was scaled up by a factor of 1,000. The selection of an appropriate strain rate is essential to ensure quasi-static conditions and maximize computational efficiency. Forterre and Pouliquen [27] state that in the quasi-static regime, the time required for macroscopic deformations is relatively large compared to that required for microscopic rearrangements, resulting in inertial number I (=γ d/ √ P/ρ P ) values lower than 10 −3 , whereγ is the strain rate, d is the particle diameter, P is the relevant stress and ρ P is the particle density. In addition to ensuring that I 10 −3 for all simulations presented here, quasi-static conditions were verified by performing a preliminary parametric study using different strain rates. A strain rate of 0.05 s −1 was then selected for this study. To quantify the relative density (RD) of the samples, the maximum and minimal void ratios should be determined. According to Yang et al. [25], the maximum void ratio can be estimated by generating samples with inter-particle friction coefficient μ = 0.5, and the minimum void ratio is obtained with μ = 0.0, both under a reference state with confining pressure = 10 kPa.

Simulation results
In three dimensional analysis, the mean effective stress ( p ) and deviator stress (q) are defined as: The term 'fabric' is frequently used to describe the microstructure of granular materials (e.g. [28]). From a microstructural perspective, for an aggregate having N measurements of the quantities that reflect the microstructure of the material, such as particle orientation, orientation of unit contact vectors, branch vectors, etc., the fabric can generally be represented by the deviatoric part F of a tensor G according to where m k is a unit vector along the k th measurement of the characteristic direction, say, the direction of contact unit normal [29], and w(m k ) is a weighting factor. If w(m k ) = 1 as in the case of a contact normal based definition of F, trG = 1.
Here F will be referred to as the fabric tensor. In order to identify the norm F and direction n F of F, the following expressions can be invoked The fabric tensor F defined in this paper is a deviatoric tensor, and has two non-trivial invariants. We used the invariants norm F and the unit direction n F to represent the tensor F, in which the norm F signifies the degree of fabric anisotropy. Therefore F = 0 refers to isotropic fabric and F > 0 means the fabric is anisotropic. A scaled version of F defined in Eq. (3) is given in Eq. (5), which is widely used in the literature [16,17,21] to quantify the fabric anisotropy of soils. It is noted that the same symbol of F is used in both Eqs. (3) and (5) to avoid introducing more variables. The scaled F has been referred to and used in the data interpretation for this paper. Figure 3 shows the macro-and micro behavior of DEM specimens sheared under constant volume (undrained) conditions. All the samples had the same initial density with relatively density D r ≈ 20% and were isotropically compressed to 200 kPa, but with different particle aspect ratios. Figure 3a illustrates the stress-strain response. At large strains (ε a ≥ 40%) the deviatoric stress reached a constant value corresponding to a critical state that was dependent on the aspect ratio (AR). All the samples exhibit a hardening behavior and the maximum deviatoric stress decreases with an increase in AR. It can be seen from Fig. 3b that the effective stress path was also affected by particle shape. There was a contractive response prior to the dilative behavior with higher effective stress towards the critical state, as observed in real sands. This has been confirmed by the trends of the pseudo pore pressure for a dry packing, equivalent to the change in σ 3 , which is held constant in real undrained tests [21]. In terms of micro-scale response, Fig. 3c illustrates the evolution of the norm F of the fabric tensor F. The norm F increases with axial strain, and reaches constant values at the critical state. Note that the fabric norm F decreases as AR increases. Figure 3d illustrates the evolutions of the product of n : n F against the deviatoric strains, representing the relative orientation between fabric direction n F and stress direction n. It is seen that unlike the fabric norm, n : n F increases with the strain and stabilizes at unity, indicating that the fabric direction is co-directional with the loading direction at the critical state. Similar results are also reported by Zhao and Guo [16], Yang and Wu [17]. Figure 4 shows the results of triaxial compression under constant cell pressure (drained) conditions. In Fig. 4a the deviatoric stress reached a peak value followed by softening behavior. Once again, at large strains the deviatoric stresses became relatively constant and this can be considered as a critical state. This was confirmed by the constant volume conditions illustrated in Fig. 4b at large strains. All the samples experience contraction followed by significant dilation. The evolution of fabric shown in Fig. 4c, d is consistent with that discussed for the constant volume conditions.

True triaxial (constant b) tests
Simulations with constant p (=200 kPa) and constant b were also carried out using different values of AR. These sets of   Figure 5a shows that both the peak and critical deviatoric stresses decrease with increasing b value. This agrees with previous DEM and experimental studies [7,11,18]. The difference in initial shear stiffness appears to be insensitive to b value. Figure 5b shows the volumetric responses, evidencing a dilative behavior for all values of b. The maximum dilation decreased as b increased as reported by Sazzad et al. [7] and Huang et al. [18]. Constant deviatoric stresses and volumetric strains at large strain ε a = 30%, indicated that the critical state was reached. Figure 5c shows that the fabric norm F increases with axial strain and it was nearly constant when ε a ≥ 30%. Although there were certain fluctuations, the norm F for b = 1 was slightly higher than for b = 0, as seen in previous studies [16,17,30]. Figure 5d shows the evolutions of product of n : n F . It became constant for all the b values, indicating that the fabric direction was co-directional with the loading direction at large strains.
Previous DEM research [10,16,18] has used spherical particles to study the influence of intermediate stress ratio on material response. We extend these studies by assessing the influence of particle shape due to different aspect ratios. In geotechnical engineering, it is common to use ϕ = sin −1 [(σ 1 − σ 3 )/(σ 1 + σ 3 )] to define the shear strength of soils. The effect of AR and b on the angle of shearing resistance ϕ at the critical state is illustrated in Fig. 8 and ϕ decreases with AR for all the b values considered. There is also a maximum value around b = 0.5. The trends obtained here agree well with test data reported by Sutherland & Mesdary [31], Barreto and O'Sullivan [10], Thornton [32], amongst others. In this case, differentiation of the effects of b and AR is challenging, but there seems to be a larger effect of AR than b on ϕ for the three sets of simulations.
Previous studies [10,30,33] show that the DEM results compare well with the Lade and Duncan [34] failure criterion defined by: in which I 1 = σ 1 + σ 2 + σ 3 is the first stress invariant and I 3 = σ 1 σ 2 σ 3 is the third stress invariant. Figure 9 presents the DEM results in the deviatoric plane together with their corresponding Lade and Duncan failure criterion. The agreement with the data is excellent. Barreto and O'Sullivan [10] reported that the Lade parameter η is a function of the inter-particle friction. As shown in Fig. 10, the Lade parameter η is also a function of AR with η = 27.09 − 22.74 AR in this study. Although very sparse data points indicate a simple linear function between AR and η, a convincing relationship could be obtained if more data points were available.
The fabric tensor can be manipulated in the same way as the stress tensor. Thornton [32] proposed a parameter η * f that determines the position of the fabric envelope in the deviatoric plane and is given by:  Figure 11 shows that the parameters defined by Eq. (7) agrees well with the data sets for different aspect ratios. As shown in Fig. 12, there is also a linear relationship between η * f and AR, with η * f = 1.863 ± 0.0554 AR, which, again, might be subjected to further verification by additional data.
The observation above is qualitatively consistent with the earlier DEM simulations on spherical particle assemblages by Zhao and Guo [16], Thornton and Zhang [30], and Thornton [32]. Recently, Li and Dafalias [35] have proved from theoretical perspective that the contour of the critical state fabric norm F c on the π-plane is approximately reciprocal to the shape of critical state yield surface. In the present study, the norm ratio of the critical fabric tensor F TC c /F TE c ranging from 1.06 to 1.10 for various AR values considered in this study, which is qualitatively found to be reciprocal to the ratio of critical stress ratios between TC and TE, i.e. previous numerical results by Zhao and Guo [16] and the recent theoretical predictions by Li and Dafalias [35].

Critical state
The classical critical state theory (CST) developed by Roscoe et al. [36] and Schofield and Wroth [37] laid the foundations for critical state mechanics. The critical state is defined as the state of shear deformation without change in shear stress and volume and can be analytically expressed by: in which p is the mean normal effective stress, s the deviatoric stress tensor, ε v the volumetric strain, e the deviatoric strain tensor, and a superposed dot denotes the rate. Upon reaching the critical state, the stress ratio η = q/ p and the void ratio e will satisfy the following conditions, η = η c = q/ p c = M and e = e c =ê( p ) Li and Wang [38] proposed a linearization approach to straighten the critical state line (CSL) in the e − ( p / p a ) ξ plane for sand, which can be expressed as, where p a is the atmospheric pressure serving as a reference pressure for normalization, e is the intersection of CSL on eaxis at p = 0, λ c is the slope of the line and ξ is a parameter used for fine tuning for optimization. The parameter ξ is assumed to be 0.7 for simplicity. In this study, additional conventional triaxial compression tests for each aspect ratio were conducted, with varying conditions such as density, confining pressure, and drainage conditions. From these results a critical state was identified  Fig. 13 Critical state line at e − p plane and the relationship between AR and the critical state can be explored further. Figure 13 shows the critical state lines obtained from the above tests in the e − p plane. It can be seen that simulation data can be fitted by the lines with the expressions given in Eq. (10) very well for each case. This demonstrates that the critical state for a given AR is independent of density, confining pressure and drainage conditions. The critical state lines converge, indicating that the influence of the particle aspect ratio is negligible when particle elongation continues to decrease. Yan and Dong [39] also discussed that the location of the critical state line in e − p space is sensitive to particle grading. By incorporating the crushable particles in DEM simulation, Muir Wood and Maeda [40] found that the changing of particle grading can cause a downward parallel shift in the critical state line in e − p space without change in slope. Recently, the effect of particle breakage has been considered in critical state framework due to changes in the particle grading, which may occur in practical geotechnical engineering, e.g. pile driving and large earth-fill dams. Through experimental tests on sands, both Bandini and Coop [41] and Ghafghazi et al. [42] confirmed that the particle breakage can cause a shift of critical state line in e − p space, while disagreed with the mode of position changes, i.e. translational offset or rotational shift. In the present study, the critical state line obtained in DEM results is lower in position than in the physical experiments. According to previous studies shown above, it can be interpreted due to the difference in grading of particles, the simplified particle shape modelled and/or due to particle crushing in the experimental tests which is ignored in these DEM simulations. Figure 14 presents the critical state line in the p − q plane. The data points can be fitted by straight lines passing through the origin, with the slope of the straight line being the critical stress ratio M c . For the sample with AR = 1, which is composed of spherical particles, M c has the smallest value = 0.770, which is consistent with the findings by Yan [43] and Gu et al. [44]. The critical state ratios for various AR are summarized in Fig. 15, and can be approximately fitted by a straight line, suggesting that samples with smaller AR (e.g. bigger elongation) have greater resistance to shear and result in large critical angles of shearing resistance. To verify this observation, the experiment data from Yang and Luo [6] are also added on this diagram, showing a similar qualitative trend between critical state ratio M c and the particle aspect ratio AR, which ranges from 0.70 to 0.85 in their experiments.
Recently, an anisotropic critical state theory (ACST) was presented by Li and Dafalias [19] to enhance the classical CST by introducing a condition that fabric must satisfy at the critical state. According to Li and Dafalias [19], the norm F of F evolves towards a unique critical state value while the direction n F evolves towards the loading direction. In addition, it has been shown that the fabric tensor F should be a per-volume measure and the critical state norm of fabric tensor F would be a function of Lode angle θ F only, according to ACST [35]. However, F defined by Eq. (5) using the measurements from DEM analysis might not be a per-volume measure and the easiest approach is to normalize F with the specific volume v of the sample, which can be written as This lead to such normalized critical state fabric norm F c is a unique value and only depends on the tensor's Lode angle θ F . It is noted that as specified in [17], the plastic part of the specific volume v p = (1+e p ) should be used to replace the total specific volume v in Eq. (11) as a rigorous normalizer. However, the elastic part of the void ratio is relatively insignificant and can be neglected for simplicity. Such treatment of using total specific volume will hardly affect the results presented in this paper. Figure 16a presents the critical values of F c for the triaxial compression simulations considering varying AR values. It is seen that F c depends on the particle shape, and the samples that composed of more elongated particles appear to have greater critical state norm, while the sample with spherical particles has lowest F c . As expected, the un-normalized critical values of F c paired with critical state e c increase with the void ratio or density of the samples, as indicated by the dashed lines, confirming that the DEM fabric tensor defined in Eq. (5) is not a per-volume measure. However, the adaption of DEM fabric tensor by normalization given in Eq. (11) gives promising results, as illustrated in Fig. 16b [16,44] are presented in Fig. 17a for the un-normalized critical state norm against e c , indicating a polynomial trend with the F c increasing with e c , which is similar and consistent with the observations shown in Fig. 16a. By applying the same normalization rule expressed in Eq. (13) for all the data from different DEM simulations indicate more scattered and biased data for AR = 1 than those for other two cases. This is thought to be the coupling effects between specific volume v and the tensor's Lode angle on the critical state norm F c [35], and the simple normalization rule expressed by Eq. (11) does not adequately tackle the adaption of DEM fabric tensor for samples involving spherical particles. Nevertheless, the particle shape dependent critical state fabric norm F c after normalization can reasonably satisfy the conditions advocated by anisotropic critical state theory.

Concluding remarks
This paper presented a systematic DEM analysis of the effect of particle aspect ratio on the mechanical behavior of granular materials. A series of constant cell pressure and constant volume conventional triaxial compression tests, as well as constant p and constant b true triaxial test numerical experiments were performed. DEM simulations considered specimens with varying aspect ratio, and a unique particle size distribution. It was shown that the aspect ratio of the particles has a significant influence on the shear strength, dilatancy and failure of granular materials. The microscopic responses were also examined in terms of the fabric evolutions during shear. Combined with the conventional critical state theory (CST) and the recent anisotropic critical state theory (ACST), the critical state behavior was discussed. The main observations can be summarized as follows: 1. In DEM, clumped particles can be used to assess the effects of aspect ratio on granular material response. It was shown that while considering different initial densities and confining pressures, a critical state and fabric evolution can be attained.

Both conventional triaxial compression tests and constant
p and constant b true triaxial tests provided evidence that particle aspect ratio can significantly affect the stressstrain behavior of granular materials. 3. The η and η * f parameters that define the shape of failure criterion in principal stress and fabric space are a function of AR. The mobilized strength is a function of both AR and b. 4. It has been established that the particle aspect ratio can influence the critical state properties. The position of critical state line in e − p plane is insensitive to particle aspect ratio. However, the critical state line in p − q plane appears to be sensitive to particle aspect ratio. 5. In ACST, the critical fabric norm should be unique, that is only dependent on Lode angle θ F of the stress path, for a given soil. The quantified fabric was adapted by a simple normalization that can reasonably satisfy this additional requirement by ACST. A unique critical state fabric norm was found whose magnitude is dependent on the particle shape but is independent of critical state void ratio. 6. Broad agreement of critical state fabric norm with those obtained from other DEM studies suggests that the simple normalization rule adopted in this study deserves further discussion when considering for samples composed of spherical particles.
Note that the particle shape in the present study is created by clumped particles with varying overlapping for simplicity. The effects of the convexity of particle shape and result-ing resistance of particle rotation are not considered and may have influence on the mechanical responses of the soils. This could be improved by employing more realistic particle shapes in the DEM simulations. Such work may involve complex algorithms and significant computational cost. Although this issue is out of the scope of this paper, it deserves further investigation.