How particle shape affects the critical state, triggering of instability and dilatancy of granular materials – results from a DEM study

Nguyen et al. (2021) studied the influence of the shape of a particle on critical state (CS) and dilatancy using discrete-element method (DEM) modelling. Triaxial compression of assemblies of spheres, ellipsoids and two clusters of spheres were investigated. One of the results is the linear stress ratio–dilatancy relationship for minimum dilatancy states. The stress ratio for zero dilatancy represents CS ratios (M_c). In the case of a small number of DEM calculations, the CS stress ratio (CS angle, ϕ_cs) may be affected by error.

The linear relationship between the stress ratio (η) and the minimum dilatancy (d_min) represents a frictional state line (FSL) for conventional drained triaxial compression and has the form (Szypcio, 2016a)

where

Here, ϕ^o is the critical frictional state angle (Szypcio, 2016a, 2016b), δε_v = δε₁ + 2δε₃ and δε_q = 2(δε₁ − δε₃)/3. The critical frictional state angle for granular materials is equal to the CS angle ϕ^o = ϕ_cs, (M^o_c = M_c; Szypcio (2016b)). Points representing minimum dilatancy states lie on the FSL. The slope of the FSL (A^o_c) in the η–d_min plane is defined by the CS friction angle ϕ^o (ϕ_cs) or M^o_c (M_c). The best approximations of stress ratios–minimum dilatancy obtained with the use of equation (13) are with M^o_c = M_c = 0·67, 0·87 and 1·17 (ϕ^o = ϕ_cs = 17·54°, 22·33° and 29·31°) for spheres, ellipsoids and clusters of spheres, and these are shown in Fig. 15.

The corresponding values of M_c (ϕ_cs) obtained by Nguyen et al. (2021) are: M_c =0·8, 0·9, 1·2 and 1·1 (ϕ_cs =20·67°, 23·04°, 30·0° and 27·7°). There are some differences between the values of M_c (ϕ_cs) calculated with the use of equation (13) and direct calculations, especially for small numbers of DEM simulations.

The special case of Rowe's equation shown in the Appendix (equation (12))

is correct only for δε_v/δε₁ = −1 (Horne, 1965; Rowe, 1969). For triaxial compression

Comparing equations (13) and (18) for δε_v/δε₁ = −1 the relationship between sin ϕ_μ and sin ϕ_cs has the simple form

Figure 16 shows the relationship between ϕ_μ and ϕ_cs calculated by equation (21). The correct values of ϕ_μ are: ϕ_μ =12·9°, 16·9°and 22·15° or, in terms of the corresponding friction coefficients, μ = 0·23, 0·30 and 0·41 for spheres, ellipsoids and clusters of spheres. The inter-particle friction angles shown in Fig. 12 (Nguyen et al., 2021) do not fulfil equation (21).

The DEM method has gained great popularity in the last decade. The stress–dilatancy relationship, equation (13), can be used to determine the CS stress ratios (CS friction angle) correctly, especially for small numbers of DEM simulations. The friction coefficient as a basic DEM modelling parameter may be calculated by use of equation (21). The discussion contributors believe that the relationships presented will be helpful for DEM modelling in the future.

The authors appreciate the discussers’ interest in their paper (Nguyen et al., 2021). The discussers show a way of approximating the CS stress ratio (M_c) from a linear approximation of the stress ratio (η) and the minimum dilatancy (d_min) relationship. This was done on the assumption by the discussers that a limited number of DEM simulations were conducted, and thus M_c was affected by the dataset, despite it having been stated in the paper that there were a large number of simulations for each particle shape, in the section ‘Effect of particle shape’; the discussers might have missed this information. However, the authors welcome this discussion to clarify the issue. The authors also note that the discussers’ approach would have some uses when independent measurement of the CS stress ratio from individual tests/simulations is not available or possible to achieve, but a dilatancy relation is, somehow, available. Therefore, these have been explained below.

The CS stress ratio (η_c = q/p′) is an independent measurement of q and p′ at CS from an experiment/simulation. The data trend from a few to many individual experiments/simulations is plotted in q–p′ space, and the slope of the trend line is taken as the M_c for that soil (see Fig. 17). Therefore, the number of data points could affect the value of M_c. In the authors’ DEM study, both drained and undrained triaxial simulations after isotropic and K₀ consolidations were considered. All these simulations gave stress ratios at CS. Thus, the corresponding M_c values were obtained from 34 simulations for spheres, 104 for ellipsoids and 46 for clusters (see Fig. 17). These numbers are certainly higher than the drained simulation data points used by the discusser in Fig. 17 and the authors’ M_c values represent a wide range of conditions. Most triaxial experiments could not afford as many data points as the authors have considered through the DEM simulations. It is well known that the measurement of M_c or a close approximation to the M_c value in q–p′ space for both experiments and simulations is much easier than for the other CS parameters associated with void ratio (e) or volume change (ε_v) in e–log(p′) space. In this respect, it is not clear to the authors why M_c has to be estimated from d_min measurements, which are associated with volume changes. The authors can think of a situation for the use of this approach – when all measurements, including ε_v, can be done reliably in triaxial experiment/simulation, but CS cannot be reached for various reasons, for example, strain limitations.

However, the authors would be cautious in using this approach even in the above hypothetical situation. This is simply because M_c is the intercept on the η-axis when d_min = 0 for $η=Mco−Acodmin$⁠. The discussers’ Fig. 17 shows no data points close to d_min = 0; to be precise most data points were d_min < −0·40, except for the ellipsoid. Thus, the best-fit linear relations were weighted with the data points away from M_c for d_min = 0 in the case of the sphere and cluster particles. However, more data points, including d_min ≈ 0, were used for the best fit for ellipsoids. Therefore, the estimated M_c of 0·87 was very close to the authors’ M_c of 0·90.

In the study under discussion, the authors did not explore any theoretical relationship between ϕ_μ and ϕ_cs. They are not clear how that justified the correctness of the discussers’ M_c or ϕ_cs. The authors explored a relationship between $ϕμ,dmin$ at d_min and d_min. The value of $ϕμ,dmin$ was calculated based on the direct measurement of major principal stresses and strains at the point of d_min and using equation (12) in the paper. The same method was done by Frossard (1979) to calculate ϕ_μ, who suggested that ϕ_μ increased with particle angularity. Fig. 12 in the original paper demonstrated the same finding. Further, a generalisation of the ϕ_μ and ϕ_cs relationship in the discussers’ equation (21) and Fig. 18 was also misleading – for example, having a CS friction angle of ϕ_cs = 0 is not realistic for soil. Furthermore, such a relationship was developed based on a liner relationship between η and d, where d_min coincides with d at η_peak. The authors found that d_min occurred well after η_peak, as shown in Fig. 18 and in other literature (Rahman et al., 2018; Nguyen et al., 2020). Such a phenomenon may be related to particle angularity, but that needs further investigation. However, this observation is not only DEM related. The experimental work by Santos et al. (2010) also showed the difference between d_min and d at η_peak.

This work was conducted at Bialystok University of Technology (Poland) and was supported by the Polish Financial Resources on Science (project no. WZ/WB-IIL/2/2022).

A^o_c

slope of the frictional state line

d

dilatancy

d_min

minimum dilatancy

M_c

critical state stress ratio

p′

mean effective stress

q

shear stress

ε_q

shear strain

ε_v

volumetric strain

ε₁, ε₃

major, minor principal strain

η

stress ratio

μ

friction coefficient

σ₁, σ₃

major, minor principal stress

ϕ_cs

critical state angle

ϕ_μ

interparticle friction angle

ϕ^o

critical frictional state angle