Exploring the micro-to-macro response of granular soils with real particle shapes by way of μCT-aided DEM analyses Open Access

For granular materials the term ‘particle shape’ informs the overall form, local angularity and surface texture of a particle in three different aspects of observation (Barrett, 1980). Particle shape can be quantitatively described by shape parameters. Krumbein & Sloss (1963) provided an intuitive shape chart based on sphericity (S) and roundness (R). Altuhafi et al. (2013) applied dynamic image analysis (DIA) to enable the calculation of sphericity (S), convexity (C_x) and aspect ratio (AR). Researchers have continued to develop new shape parameters to better capture the role of particle shape in geomechanics, including overall regularity (OR) (Yang & Luo, 2015) and shape angular group indicator (SAGI) (Altuhafi et al., 2016). Angelidakis et al. (2022) proposed a new interpretation of elongation, flatness and compactness to better describe the particle form in a three-dimensional (3D) perspective. In this contribution, sphericity (S), convexity (C_x), AR and OR are employed and their definitions are given in the Appendix, where greater values of these particle shape parameters indicate a more sphere-like/less irregular shape, whereas lower values indicate the opposite.

Within a small-strain range (typically ≤10⁻⁵), small-strain shear stiffness (G₀) can be expressed as a function of mean effective stress (p′) and void ratio (e) as

where a and b are material-dependent constants that determine G₀, and F(e) represents the void ratio correction function proposed by Hardin & Richart (1963). For an extensive deformation where the intermediate principal stress is equal to the minor principal stress under dry and drained triaxial compression, granular materials exhibit a unique critical state characterised by two equations:

where q is the deviator stress q = σ′₁ − σ′₃, in which σ′₁ and σ′₃ are major and minor principal stresses, respectively; M, e^Γ, λ are critical state parameters.

The significance of particle shape in the above-mentioned soil characteristics from the small-strain regime to the critical state is widely recognised by experimentalists. Cho et al. (2006) and Altuhafi et al. (2016) revealed that more irregular particles can attain higher maximum (e_max) and minimum void ratios (e_min). Cho et al. (2006) reported that increasing particle irregularity reduces the small-strain stiffness (lower a) while increasing stress sensitivity (greater b). However, Altuhafi et al. (2016) found that G₀ normalised by the coefficient of uniformity (C_u) tends to increase with irregularity. The latter observation was supported by experiments on granular mixtures composed of different fractions of particle shapes (Shin & Santamarina, 2013; Liu & Yang, 2018). In addition, Chang et al. (1991) proposed a value of b = 1/3 for Hertzian contacts between spheres, while experimental results gave a range of b from 0·32 to 0·72 (Cho et al., 2006; Altuhafi et al., 2016). A consensus on the effect of particle shape on G₀ has not been achieved because of these contradictory findings.

In relation to critical state parameters, Cho et al. (2006) reported that increasing particle irregularity leads to decreases in e^Γ and λ; however, Altuhafi et al. (2016) found the same trend for e^Γ whereas a less obvious relation between λ and particle shape. Yang & Luo (2015) mixed Fujian sands and glass beads and found decreasing trends of M and λ as particles become more irregular. Xiao et al. (2019) mixed two different-shaped glass beads and found increases in critical stresses for more irregular particles. Wu et al. (2021) tested clinker ash in drained triaxial tests and observed decreasing trends of e^Γ and λ as particles become less irregular. Overall, the sensitivity of M, e^Γ, λ to particle shape has become clearer with the advancement of experimental technology.

In contrast, discrete-element method (DEM) studies that use discs or spheres lack natural particle rolling resistance and particle interlocking which crucially affect the mechanical characteristics of real granular materials (Zhao et al., 2023). To better represent realistic soil behaviour, non-spherical particles such as super-ellipsoids (e.g., Zhao et al., 2017), sphere-clusters (e.g., Nguyen et al., 2021), sphere-clumps (Katagiri et al., 2010), level-set functions (Kawamoto et al., 2016, 2018) and spherical harmonics (e.g., Xu et al., 2021) have been used in DEM simulations. Zhao et al. (2017) found that the densest packing density appears at AR = 0·3, but not in the case of spheres (AR = 1). Nguyen et al. (2021) examined the sensitivity of critical state parameters to particle shape and found decreasing trends of M and λ with S, in line with Yang & Luo (2015); while there was an increasing trend of e^Γ with S, different from that in Cho et al. (2006). Xu et al. (2021) observed an increase in shear strength in both peak and critical states with increasing particle irregularity. They also found that more irregular materials tend to develop more significant stress-induced fabric anisotropy during the loading process.

X-ray micro-computed tomography (μCT) enables the acquisition of detailed images of granular materials, allowing for the replication of complex shapes in simulation. Zhang et al. (2021) generated non-spherical particles based on Ottawa sands and found that particles made of more irregular materials can achieve looser packing conditions (i.e., higher e) during isotropic compression. Nie et al. (2021) simulated Leighton Buzzard sand (LBS) and reported that the mean coordination number (⁠$Z¯$⁠) reduces while the mean contact force grows with the decrease of OR. Otsubo et al. (2023) used a μCT image of an LBS particle to investigate systematically the influence of grain shape and surface asperity on G₀. They revealed that particle morphology imposes a difference in total contact area between particles, leading to a variation of G₀. However, these DEM studies did not consider a wide variability of particle shapes typically observed in geotechnical experimentation. To this end, practical protocol of simulation campaigns needs to be established to investigate the universality of particle shape effects on the mechanical behaviour of granular materials.

In the current paper, μCT is applied to acquire 3D images from a wide selection of real granular materials. Non-spherical clumps are generated to provide multi-sphere representations for these real particles. DEM simulations of drained triaxial compressions are conducted from small to large strains using the generated clumps, to investigate how particle shape affects the macroscopic mechanical response and microscopic soil fabric evolution throughout the triaxial loading. Finally, the aim of this contribution is to establish a micro-to-macro connection between critical deviator stress ratio and internal fabric anisotropy, linking experimental and simulation observations which incorporate particle shape.

Spherical, deformed, clumped and angular glass beads (SGB, DGB, CGB and AGB), silica sand (SS) and Leighton Buzzard sand (LBS) were scanned; the abbreviated names were used to label and distinguish different materials in the laboratory experimentation, as outlined in Li et al. (2024). The resolution of the images was 17·62 μm for glass beads, 11·65 μm for SS and LBS in voxel sizes, respectively. This means that a sand grain with a diameter of 1 mm is represented by approximately 86 voxels across its diameter. Owing to the high shape variability of CGB and AGB, two grains were selected for each, noted as CGB1, CGB2 and AGB1, AGB2, respectively. Therefore, eight different shapes were selected to investigate the effects of particle shape on micro-to-macro responses in this contribution.

Figure 1 shows the main image processing operations for a slice of the silica sand particle. The 3D μCT image was first binarised by identifying the solid phase of each material using the thresholding technique in Otsu (1979). To separate neighbouring touching particles, the watershed segmentation algorithm was adopted in three dimensions as detailed in Nadimi & Fonseca (2018). After segmentation, each particle was marked with a unique label with individual particle geometry information. The labelled μCT images of the materials of interest were analysed to identify all the different particles and approximate their morphology using multi-sphere representations – that is, sets of rigidly connected overlapping sub-spheres. The use of this method feasibly preserved the concavity of non-spherical particles to capture the particle interlocking. This operation was performed using the open-source software Code Library to generate universal multi-sphere particles (CLUMP) (Angelidakis et al., 2021a). The multi-sphere clumps were generated based on the Euclidean distance transform of 3D images. Each sub-sphere was generated at a location where the mass of the particle was most represented in each iteration of the procedure. According to the strategy, the first sub-sphere was always the largest possible inscribed sphere, while each subsequent one was smaller than its previous iteration. With this method, particles can be generated at predefined fidelity levels with a given overlap ratio between sub-spheres, which was adjusted as 0·2 to 0·45 in this paper. This overlap ratio determined the ratio between the overlap distance and the radius of the later generated sub-sphere for two connected sub-spheres in every iteration based on the Euclidean distance transform.

A surface extraction algorithm within the CLUMP software allows for the tessellation of the surface of the generated multi-sphere particles, which can be used to quantify their shape parameters. Utilising the labelled particle images of the analysed materials, shape characteristics of the generated clumps and original particles from the μCT were calculated using the open-source software SHAPE (shape analyser for particle engineering) (Angelidakis et al., 2021b). Here, the minimal-volume oriented bounding box of each particle was used to calculate its main particle dimensions. It becomes apparent in the literature that the same shape parameter may be defined in many different ways. To this end, the formulae used to calculate S, C_x, AR and OR in this contribution are listed in the Appendix.

To achieve an optimal balance between simulation efficiency and morphological fidelity of the non-spherical particles in DEM simulations, the number of sub-spheres per particle (N_s) needs to be decided so as to minimise the difference in shape parameters between the clump and real particles.

Figure 2 demonstrates the variation of S, C_x, AR and OR for generated clumps with N_s = 1, 2, 3, 5, 10, 25, 100, 200, 500 using an SS particle example, where the images of clumps are provided to visualise how particle morphology evolves with N_s. The generated clump fills more volume enclosed by the surface mesh of the original SS particle with the increase of N_s. All shape parameters drop markedly with the increase of N_s; in particular, S, C_x and OR decrease until N_s = 25, followed by a plateau for larger N_s; AR noticeably reduces from N_s = 1 to N_s = 2 and slightly increases with increasing N_s. After similar comparisons, the other considered shapes present stable values of shape parameters no later than N_s = 25, where only a small relative variation can be observed. Table 1 lists the shape parameters of the reconstructed clumps with N_s = 25, along with those for μCT images of the corresponding real particles. All parameters of clumps, except AR, have similar values to those of μCT images within a relative difference smaller than 7%, revealing a reliable representation of original particle morphology. It appears that more satellite sub-spheres may be needed for CGB2 and AGB2 to accurately reconstruct clump dimensions closer to the real particle in the three principal directions, which leads to AR values closer to those of real particles.

In addition, Otsubo et al. (2023) applied the same algorithm to generate clump LBS particles and found that the stress responses become similar for simulations throughout the entire triaxial loading with N_s ≥ 10. Katagiri (2019) conducted a series of isotropic compression tests using clumps to simulate Toyoura sand and found that the void ratio and mean coordination number in the isotropic state do not change significantly beyond a threshold of N_s = 14. According to the parametric study and relevant literature, N_s = 25 was decided as a sufficient number of sub-spheres on the safe side to generate representative non-spherical clumps for subsequent simulations, specifically for the particle shapes considered in this contribution.

Figure 3 presents non-spherical clumps of the studied materials with N_s = 25 sub-spheres, except for SGB, which was represented by a single sphere. The clump generation technique of Favier et al. (1999) was used to create the clumps for DGB because of its axisymmetric shape. Referring to Table 1, SGB takes a value of 1·0 for all shape parameters, while CGB2 shows the lowest C_x value because of its special cluster formation; AGB1 and AGB2 have the lowest AR due to their planar shapes.

This contribution adopted a modified version of LAMMPS (Plimpton, 1995; Nguyen & Plimpton, 2019), which allowed for a robust implementation of the simplified Hertz–Mindlin (HM) contact model (Itasca Consulting Group, 2007). In comparison to the linear contact model, the HM contact model can capture more realistic non-linear contact behaviours (Johnson, 1985). The simplified HM contact law was applied to each sub-sphere in non-spherical clumps; this scheme enables the local asperities of the non-spherical clumps to be reflected in the contact behaviours. Specifically, the smaller sub-spheres in contact exhibit a softer contact stiffness with a given contact force. Therefore, the use of clumps facilitated a more systematic investigation on how the mechanical response evolves from a sphere to a non-spherical clump. The material properties were a material density, ρ_p = 2500 kg/m³, Young's modulus, E_p = 64·5 GPa and Poisson's ratio, v_p = 0·23. Following O'Sullivan & Bray (2004), a density scaling factor of 1000 was applied to enhance computational efficiency. The material properties were uniformly assigned to each sub-sphere, where the mass and moments of inertia of overlapped regions were excluded. A linear particle size distribution (PSD) was identically assigned to every material as 0·6–1·2 mm. The minimum volumetric ellipsoid method was applied to generate packings made of non-spherical clumps (Zheng & Hryciw, 2017), where circumscribed spheres were generated using the given PSD, which were then replaced by the non-spherical clumps used in each simulation scenario.

Periodic boundaries were employed in the three principal directions of the simulated samples. Using periodic boundary conditions, when a particle crosses one boundary of the simulation domain, it reappears on the opposite boundary with the same velocity and material properties. These boundary conditions offer the advantage of avoiding initial heterogeneity in soil fabric such that the number of particles for a sample can be efficiently optimised (Thornton, 2000). The sample was composed of 5000 spheres for SGB and 5000 non-spherical clumps (i.e., 1·25 × 10⁵ nominal sub-spheres). All clumps were generated in a cubic space without initial particle contacts in the absence of gravity. Subsequently, the cloud of particles was compressed under an isotropic stress of 50 kPa using a servo-control scheme, which ensured a stable equilibrium for the generated non-spherical clumps following Morimoto et al. (2021). All granular samples were compressed with the interparticle friction coefficient, μ_iso = 0·001, which ensured the densest condition in the DEM environment, serving as an analogy of relative density = 100% for a given shape. A representative cubic sample is given in Fig. 4, which consists of the non-spherical clump made of N_s = 25 using the imaging data of an SS particle.

After undergoing isotropic compression, all samples were loaded in the vertical direction (1-axis in Fig. 4) under dry and drained conditions, considering an interparticle friction coefficient, μ_shear = 0·35. The triaxial loading was performed at a constant strain rate of 0·002 (s⁻¹), while maintaining the minor principal stress as σ′₂ = σ′₃ = 50 kPa. This loading rate was employed to ensure quasi-static loading conditions as in Lopera Perez et al. (2016). Local damping (ζ_local) and viscous damping (ζ_vis) were omitted in the isotropic compression; however, ζ_local = 0·1 and ζ_vis = 0·05 were applied during triaxial compression to suppress kinetic energy and minimise dynamic effects. In addition, ζ_local and ζ_vis were deactivated during small-strain probing tests to evaluate G₀. According to Otsubo et al. (2017), a timestep (Δt) of 2 × 10⁻⁶ s was selected for simulations with SGB, and the same was used for non-spherical clumps.

Figure 5 compares the deviator stress (q) and volumetric strain (ε_v) for all tested materials with axial strain (ε_a) under the densest condition; the images of non-spherical clumps are annotated for reference. Overall, granular samples demonstrate significantly different responses due to different particle shapes. SGB shows the lowest stress level and minimal dilation. In contrast, DGB initially experiences a noticeable increase in q, reaching its peak value the earliest among all particles, while displaying the second lowest stress level in the critical state. AGB1 and AGB2 show less significance of strain softening, even under the densest conditions. ε_v for AGB2 gradually develops to show the latest onset of the phase transformation point, but demonstrates the second most significant dilation at ε_a = 40%. CGB2 gives the highest stress levels, particularly in the peak state and the most marked dilation upon shearing, while LBS and SS show intermediate stress levels and volumetric changes among these materials. It is noted that the initially dramatic increase in q (shown in the inset) that is different from conventional laboratory experiments is attributed to the manual adjustment of μ_shear, as explained by Li et al. (2021).

Figure 6 explores the changes in void ratio (e) and mean coordination number (⁠$Z¯$⁠) with ε_a. e continuously increases, while $Z¯$ dramatically drops initially and stabilises with straining. This initially substantial reduction in $Z¯$ results from the rapid collapse in horizontal contacts in the granular assembly, while the subsequent stable $Z¯$ value mainly reflects the particle contacts in the loading direction to form force chains, as elucidated by Li et al. (2022a). CGB2 and AGB2 exhibit much higher e values throughout the loading; however, they show the greatest stress levels in the peak and critical states. This highlights the more dominant role of the particle shape compared with initial density in the stress response while both are interplaying. DGB and CGB1 attain a lower e value than SGB in the isotropic state, but end up with an e value similar to SGB in the critical state. SGB and DGB display a markedly lower $Z¯$ than other particles at later stable stages after the initial reduction. For irregular particles, CGB2, AGB1 and AGB2 show a higher $Z¯$ compared with LBS and SS.

Figure 7(a) reveals the relations between peak/critical deviator stress ratios ((q/p′)_peak, (q/p′)_cri) and shape parameters (S, C_x, AR and OR). (q/p′)_peak and (q/p′)_cri increase as the employed shape parameters reduce, indicating that more irregular particles have higher shear resistance in the peak and critical states. Indicated by the semi-transparent shading, (q/p′)_cri can be better correlated than (q/p′)_peak with shape parameters; specifically, AGB gives a lower (q/p′)_peak than CGB2 despite having similar OR. This indicates that peak stress is more related to the particle interlocking or sliding/rolling motions. It is noted that ε_v,cri with shape parameters is not provided, as its trend is similar to that of void ratio with shape parameters under dry and drained conditions.

Figure 7(b) compares the relations between isotropic/critical void ratios (e_iso and e_cri) and shape parameters. Samples containing more sphere-like particles can achieve denser packings in the isotropic and critical states. e_iso here can be considered as an analogy of densest attainable void ratio (e_min) in the laboratory. This decreasing trend is consistent with experimental results in Altuhafi et al. (2016). As shown in Fig. 7(b) and Table 2, e_iso values for DGB and CGB1 are lower than that of SGB, e_cri for DGB is also lower than that of SGB. Compared with SGB, DGB and CGB1 have similar S and C_x, but a lower AR. This is in agreement with Salot et al. (2009), who reported that slightly elongated particles can achieve denser packings than idealised spheres.

Figure 7(c) explores the changes in mean coordination numbers in the isotropic and critical states (⁠$Z¯iso$ and $Z¯cri$⁠) with shape parameters. A two-phase variation of $Z¯iso$ with shape parameters can be observed: CGB2, LBS and AGB1, AGB2 demonstrate $Z¯iso$ being around 12 despite variations in their shape parameters, whereas others show a decreasing trend in $Z¯iso$ for more sphere-like materials. This is because asperity-to-asperity contacts were quantified rather than clump-to-clump contacts – in other words, identifying contacts between sphere elements of a clump were identified instead of those between non-spherical clumps; Otsubo et al. (2023) reported that $Z¯iso$ of LBS particles is around 12 for asperity-to-asperity contacts, but around 8 for clump-to-clump contacts at their densest packing configuration.

Small-strain shear stiffness (G₀) was estimated under isotropic stresses of 50, 100, 200 and 400 kPa using static probing – that is, the linear stress–strain relation within a strain range of 10⁻⁵ (Kuwano & Jardine, 2002). According to Hardin & Richart (1963), the void ratio correction function in equation (1) was interpreted as F(e) = (B − e)²/(1 + e). Following Iwasaki & Tatsuoka (1977), Kuwano & Jardine (2002) and Liu & Yang (2018), this contribution used B = 2·17 for all tested materials.

According to Fig. 8, a in equation (1) decreases with the increase of shape parameters, which is contradictory to the experimental observations summarised in Cho et al. (2006). Liu & Yang (2018) used materials with identical uniformity coefficients (C_u) and argued that the effect of C_u on G₀ should not be neglected, as also pointed out by Wichtmann & Triantafyllidis (2009) and Li et al. (2022b). In the present paper, identical C_u was used for all materials so that the reduction in a is consistent with that in Liu & Yang (2018). However, b values (equation (1)) are found to be insensitive to shape parameters, with b ≈ 0·35 (Table 2) for all materials. This value is close to the theoretical value of 1/3 based on the Hertzian contact law and effective medium theory (Chang et al., 1991), irrespective of particle shape. Otsubo et al. (2023) found that b is sensitive to N_s, implying that surface topography is governing the stress sensitivity of G₀, which was effectively supported by Otsubo & O'Sullivan (2018). In the current paper, a consistent N_s = 25 was used for non-spherical clumps, suggesting that varying N_s is needed to capture b of experimental data.

The critical state parameters can be estimated according to equations (2) and (3) using the samples mentioned above, which were sheared to the critical state (ε_a = 40% in this contribution). The critical state parameters are plotted against shape parameters in Figs 9 and 10, along with reference data from the literature (Xie et al., 2017; Jiang et al., 2018; Nguyen et al., 2021; Nie et al., 2021). In Figs 9 and 10, the equations of the linear trendlines are provided with the determination of coefficient of determination (R²) to indicate the evolving trends with reliability. M decreases with the increase of shape parameters, in agreement with Nie et al. (2021) and Nguyen et al. (2021). Yang & Luo (2015) experimentally reported a qualitatively decreasing trend for M using Fujian sands. λ and e^Γ also show decreasing trends as particles become less irregular. This trend aligns with the results observed in Nie et al. (2021) while Nguyen et al. (2021) displayed a contradictory trend in their DEM simulations. In contrast, e^Γ gives a good comparison between the present paper and past studies. In terms of R², it appears that C_x gives the lowest value, which indicates that the convexity conditions of the particles may not dominantly govern the mechanical responses in the critical state, but contribute to the particle interlocking, as more evidently reflected in the peak stress strength.

According to Fig. 6, $Z¯$ approaches an identical value for each particle shape in the critical state. Similarly, equation (4) is proposed to capture a unique relation between $Z¯$ and p′:

where $Z¯Γ$ and $λZ¯$ are the intercept and slope in the $Z¯$–p′ relation with varying confining stresses in the critical state. Figure 11 portrays the variations of $Z¯Γ$ and $λZ¯$ with shape parameters. A lower $Z¯Γ$ is achieved for more sphere-like materials for a given stress level in the critical state. $λZ¯$ decreases with the increase of particle shapes, indicating that more irregular materials have a higher sensitivity of $Z¯$ to stress level in the critical state.

The microstructure network of granular samples is described by various directional fabric indices, including contact normal (CN), particle orientation (PO) and branch vector (BV). CN refers to the vector normal to the contact plane for asperity-to-asperity contacts; PO describes the non-spherical particle's major axis orientation; BV is defined as the vector connecting the centroids of two contacting particles.

Following Fonseca et al. (2013), ψ is defined as the inclination between the unit vectors of the above-mentioned soil fabric indices and the horizontal plane (2–3 plane), as depicted in the insets in Fig. 12; further, ψ^CN, ψ^PO and ψ^BV are the angles relating to CN, PO and BV, respectively. Fig. 12 shows the normalised frequency distribution of ψ^CN, ψ^PO and ψ^BV at selected ε_a throughout loading for SS. In this context, the vertical axis refers to the normalised frequency, indicating the ratio between the number of vectors within an examined angle range and the total number of vectors. The theoretical normalised frequency distribution is derived from an enumeration of the unit vectors that describe all the orientations originating from the centre of a unit sphere. By comparison, the distributions in the isotropic state closely align with the theoretical distribution, confirming the homogeneity of samples in terms of soil fabric. When load is applied, particle contacts are progressively dominant in a ψ^CN value between 30° and 40°. Similarly, ψ^BV distributions gradually show the peak between 15° and 25° with loading. Differently, ψ^PO shows a rather monotonic distribution with fewer and fewer particles orientated with a larger ψ^PO. Throughout loading, an increasing number of particles align their major axes with ψ^PO < 30°.

A second-order tensor is utilised to quantify soil fabric (Oda, 1982):

where N is the total number of vectors existing in the granular system; and n^k_i and n^k_j are the unit orientation vectors alongside the i and j directions, respectively. Referring to Fig. 4, Φ₁₁ represents the component of the soil fabric tensor in the loading direction while (Φ₂₂ + Φ₃₃)/2 stands for that perpendicular to the loading direction where Φ₂₂ ≈ Φ₃₃ was confirmed.

Figure 13 visualises the 3D distributions of the three normalised soil fabric indices for SS at sequential loading stages, following the interpretation by Zhao & Zhou (2017). At ε_a = 0%, CN, PO and BV display approximately symmetric polar distributions, signifying the fabric isotropy prior to the initiation of triaxial compression. Subsequently, CN becomes more pronounced in the vertical direction, and its orientation remains relatively stable in the later stages of loading. In contrast, PO exhibits more noticeable distributions in the horizontal plane. Interestingly, BV initially evolves to be more vertically orientated, but this anisotropic tendency diminishes as the distribution turns to be less significant after ε_a = 13·0%.

To quantify the evolution of soil fabric indices, the fabric ratio is calculated as Φ = Φ₁₁/[(Φ₂₂ + Φ₃₃)/2] to represent the fabric anisotropy of the granular samples during shearing. Φ^CN, Φ^PO and Φ^BV are defined as fabric anisotropy indices for CN, PO and BV, respectively. Fig. 14 illustrates the variations of Φ^CN, Φ^PO and Φ^BV against ε_a upon loading using the examples of SS and CGB2. Φ^CN increases and subsequently tends to be stable for both SS and CGB2. In contrast, Φ^PO continuously decreases upon loading and tends to be stable in the critical state, giving a value consistently lower than 1. As non-spherical clumps were randomly created with an isotropic Φ^PO (i.e., Φ^PO = 1), this discloses that non-spherical particles undergo particle rearrangement during shearing and tend to align their major axes perpendicular to the loading direction. Interestingly, the Φ^BV of SS and CGB2 shows different variations: Φ^BV of SS gradually increases and then gently decreases, supporting the visualisation in Fig. 13, whereas that of CGB2 initially increases but reduces markedly with further shearing and appears to evolve below the fabric isotropy line, which indicates that BV of CGB2 is first vertically significant but subsequently horizontally significant with shearing. This can be attributed to the fact that BV is not influenced only by contact condition but also by particle shape. More irregular materials with a low AR may exhibit horizontally significant fabric anisotropy concerning BV as loading continuously induces particle rearrangement.

The deviator part of the fabric tensor can be defined as

A scalar parameter can then be introduced according to the invariant of Φ′_ij as follows:

Furthermore, a^CN_c, a^PO_c and a^BV_c are defined as anisotropy degrees for CN, PO and BV, respectively.

Figure 15 provides the variations of a_c for CGB2 and SS in terms of CN, PO and BV during loading. In comparison to Fig. 14, the variations of a^CN_c are similar, gradually increasing and eventually stabilising. The variations of a^PO_c increase progressively prior to a relatively mild trend. Especially, a^BV_c for SS first increases and then decreases as similar to the variation of Φ^BV; however, a^BV_c for CGB2 shows special transitional variations where a secondary increase in a^BV_c is observed after about ε_a > 28%. Referring to equation (7), a_c is a scalar parameter where a_c = 0 means fabric isotropy. In this sense, a_c does not convey spatial information regarding whether the soil fabric vectors are more oriented in the vertical direction or the horizontal direction. For instance, particle orientation is more significant in the horizontal direction as the non-spherical particles align their longer axis in the horizontal plane during triaxial compression. The later increasing part of the transitional point for a^BV_c of CGB2 reveals that the branch vector starts to become more significant in the horizontal plane.

The values of Φ^CN, Φ^PO and Φ^BV at ε_a = 40% are termed as critical fabric anisotropy (Φ^CN_cri, Φ^PO_cri and Φ^BV_cri), and their variations are correlated with OR in Fig. 16. As SGB has a completely symmetric shape, so Φ^CN and Φ^BV are collinear, whereas Φ^PO can always be treated as isotropic. The evolution of Φ^CN_cri shows a scattering feature, being not clearly related to OR. In particular, the significantly high value of Φ^CN for AGB2 is due to the identification of asperity-to-asperity contacts for the CN estimation, where multiple asperity-to-asperity contacts were counted for two contacting non-spherical clumps due to its planar shape, when their bodies are contacting other than at their tips. Nie et al. (2021) reported that the fabric anisotropy in terms of CN is more significant for more irregular materials, in which a more systematic alteration of particle morphology was considered. By contrast, more irregular materials exhibit more pronounced Φ^PO_cri during loading, being more significant in the directions perpendicular to the loading direction. Φ^BV_cri also displays an increasing trend with OR with a similar slope to that of Φ^PO_cri plotted against OR, while crossing the fabric isotropy line. This means that Φ^BV_cri can be more significant either in the vertical or horizontal direction. More irregular materials tend to distribute more BV in the horizontal directions, while more sphere-like materials show an opposite effect.

Figure 17 compares the critical stress ratios (q/p′)_cri of simulation results obtained from this contribution to experimental results originally from Li et al. (2024) extrapolated to (q/p′)_cri following Yilmaz et al. (2023). The simulation and experimental results are qualitatively comparable to reveal that more irregular particles lead to a higher stress level in the critical state. The quantitative differences in this comparison can be attributed to the following reasons. (a) The samples had different shapes for each individual particle, leading to distributions of shape parameters like PSD (Li et al., 2024); whereas the DEM simulation considered the mono-dispersed shape distribution in each DEM sample. To highlight this, Fig. 17 provides horizontal bars to represent ±1 standard deviation of the particle shape distribution as outlined in Li et al. (2024). (b) Boundary conditions were different, with periodic boundaries in this study, but lateral membrane conditions with rigid top and bottom walls in the laboratory (Li et al., 2021).

A shape-weighted fabric anisotropy index (ORΦ^PO_cri or ORΦ^BV_cri) is newly proposed by combining OR with fabric anisotropy in the critical state. Fig. 18 shows that (q/p′)_cri is well correlated by the proposed index, indicating an essential role of shape-incorporated soil fabric in the stress response. Fig. 16 demonstrate that the microstructure of soils in the critical state is significantly affected by particle shape. Meanwhile, Li et al. (2021) revealed that elastic wave velocities in the critical state are influenced by the critical soil fabric at least for spherical particles. Wen & Zhang (2022) discovered the existence of a fabric surface in the critical state regardless of the shearing mode, where (q/p′)_cri is linearly correlated with strong-contact soil fabric. It can be concluded that the macroscopic response of soils is internally governed by microscopic parameters that are influenced by particle shape in the critical state. Therefore, incorporating the evolution of particle shape is essential when establishing constitutive models (e.g., Nguyen et al., 2020). In practice, since it is trivial to determine (q/p′)_cri physically and numerically, Fig. 18 builds the bridge towards predicting the microscopic critical soil fabric based on the macroscopic stress ratio in the critical state.

This contribution has systematically investigated the micro-to-macro response of granular materials with different particle shapes under the densest possible packings and monotonically sheared to the critical state in a DEM virtual environment. Based on the simulation results, the following conclusions can be drawn.

• μCT is a powerful tool to reconstruct non-spherical particles for use in DEM simulations. With the given imaging data obtained from the granular materials tested in Li et al. (2024), N_s = 25 is found to be a sufficient number to preserve 93% fidelity of S, C_x and OR, whereas a greater number of N_s is required to obtain AR closer to that of the originally irregular particle.

• Particle shape has a profound effect on the mechanical response of granular materials from small to large strains. As the particle becomes more irregular, q_peak and q_cri, e_iso and e_cri, $Z¯iso$ and $Z¯cri$ exhibit increasing trends.

• Small-strain shear stiffness constant a decreases with increasing values of the monitored shape parameters; however, b seems to be less affected by particle shape but is influenced more by surface topography. M, λ and e^Γ are affected by particle shape and give lower values with the increase of shape parameters; $Z¯Γ$ and $λZ¯$ are also found to be shape sensitive, with a decreasing tendency as shape parameters increase.

• The evolution of soil fabric demonstrates different behaviours for different particle shapes during shearing. CN progressively becomes more significant in the loading direction while PO shows the opposite trend during loading. By contrast, BV initially increases while stabilising or decreasing with further shearing, influenced by particle shape. In the critical state, samples with more sphere-like particles present less significant fabric anisotropy in terms of PO; BV demonstrates a transition from a more vertically to a more horizontally significant fabric response as particles become more irregular.

• (q/p′)_cri from DEM simulations captures the laboratory results that used the granular materials for the non-spherical clump generation. Furthermore, the shape-weighted fabric indices (ORΦ^PO_cri or ORΦ^BV_cri) are newly introduced to capture linearly the (q/p′)_cri values, which is expected to predict microscopic soil fabric in the critical state using macroscopic stresses.

The aim of this contribution was to explore the role of particle shape on the micro-to-macro characteristics of granular soils, with a particular focus on a more realistic reflection of their responses. For this, μCT was utilised to acquire high-resolution imaging data, which ensured accurate 3D shape characterisation and construction of non-spherical particles while preserving their morphological fidelity. A vastly wide selection of particle shapes used in laboratory experimentation was considered for a more general representation of granular materials in real-world conditions. In addition, the simplified HM contact model, which more accurately captured the non-linear contact characteristics, was applied considering the local asperities of the multi-sphere clumps. Several future directions are envisioned. First, there is a need to develop a constitutive framework that implicitly incorporates particle shape terms (e.g., Vu et al., 2024). This contribution empirically proposed a correlation between (q/p′)_cri and the shape-weighted fabric indices; this can be further refined by introducing a shape variable that governs the critical state parameters into the constitutive law. Second, analytical studies (e.g., Rothenburg & Bathurst, 1989; Seyedi Hosseininia, 2013) have proposed the intrinsic relation between macroscopic deviatoric stress and microscopic shear-induced fabric anisotropy using idealised particle shapes. It would be interesting to numerically validate their theoretical relations using the more complex particle morphology adopted in this contribution. Another promising direction involves exploring the combined effect of particle shape and surface roughness. Naturally occurring soils have multiscale particle morphology (Barrett, 1980), while this contribution has focused on the coarsest scale particle shape. Future work could use μCT and spherical harmonics to broadly assess particle morphology from particle shape to surface roughness (e.g., Zhou et al., 2015; Nie et al., 2024).

This contribution was funded by the Japan Society for the Promotion of Science (22K14322 and JPJSBP120195701), the UK Royal Society (IEC\R3\183026) and partly by the UK Engineering and Physical Sciences Research Council (EPSRC) grant number EP/V053655/1. DEM simulations were performed using the Fujitsu Primergy CX400M1/CX2550M5 (Oakbridge-CX) in the Information Technology Centre at the University of Tokyo, Japan.

a, b

material-dependent constants to determine G₀

a_c

anisotropy degree

C_x

convexity

e

void ratio

e_iso, e_cri

e in the isotropic state and critical state

G₀

small-strain shear stiffness

M, e^Γ, λ

critical state parameters

N_s

the number of sub-spheres to form one single clump

p′

mean effective stress = (σ′₁ + 2σ′₃)/3

q

deviator stress = σ′₁ − σ′₃

q/p′

deviator stress ratio

q_peak, q_cri

q in the peak state and critical state

S

sphericity

$Z¯$

mean coordination number

$Z¯iso$⁠, $Z¯cri$

$Z¯$ in the isotropic state and critical state

$Z¯Γ$⁠, $λZ¯$

proposed critical state parameters for $Z¯$

ε_a

axial strain

ε_v

volumetric strain

ζ_local, ζ_vis

local and viscous damping coefficients

μ_iso, μ_shear

interparticle friction coefficient during isotropic compression ( = 0·001) and triaxial compression ( = 0·35)

σ′₁, σ′₃

major and minor principal stresses (σ′₃ = 50 kPa)

Φ

fabric ratio = Φ₁₁/[(Φ₂₂ + Φ₃₃)/2]

Φ₁₁, Φ₂₂, Φ₃₃

second-order fabric tensor in the three principal directions

ψ

inclination of the soil fabric vectors relative to the horizontal plane

This appendix provides the formulae for calculating the shape parameters (S, C_x, AR and OR) especially used in this contribution based on Fig. 19.

where V_A and V_C are the volumes of a particle and its convex hull, respectively. A is the surface area of the particle. S_l, I_l and L_l are the short, intermediate, respectively. For a 3D object, flatness (S_l/I_l) and elongation (I_l/L_l) can describe a particle form in more spatial aspects; here the AR is considered as S_l/L_l = (S_l/I_l)(I_l/L_l).