Coupled effects of particle overall regularity and sliding friction on the shear behavior of uniformly graded dense sands

Jiyn Nie, Shiwei Zho, Yifei Cui,*, Yu Wng

a State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing,100084, China

b South China Institute of Geotechnical Engineering, South China University of Technology, Guangzhou, 510640, China

c Key Laboratory of Mountain Hazards and Surface Process,Institute of Mountain Hazards and Environment,Chinese Academy of Sciences,Chengdu,610041,China

d University of Chinese Academy of Sciences, Beijing,100049, China

Keywords:Discrete element method (DEM)Overall regularity (OR)Particle sliding friction Peak and critical friction angles Fabric anisotropy Dense sands

A B S T R A C T

1. Introduction

Particle morphology has been regarded as an important factor affecting mechanical behaviors of granular soils,such as the stressstrain relationship and shear strength, which has been proved through a great number of experimental studies (Cho et al., 2006;Yang and Luo, 2015; Alshibli and Cil, 2018; Xiao et al., 2019; Shi et al., 2020a). With the rapid development and increasing popularization of X-ray micro-computed tomography (μCT) in geotechnical engineering,different particle shape aspects(e.g.overall form,local roundness)proposed by Barrett(1980)can be quantified very well for real three-dimensional (3D) particles except for the roughness due to the limited CT scanning resolution (Zhao and Wang, 2016; Zhou et al., 2017). However, most of the current experimental researches adopted scanning electron microscopic(SEM) tests or dynamic image processing system to acquire the equivalent two-dimensional (2D) shape parameters for their high efficiency and relatively low cost,particularly in statistical analyses which require a larger number of particles(Su and Yan,2020).The corresponding research results may deviate from real situations since the 2D shape parameters determined from a projected 2D image along a pre-defined or a random direction cannot reflect the authentic morphological characteristics of the original 3D particles.To this end, a few of attempts have been made to discuss the correlations between 2D shape parameters and their 3D counterparts for regular geometries or natural aggregates particles (Cavarretta et al., 2009; Kutay et al., 2011; Fonseca et al., 2012; Alshibli et al.,2015; Cepuritis et al., 2017; Yan and Su, 2018; Su and Yan, 2020).Nevertheless, it is necessary to use the 3D shape parameters directly to describe the effect of particle shape on soils’responses to avoid the unnecessary dispersion of these correlations in the future work.

Alternatively, discrete element method (DEM) has become one of the most useful numerical tools to explore the effect of particle morphology on the shear behaviors of granular soils in 3D cases from different shape aspects. For example, Azéma et al. (2013)studied the influence of particle’s non-convexity on the rheological properties of granular materials. Zhao et al. (2015) and Zhao and Zhou (2017) used an open-source non-spherical DEM code,SudoDEM, to investigate the effects of particle angularity and asperity on the direct shear and triaxial shear behaviors of granular soils,respectively.Gong and Liu(2017)adopted the popular PFC3D(Itasca, 2015) to discuss the effect of particle aspect ratio on the triaxial shear behaviors of multi-sphere ellipsoid assemblies. Xie et al. (2017) also explored the effect of the particle aspect ratio on the shear behaviors of sandy soils using the open-source DEM software, Yade. In addition, Kodicherla et al. (2019) discussed the influence of particle elongation on the direct shear behaviors of granular materials using DEM. Recently, Xu et al. (2021) numerically investigated the effect of irregularity and elongation on the shear behaviors of granular assemblies combined with Fourier shape descriptors. These studies highlighted the influence of different particle shape aspects on the behaviors of granular soils.However, a more robust shape descriptor needs to be considered when conducting shape-related DEM studies so that the obtained research findings can better serve for the natural sands.Recently,a well-defined shape descriptor named overall regularity (OR) has been adopted by many researchers to quantify the effect of the particle overall irregularity on the mechanical behaviors of sands both experimentally and numerically (Yang and Luo, 2015, 2018;Xiao et al.,2019;Nie et al.,2021a,b).It is defined as the average of the particle aspect ratio,sphericity and convexity,and can quantify both the particle overall form and local roundness characteristics synthetically.It should be noted that this shape parameter ignores the roughness effect, and all the related research is case-studied and cannot reach the consensus for the distinct particle roughness characteristics of different types of sands.

On the other hand,the particle roughness is generally observed through SEM and quantified by the root mean square of the asperity heights(Ren et al.,2021).Empirical formulae between the particle sliding friction coefficient and particle surface roughness and other material properties have also been built to quantify the effect of roughness on the friction behavior of contacting bodies(Sandeep and Senetakis,2018).In order to incorporate the particle roughness’s effect into DEM modeling,sliding friction coefficient is generally assigned onto DEM particles.Based on this,a series of 3D numerical studies on the influence of the sliding friction coefficient on the shear behaviors of spherical or ellipsoidal particle assemblies have been conducted (Yang et al., 2012; Huang et al., 2014;Gong et al.,2019a).Research results showed that the particle sliding friction had a great influence on the peak and critical state friction angles of granular soils, and the peak friction angle increases with the sliding friction,while the critical friction angle firstly increased and then saturated at a constant value. Moreover, the sliding friction thresholds above the critical friction angles were different comparing the results of both previous spheres and ellipsoids DEM studies. Different physical mechanisms have been explored to reveal these interesting phenomena from perspectives of the buckling load of a simple particle column model and geometrical and mechanical fabric anisotropies, respectively (Huang et al.,2014; Gong et al., 2019a). At this point, it is an open issue to understand the underlying mechanism for more realistic particles with concave shapes. Therefore, a systematic numerical investigation regarding the coupled effects of particle OR and sliding friction is needed to lead to a better understanding of the role of the particle shape on the mechanical behaviors of the granular soils from different shape aspects.

Motivated by the above mentioned research purposes, a series of 3D drained triaxial compression tests on representative volume elements(RVEs)with different OR and sliding friction is conducted using 3D DEM. The macro- and micro-mechanical behaviors of these assemblies during triaxial shear tests are explored to get a cross scale understanding of the soil responses governed by different particle shape aspects.The remaining part of this paper is organized as follows. Firstly, the adopted numerical model is introduced, followed by the analysis of macro-mechanical results.Then, the micro-mechanical behaviors are studied to explain the observed OR- and sliding friction-dependent shear strengths.Finally, the major conclusions in this study are obtained.

2. Numerical model

Drained triaxial shear tests were numerically conducted in this study to explore the influence of both the particle OR and sliding friction on shear behaviors of the uniformly graded dense sands.To consider the morphology characteristics of natural sand particles more reasonably in DEM,3D clump technique was adopted in this study to simulate the concave-convex profiles of 5000 Leighton Buzzard sand(LBS)particles.The distinct LBS particle morphologies were firstly generated randomly based on a limited number of Xray μCT scanned LBS particles (i.e.100 particles or so) through our previously proposed random generation method considering both grading and shape characteristics(Nie et al.,2020,2021a).Then,the 3D clumps were automatically generated within PFC3D based on the “Bubble Packing” algorithm (Taghavi, 2011), taking the imported LBS particle morphologies as templates. This algorithm has been adopted by many researchers to consider realistic shape particles in DEM(e.g.Wei et al.,2018;Shi et al.,2020b,Zhang et al.,2021a). In DEM, the translational and rotational kinematics of the irregular clump was respectively determined by the generalized Newton’s and Euler’s equations asFi=mdvi/dtandMi=Iidωi/dt-(Ij-Ik)ωjωk,whereFiandMiare the resultant force and torque acting on or around the clump centroid, respectively;mis the clump mass;Ii,IjandIkare the principal moments of inertia along theith,jth, andkth principal directions, respectively; and viand ωiare the translational and angular velocities along or around theith direction, respectively. For the ideal sphere, the Euler’s equation would be simplified asMi=Iidωi/dtdue toIj=Ik. To calculate the mass centroid position, mass and inertia tensors of these irregular clumps, the following strategy is adopted in PFC3D(Itasca, 2015). When the clump is generated and the spatial distribution of its internal pebbles is determined,the space containing the clump is recursively subdivided into a certain number of cubes.The volumes and inertia tensors of all cubes that fall completely within any pebble are accumulated. The parallel axis theorem is used to accumulate the inertia tensors. The iteration with the increasing number of cubes would cease if the relative error of the estimated clump volumes between any two successive iterations is within the desired tolerance (e.g. 0.005 in this study). The mass centroids of all cubes are finally averaged to acquire the final mass centroid position of the clump assuming that the material density is uniform within the clump. In addition, only the contact forces between particles are considered to contribute to the resultant force and torque applied on the clump,since the gravity is ignored in this study.

Fig.1a presents a typical clump morphology used in this study although it may somewhat deviate from that of natural sand particles especially at the scale of particle corner and edge.In order to quantitatively describe the particle morphology characteristics,different shape indices, e.g. equivoluminal sphere diameter (de),aspect ratio (A), sphericity (S), and convexity (C) are employed to quantify the multi-scale characteristics of the particle morphology according to the previous works (Yang and Luo, 2015, 2018; Zhao and Wang, 2016). Among these,derepresents the particle size while the latter three indices reflect the overall form and local roundness characteristics according to the three aspects of particle morphology in Barrett(1980).It is noted that more realistic particle morphology with advanced contact detection algorithm is needed for a refined DEM study of the specific natural sand(e.g.Kawamoto et al., 2016; Wang et al., 2021). Fig.1 gives the definitions of these morphology indices and their cumulative probability distribution(CDF) curves. It is seen that these 3D clumps can be used to simulate the shape and size poly-dispersity characteristics of coarse sub-rounded sand particles in nature.The grading characteristics of these numerical samples are uniform with a uniformity coefficient of 1.27 determined based onde. Noted that Monte Carlo sampling method is adopted to compute the surface area(Ac)and volume(Vc)of these clumps,while three principal dimensions,P1,P2andP3,and convex hull volume,Vch,are calculated by surface meshing method according to our previous work(Nie et al., 2021a).

To prepare granular assemblies with different particle OR characteristics, a certain number proportions (0%, 25%, 50%, 75%,and 100%) of 3D clumps are randomly replaced by equivoluminal spherical particles.The particle OR(ρo)for each granular assembly is determined as ρo= ρocωc+ ρosωsaccording to (Yang and Luo,2015, 2018), where ρocand ρosare the representative ORs of the clumps and spheres, respectively; ωcand ωsare the mass proportions of clumps and spheres, respectively. ρocis defined as the average of aspect ratio, sphericity and convexity corresponding to 50%cumulative probabilities for 3D clumps,as shown in Fig.1b.ρosis defined for spheres similarly with ρoc,and equals 1 since all shape parameters of different shape aspects are 1 for spheres following their definitions. In this study, the granular assemblies with 5 different OR values are prepared according to this strategy,and the values of OR are 0.896, 0.922, 0.948, 0.974 and 1, respectively,representing the granular assemblies with low to high mass proportions of spheres. A larger OR value means that the particle morphology characteristic of the granular assembly is more regular.Table 1 summarizes the representative OR values and mass proportions of the clumps and spheres for the 5 granular assemblies.

For the triaxial shear simulations, firstly, 5000 particles consisting of the clumps and ideal spheres of different mass proportions were randomly generated into a cubic shear box with no initial overlap. The cubic shear box is made of 6 frictionless rigid walls. Then, an isotropic compression procedure is initiated to simulate the ‘consolidation’ process in the laboratory. During this process, the particle sliding friction is assigned to 0 to ensure that the final mechanically stable assemblies can achieve a densest state with a given confining stress σ3= 100 kPa (Deluzarche and Cambou, 2006; Abbireddy and Clayton, 2010). Hence, it can be assumed that the initial relative densities of these 5 types of granular assemblies are consistent with each other,which ensures that their peak shear strengths are mainly affected by the particle OR and sliding friction.The corresponding relationships built in this study can provide a reasonable reference for the future particle shape related DEM studies of very dense sands. The granular assembly is assumed to reach the mechanical equilibrium state once the ratio of the mean unbalanced force to the mean contact force is less than 1 ×10-5(Masson and Martinez, 2001; Zhou et al., 2018;Gong et al., 2019b).

Fig.2a shows the‘consolidated’assembly filled with about 50%mass proportions of the clumps and spheres, where the final size(L0,W0andH0)of the consolidated sample is about 33 mm,11 times of the maximum particle diameterdmax, meeting the basic requirement for RVE in DEM to avoid the rigid boundary effect(Jamiolkowski et al., 2004). After the consolidation, a constant loading rate() is applied on both the top and bottom walls along thezdirection to compress the whole sample while the movements of other 4 boundary walls inxandydirections are constantly adjusted under a numerical stress-servo mechanism, in order to maintain a constant confining stress of σ3= 100 kPa during the whole compression process. The loading rate () is chosen as 0.01 s-1in this study to ensure that the inertial numberis less than 0.001 and the whole shear process is quasi-static(MiDi,2004;Gu et al.,2020),where ρ is the density of granular particle andis the average particle diameter. It is noted that the particle sliding friction coefficient during shearing is reassigned from 0.001 to 1 in order to investigate its influence on the shear behaviors of irregular concave particles.The compression is terminated at the axial strain of 40%, where an obvious critical state is observed, and the shear stress and volume of the whole sample are not changed with the continuous shearing.

During shearing, the stress tensor σijcan be written as(Christoffersen et al.,1981):

wherefciandbcjaretheithcomponentofthecontactforcevectorfandjth component of the branch vectorb,respectively,at the contactc.

The branch vectorbjoins the centroids of two contacting particles. Based on the stress tensor σijin Eq. (1), the effective mean stresspand deviatoric stressqcan be determined asp=，andq=, where σ′ijis the deviatoric part of σij, i.e. σ′ij=σij-pδij, and δijis the Kronecker delta tensor.

The axial strain εaand volume strain εvare derived from the changes of sample heighthand sample volumevduring shearing,which are calculated as follows (Zhao and Zhao, 2019):

whereHandH0are the heights of sheared sample and initial sample, respectively;VandV0are the volumes of sheared sample and initial sample, respectively; and a positive εvmeans a compressive deformation following the convention of traditional soil mechanics.

The simple linear contact model is used in this study, which is much more computationally efficient than nonlinear contact models(e.g.Hertz-Mindlin)and can yield comparable results(Zhao et al., 2018). Table 2 summarizes the corresponding model parameters. The density scaling strategy is adopted to acquire a relatively large critical time step Δtcrto improve the computational efficiency according to the study of O’Sullivan (2011). The computational time step Δtin this study is finally set to be 1×10-5s for all simulations.The normal(kn)and tangential(kt)contact stiffnesses are set askn/r* =kt/r* = 100 MPa, wherer*is the average equivoluminal particle radius of the 5 assemblies(Guo and Zhao,2013;Zhao and Zhao, 2019). In order to accelerate the dissipation of the kinematic energy of the whole granular system, a local damping system with the damping ratio of α = 0.3 is adopted.

Table 1 Representative OR values and mass proportions of clumps and spheres for 5 types of granular assemblies.

Table 2 DEM model parameters for the drained triaxial shear simulations.

Fig.1. Particle shape simulation and quantification: (a) Definitions of particle size and shape indices; (b) CDFs of particle shape indices of realistic shape clumps; and (c) CDF of particle size of realistic shape clumps (modified from Nie et al., 2021b). Cu is the coefficient of uniformity;de60, de50 and de10 are the de values corresponding to 60%, 50% and 10%cumulative probability, respectively; and A50, S50 and C50 are the A, S and C values corresponding to 50% cumulative probabilities, respectively.

3. Macro-mechanical properties

3.1. Stress-strain characteristics

Fig. 3 shows the stress-strain curves for dense samples with different sliding friction coefficients when ρo=0.948 under drained triaxial shear conditions. It can be seen that the particle sliding friction has a noticeable effect on the deformation characteristics,shear strength and volume deformation of dense granular soils.Interestingly,for dense samples with ρo=0.948 and μs≤0.01,the deviatoric stressqincreases with the axial strain until stabilizing at a constant value with the slight fluctuation, which is the so-called‘strain-hardening’. In addition, the volume deformation of these dense samples during shearing is very limited(e.g.within 1%)and the dilation will diminish when μsapproaches 0. In contrast, the typical‘strain-softening’behaviors are observed for dense samples with ρo=0.948 when μs≥0.05,that is,the deviatoric stress q firstly climbs to a maximum before decreasing and stabilizing. Meanwhile, all samples experience significant volume dilation. Moreover,the shear strength and volume dilation are enhanced with the increase of the particle sliding friction.Similar phenomena are also found in other cases with a given ρoand different μsvalues. Fig. 4 presents the drained shear behaviors for the dense samples with different ρowhen μs=0.5.Typical drained shear responses of dense sands in the laboratory are found for all samples,and both the shear strength and volume deformation decrease with the particle OR.As revealed from Figs. 3 and 4, it can also be seen that the deviatoric stressqand volumevfor all numerical samples are kept relatively steady when the axial strain is beyond 30%, suggesting that all samples almost reach the critical state.

3.2. Peak and critical friction angles

In order to quantitatively discuss the effects of the particle OR and sliding friction on the shear strength of granular soils,the peak and critical friction angles, i.e. φpand φcfor all the samples with different combinations of ρoand μsare determined as φp= sin-1and φc= sin-1,where σ1pand σ1care the peak and residual axial stresses along thezdirection, respectively; σ3is the confining stress along thexorydirection, which equals 100 kPa herein.It should be noted that in order to alleviate the fluctuation influence of the axial stress at the critical state,the average of the axial stresses at the axial strain ranging from 30%to 40% is chosen as the final σ1cto determine the critical friction angles. Figs. 5 and 6 present the relationships between φp, φcand μsand ρo,respectively.It can be seen that φpincreases with μs,while decreases with ρomonotonously.However,φcfirstly increases with μsand then stabilizes when μsexceeds a threshold,which seems to be OR-dependent and increases with ρo.Similar to φp,φcdecreases with the increase of рo due to the weak interlocks between irregular particles. Similar trend has also been observed in the experimental work of Yang and Luo (2015), which is also presented in Fig.6b.It can be speculated that the particle sliding friction for the Fujian sands and glass beads used in their experiments due to the surface roughness may fall in between 0.2 and 0.3, which is also consistent with the micro-mechanical experiment values (e.g.0.15-0.23) on common quartz sands (Sandeep and Senetakis,2018). A slight deviation may be attributed to the definition differences of the particle OR in both studies, and the 2D shape parameters are adopted in Yang and Luo(2015)while in this study,3D shape parameters are used. Hence, Fig. 6 provides a reasonable reference for the critical friction angles of uniformly graded sands with different particle ORs and sliding frictions. Hence, different mechanisms are involved in the influences of the particle OR and sliding friction on the critical friction angles,which will be explored from the perspective of micro-mechanical behaviors in the following sections.

Fig. 2. Numerical modeling configuration (a) before and (b) after shearing(L0 ≈W0 ≈H0 ≈11dmax).

4. Micro-mechanical behaviors

4.1. Contact network characteristics

To explore the microscopic origin of the particle OR-and sliding friction-dependent peak and critical friction angles, the contact network characteristics of the granular samples with different sliding frictions (ρo= 0.948), and different particle ORs (μs= 0.5)are analyzed. Figs. 7 and 8 give the evolutions of the mean coordination numberand mean normal contact forcen.is defined as the average of the contact numbers over all the particles, whilenis the average of the normal contact forces over all the contacts.Due to the very dense characteristic of these granular samples,decreases upon shearing due to the instant dilation and finally declines to a critical value for all the samples.For the samples with μs≤0.01,changes slightly before and after shearing.In addition,evolutions ofnduring shearing are similar to those of the deviatoric stress, as shown in Figs. 3a and 4a. Moreover, similar trends are observed forwith μsand ρo,i.e.with the increase of μsand ρo,decreases, and hencefnincreases, suggesting that ignoring the concave shape characteristics of the particles in DEM may underestimateor overestimaten.Actually,the phenomenon of largerfor the granular assembly consisting of more irregular particles has also been observed by the experimental studies aided by the X-ray CT scanning and image processing techniques (Fei and Narsilio,2020; Zhang et al., 2021b). It further suggests that the coordination number can be chosen as a microscopic index for quantifying the accuracy of DEM studies considering the realistic particle shape’s characteristic other than the common macroscopic soil properties.Fig.9 shows the evolutions of their corresponding mean mobilized friction coefficientm,which is defined as the average of the mobilized friction coefficients μmover all the contacts, and μm=ft/fn, whereftis the tangential contact force.malso has a similar trend to the deviatoric stress, andmincreases with μswhile decreases with ρo. This phenomenon explains the OR- and sliding friction-dependent peak and critical friction angles drawn from Figs. 5 and 6, except that the critical friction angles will be saturated when μsis larger than a certain OR-dependent threshold.Since from Fig.9b,we can see thatmincreases monotonously with μs for samples of ρo= 0.948 at the critical state. Therefore, the micro-mechanisms of the OR- and sliding friction-dependent critical friction angles need to be further explored in the following section.

4.2. Anisotropic fabric characteristics

Granular soils resist the external load through the inter-particle contact network, which means that during loading, the soil particles will adjust themselves to form a well-organized internal fabric structure to support the loading better. Fig. 10a shows the typical contact force networks for the granular sample of ρo= 1 and μs=0.5 at three characteristic states:initial,peak and final states.At the initial isotropic stress state after‘consolidation’,the contact force uniformly distributes in the sample.It is noted that due to the rigid boundary effect,the stiffness of‘wall-particle’contacts is larger than that of ‘particle-particle’ contacts, and thus some sporadic strong contact force chains develop around the confining rigid walls.At the peak state under shearing, it can be seen that the contact force network behaves a distinct bimodal characteristic, that is, strong contact force chains pass through the top and bottom walls to support the axial loading. Meanwhile, there are more weak-contact force chains intersecting or even being perpendicular to the strong contacts in the sample to stabilize the strong contact force chains.With further shearing and sample dilating,the contacts dissolve and hence the contact force network seems to be sparser,but an obvious anisotropic and bimodal contact networks can be still seen with a less mean normal contact force.The right subfigure of Fig.10a gives the corresponding contact force network at the critical state.

In order to visualize the spatial distribution of normal contact forces, Fig. 10b shows the 3D histograms of local average normal contact forcenshown in Fig. 10a at three characteristic states referring to Zhao and Zhou (2017). Each bar represents the normalized local average normal contact forcen/0in that direction,wherenis defined as the average of the magnitudes of the normal contact forces that fall in the direction interval of the corresponding bar, and0is the average ofnover all bars. From Fig. 10b, it can be easily seen that the spatial distribution of the normal contact forcefnis almost isotropic at the initial state due to the isotropic compression, and then becomes elongated along the loading direction during the following shearing, corresponding to the strong contact force chains forming in the loading direction as observed in Fig.10a. Moreover, a larger shear strength at the peak state corresponds to a more elongated histogram than that at critical state. This means that the shear strength of the granular soils comes from the fabric anisotropy within the sample. Fig.10c and d presents the corresponding 3D histograms of the sample with ρo= 1 and μs= 0.05 and ρo= 0.896 and μs= 0.5 at three characteristic states, respectively. Comparing those with Fig. 10b,we can find that both the particle OR and sliding friction affect the spatial distribution of the normal contact forces,i.e.the anisotropy of the normal contact force fabric. With this in mind, the micromechanisms of the particle OR- and sliding friction-dependent critical friction angles are explored in the following from the perspective of the anisotropic fabric.

Fig. 3. Stress-strain relationships for samples with different μs when ρo = 0.948 under drained shear: (a) q-εa and (b) εv-εa.

Fig. 4. Stress-strain relationships for the samples with different рo when μs = 0:5 under drained shear: (a) q-εa and (b) εv-εa.

Fig. 5. Relationships between φp and (a) μs and (b) ρo.

Fig. 6. Relationships between φc and (a) μs and (b) ρo.

Fig. 7. Evolutions of (a)Z and (b)fn of the samples with different μs when ρo = 0.948.

Fig. 8. Evolutions of (a)Z and (b)fn for the samples with different ρo when μs = 0.5.

Fig. 9. Evolutions ofm for samples with (a) different ρo when μs = 0.5 and (b) different μs when ρo = 0.948.

To further characterize the anisotropic fabric and quantify its anisotropy degree,there are commonly 5 types of fabric tensors to describe the internal geometrical or mechanical fabric structures of the granular assembly, i.e. fabric tensor of contact normal,normal and tangential branch vectors,, and normal and tangential contact force vectors,(Oda,1972,1982;Iwashita and Oda,1999), which can be determined based on DEM data (Zhao and Zhou, 2017; Wang et al., 2020; Nie et al., 2021a):

whereniis theith component of unit contact normal vectorn,whiletiis theith component of unit contact tangential vectortat the contactc;bnandbtare the projections of branch vectorbon thenandtat the contactc,respectively;fnandftare the projections of branch vectorfon thenandtat the contactc, respectively;Ncis the total number of the contacts within the numerical assembly;Acklis a specific deviatoric fabric tensor characterizing the anisotropy of fabric tensor of contact normal vector,i.e.the anisotropy tensor of the contact normal vector,which is determined as follows(Zhao and Zhou, 2017; Nie et al., 2021a):

where φc′

klis the deviatoric part of， in which δklis the Kronecker delta tensor.

Likewise, the anisotropy tensors of the normal and tangential branch vectors,, and the anisotropy tensors of the normal and tangential contact force vectors,are defined as (Zhao and Zhou, 2017; Nie et al., 2021a):

In general, a scalar value,a*is used to quantify the anisotropy degree of the anisotropy tensor, i.e. the deviatoric invariant of the tensor (Zhao and Zhou, 2017; Nie et al., 2021a):

where the subscript or superscript“*”represents the contactc,bn,bt, fn, and ft in Eqs. (9)-(11),σ′ijis the deviatoric part of the stress tensor σij, and has been defined in Eq. (1). The positive value of sign(A*ij) indicates that the principal direction of the anisotropy tensorA*ijis consistent with that of σ′ij, and vice versa. The definitions of different types of the fabric tensors can be referred to the relevant literature (Oda, 1972, 1982; Iwashita and Oda, 1999) for more details.

Fig.10. Illustration of the contact force networks and 3D histograms of local average normal contact force at the initial,peak and final states for the typical samples:(a,b)Sample with ρo = 1 and μs = 0.5; (c) Sample with ρo = 1 and μs = 0.05; and (d) Sample with ρo = 0.896 and μs = 0.5.

The stress-force-fabric relationship proposed by Ouadfel and Rothenburg (2001) provides a feasible solution to approximate the critical friction angle from the fabric anisotropies, as shown in the following equations (Zhao and Zhou, 2017; Nie et al., 2021a):

whereMcis the critical deviatoric stress ratio;qcandpcare the deviatoric stress and effective mean stress at the critical state,respectively;acis the anisotropy degree of the contact normal fabric;abnandabtare the anisotropy degrees of the normal and tangential branch vector fabrics, respectively;afnandaftare the anisotropy degrees of the normal and tangential contact force fabrics, respectively. A largera*means that the corresponding fabric is more anisotropic and contributes to a higher critical friction angle.In addition,Eq.(13)has been validated based on the samples with 5 OR values and μs=0.3 under the same simulation conditions in our previous work (Nie et al., 2021a). On the other hand, Eq. (14) reveals that φcincreases withMcmonotonously.Therefore, herein these 5 types of fabric anisotropies are directly compared between samples with different рo and μs so as to give some insights into the micro-mechanisms of the particle OR- and sliding friction-dependent critical friction angles.

Figs.11-14 give the relationships betweenac,afn,aftandabnand рo,μs at the critical state.It can be seen from Figs.11a,12a and 13a and 14a that with the increase of μs,acfirstly increases with μs and then reaches a plateau when μs is larger than a certain particle ORdependent threshold for the samples of ρo= 1 and ρo= 0.974.However,when рo further decreases,the evolution ofacwith μs at the critical state will experience a‘softening’trend,that is,acfirstly climbs up to a peak point and then declines with μs, which is different from the previous observation for irregular elongated clump samples in Gong et al.(2019a).And the result of Gong et al.(2019a) is the same as that of the samples with ρo= 1 and ρo= 0.974. One possible reason for this difference may be the concave morphology characteristics considered in this study,which are more close to the realistic granular soils.For the evolution ofafnwith μs, a similar softening trend can be obviously seen for all samples,whileaftincreases monotonously with μs.In addition,abnkeeps at a constant level with a slight fluctuation regardless of μs.Sinceabtis rather limited compared with other four fabric anisotropies, it is not presented herein. Based on the evolutions of these fabric anisotropy degrees with μs, it is concluded that the saturation of the critical friction angle φcwith μs comes from the decrease ofacorafnbeing compensated by the increase ofaftfor irregular concave particles.

Fig.11. Relationships between ac at critical state and (a) μs and (b) ρo.

Fig.12. Relationships between afn at critical state and (a) μs and (b) ρo.

Fig.13. Relationships between aft at critical state and (a) μs and (b) ρo.

Fig.14. Relationships between abn at critical state and (a) μs and (b) ρo.

On the other hand,from Figs.11b,12b and 13b and 14b,it can be seen that these 4 fabric anisotropies increase or decrease with the particle OR monotonously,and the total decrease ofac,aft,andabnis much larger than the increase ofafn, and hence the critical friction angle φcdecreases with ρo,as shown in Fig.6b.The microscopic origin of the OR- and sliding friction-dependent critical friction angles has been successfully revealed from the fabric anisotropies so far.

5. Conclusions

This study investigates the coupled effects of the particle OR and sliding friction coefficient on the shear behaviors of the uniformly graded dense sands systematically based on the 3D DEM simulations.Five types of RVEs with different particle ORs are prepared by mixing irregular concave clumps with ideal spheres of different mass proportions but the same grading. The stress-strain relationships and internal friction angles of these assemblies under the triaxial drained shear conditions with different sliding frictions are analyzed. The corresponding micro-mechanical behaviors are also explored to explain the observed OR- and sliding frictiondependent macroscopic phenomena. The major conclusions are as follows:

(1) The particle OR and sliding friction jointly affect the stressstrain relationships and internal friction angles of dense sands. The peak friction angle increases with the sliding friction,while the critical friction angle firstly increases with the sliding friction and then stabilizes at a particle ORdependent value. The sliding friction threshold above the critical friction angle increases with the particle OR. On the other hand, the particle OR weakens the peak and critical friction angles given a specific sliding friction. A more representative shape descriptor considering both effects of the OR and sliding friction is recommend when building the relationships between particle properties and macroscopic soil properties in the future studies.

(2) During the shearing, the mean coordination number decreases while the mean normal contact force increases with the particle sliding friction and OR. The mean mobilized friction coefficient increases with the particle sliding friction while decreases with the particle OR, which partly explains the sliding friction- and OR-dependent peak and critical friction angles, but fails to interpret the saturation of the critical friction angles when the sliding friction is larger than an OR-dependent threshold.

(3) From the perspective of microscopic fabric anisotropies, the saturation of critical friction angle with the particle sliding friction results from the decreasing anisotropy degrees of the contact normal or normal contact force fabrics being compensated by the increasing one of the tangential contact force fabric.On the other hand,the monotonously decreasing critical friction angle with the OR is due to the larger total decrease of the anisotropy degrees of contact normal,tangential contact force, normal branch vector fabrics than the increase of normal contact force fabric anisotropy degree.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This study was supported by the National Natural Science Foundation of China(Grant Nos.42077238 and 41941019),and the Guangdong Basic and Applied Basic Research Foundation, China(Grant No. 2020A1515011525). All financial support is gratefully acknowledged.

List of symbols

εa, εvAxial and volumetric strains

φp, φcPeak and critical friction angles

ωs,ωcMass proportions of spheres and clumps

ρ Particle density

ρo,ρos,ρocRepresentative overall regularity of particles, spheres and clumps

μcParticle sliding friction during consolidation

μsParticle sliding friction during shearing

μmMobilized friction coefficient

μmMean mobilized friction coefficient

σ1p,σ1cPeak and residual axial stresses

σ3Confining stress

a*The anisotropy degree of the anisotropy tensor

bn,btProjections of branch vector on the unit contact normal and tangential vectors

deEquivoluminal sphere diameter

de10,de50,de60Equivoluminal sphere diameters corresponding to

10%, 50% and 60% cumulative probability

dmaxMaximum equivoluminal sphere diameter

dThe average particle diameter of the granular assembly

fn,ftNormal and tangential contact forces

fnMean normal contact force

fn,maxMaximum normal contact force

kn,ktNormal and tangential contact stiffnesses

mParticle mass

p,qEffective mean stress and deviatoric stress

pc,qcEffective mean stress and deviatoric stress at the critical state

r* Average equivoluminal sphere radius

vi,ωiThe translational and angular velocities of particles along or around the i th direction

Ac,VcSurface area and volume of clump

A50,S50,C50Aspect ratio,sphericity and convexity corresponding to

their 50% cumulative probabilities

E,F,AElongation, flat and aspect ratios

Fi, MiThe resultant force and torque acting on or around the particle centroid

Ii,Ij,IkThe principal moments of inertia along the i th,j th,and k th principal directions

IInertial number

L0,W0,H0Initial length, width and height of the sample

L,W,HCurrent length,width and height of the deformed sample after shearing

McCritical deviatoric stress ratio

NcTotal contact number

P1,P2,P3First,second and third principal dimensions of the clumps

S,CSphericity and convexity of clump

VchConvex hull volume

ZMean coordination number

ΔtcrDEM critical time step

ΔtDEM computational time step

δ Kronecker delta tensor

σ Stress tensor

σ′Deviatoric stress tensor

AcAnisotropy tensor of contact normal vector

Afn,AftAnisotropy tensors of normal and tangential contact force

vectors

Abn,AbtAnisotropy tensors of normal and tangential branch vectors

ΦcContact normal fabric tensor

ΨbnNormal branch vector fabric tensor

ΨbtTangential branch vector fabric tensor

ΩfnNormal contact force fabric tensor

ΩftTangential contact force fabric tensor