Investigations of the effects of particle morphology on granular material behaviors using a multi-sphere approach

Shiv Prshnth Kumr Kodiherl, Guobin Gong, Lei Fn, Stephen Wilkinson,Chrles K.S. Moy

a Department of Civil Engineering, Xi’an Jiaotong-Liverpool University (XJTLU), Suzhou, 215123, China

b Department of Civil Engineering and Industrial Design, University of Liverpool, Liverpool, UK

c Department of Civil Engineering, University of Wollongong, Dubai, UAE

Keywords:Discrete element method (DEM)Particle morphology Granular materials Triaxial compression Fabric

A B S T R A C T This article studies the influences of particle morphology on the behaviors of granular materials at both macroscopic and microscopic levels based on the discrete element method (DEM). A set of numerical tests under drained triaxial compression was performed by controlling two morphological descriptors,i.e.ratio of the smallest to the largest pebble diameter,ξ,and the maximum pebble-pebble intersection angle, β. These descriptors are vital in generating particle geometry and surface textures. It was found that the stress responses of all assemblies exhibited similar behavior and showed post-peak strainsoftening. The normalized stress ratio and volumetric strains flatten off and tended to reach a steady value after an axial strain of 40%. While the friction angles at peak state varied with different morphological descriptors,the friction angles at critical state showed no significant variation. Moreover,evolution of the average coordination numbers showed a dramatic exponential decay until an axial strain of about 15% after which it stabilized and was unaffected by further increase of axial strain. In addition,stress ratio q/p and strong fabric parameter φsd/φs m were found to follow an approximately linear relationship for each assembly. These findings emphasized the significance of the influences of particle morphology on the macroscopic and microscopic responses of granular materials.

1. Introduction

It is generally accepted that particle morphology is an indispensable factor in affecting mechanical behaviors of granular materials(Zhao and Zhao,2019).Many researchers have attempted to study granular materials at different scales including aspects of shear strength, crushability, dilatancy, shear-induced localization and instability(Sukumaran and Ashmawy,2003;Guo and Su,2007;Tsomokos and Georgiannou,2010).However,it is a challenging task to quantitatively relate the morphological features of soil particles to the mechanical behaviors of granular materials using physical experimentation. Instead, micromechanics based numerical models,such as discrete element method(DEM),have been widely applied across many disciplines for capturing the microscopic insights of granular materials (Radjai et al.,1998; Thornton, 2000;Nicot and Darve,2007;Wang et al.,2012;Guo and Yu,2015;Nicot et al., 2017; Jiang et al., 2018; Shen et al., 2019).

Over the last two decades, developments in X-ray microcomputed tomography (μCT) technology have allowed threedimensional (3D) visualization and characterization capabilities at the micro-scale, enhancing significantly the insights into micromechanical behaviors of granular materials. Nevertheless, the number of scanned particles is usually restricted because of limitations in costs and resolutions of μCT scanning. The DEM framework allows reconstruction of particle surfaces and realistic particle shapes that can be modeled using advanced clump template logics(Ferellec and McDowell,2008;Gao et al.,2012;Zheng and Hryciw,2015).Although the limitation of the number of particles also exists for DEM simulations, the specimens for describing the particle morphologies of DEM are repetitive and can represent the realistic particles to the least, qualitatively.

To date, many numerical simulations have been performed employing irregular particle shapes,such as clumped discs(Maeda et al., 2010; Saint-Cyr et al., 2011), polygons (Tillemans and Herrmann, 1995; Hosseininia, 2012), spheropolygons (Alonso-Marroquin and Wang, 2009) in two dimensions, and polyhedrons(Zhao et al., 2015), sphero-cylinders (Pournin et al., 2005), and superellipsoids (Wellmann et al., 2008; Cleary, 2010) in three dimensions.However,several studies indicate that ellipsoids are able to provide a sufficiently realistic behavior for most granular material applications (Lin and Ng, 1997; Ng, 2004; Bagherzadeh-Khalkhali and Mirghasemi, 2009; Zhou et al., 2013). On the other hand,some researchers argued that grain shape underpins specific examples where having a realistic particle shape is important and should be thoroughly accounted for in their modeling (Shin and Santamarina, 2013; Payan et al., 2016; Alshibli and Cil, 2018; Shi et al., 2018; Xiao et al., 2019). In addition, particle irregularity and roughness are also found to have an important influence on the shear strength and volume change of granular materials (Kozicki et al., 2012, 2014; Araujo et al., 2017). However, the effects of particle morphology on macroscopic and microscopic responses of granular materials are not fully explored. The main focus of this article is on two morphological descriptors, i.e. the ratio of the smallest to the largest pebble diameter, ξ, and the maximum pebble-pebble intersection angle, β. A popular commercial DEM package PFC3D(Itasca Consulting Group, 2018) was adopted to perform the simulations in this context. Considering the broad range of variations in particle morphology and computational limitations, a series of DEM simulations at the macroscopic level was explored initially.Then,the microscopic parameters in terms of the average coordination number,percentage sliding fractions and fabric anisotropy measurements were evaluated. Finally, correlations between macroscopic and microscopic quantities of granular materials were established.

2. Methodology

In this study,the elliptical particles were modeled using a multispherical (MS) method. In three dimensions, the Voronoi diagram of a closed surface may be approximated by calculating a constrained Delaunay tetrahedralization of the ellipsoid mesh to represent the surface (Chew, 1989; George, 2003; Dey and Zhao,2004). Given an ellipsoid surface enclosing a volume, a Delaunay tetrahedralization is initially constructed and then for each tetrahedron,the center and radius of its circumscribed sphere are noted as the balls or pebbles of a clump. This approach can be implemented in the Kubrix automatic mesh generation software package(Taghavi, 2000; Simulation Works Inc., 2009) as an option called BubblePack which automatically generates clump templates for PFC3D.In general,the clump template is a rigid collection of pebbles to create an irregular particle shape by joining and overlapping different pebbles at different coordinates. These clumps can translate and rotate as spheres and obey the equation of motion.In addition, a clump template may hold the surface description that can be used for calculation of the inertial parameters and clump visualization. To use the BubblePack, an ellipsoid surface in the form of an STL file(see Fig.1),and two morphological descriptors,ξ representing the particle geometry and β denoting the surface texture of the particle,should be identified.The angular measure of smoothness (β) varies between 0°and 180°, and the ratio of the smallest to largest pebble diameter(ξ)kept in the clump template varies in 0-1, as suggested by Taghavi (2011). The BubblePack outputs a file containing a surface description of the clump template that is compatible with PFC3Das shown in Fig.2.To optimize these two parameters, Fig. 3 shows the effects of ellipsoid clump templates for different values of ξ and β. The present DEM model can model the total range of ellipsoids as controlled by two morphological descriptors. However, due to the computational limitations and associated costs, we herein focus on the following cases:β=150°with ξ=0.4,0.6 and 0.8;and ξ=0.4 with β=100°,130°and 160°.

Fig.1. A clump representing an ellipsoid comprised of three spheres (pebbles).

Fig. 2. A surface representing an ellipsoid.

3. DEM simulations of drained triaxial compression

In 3D DEM, a triaxial test specimen was created as a cubic box bounded by six rigid walls. A random seed generator was used to mimic the realistic distribution of particles and no gravity was applied. It is presumed that the packing with zero inter-particle friction coefficient during specimen generation may achieve the densest state of assembly(Abbireddy and Clayton,2010;Zhao and Zhao, 2019). Each individual assembly is restricted to a cloud of 2845 ellipsoids which means a total of 25,605 individual rigid spheres (pebbles) at most. This number was chosen based on the previous research concerning the potential effects of particle numbers on the simulation results (e.g. Ng, 2004; Wang and Gutierrez, 2010; Huang et al., 2014). Fig. 4 illustrates the particle size distribution of the ellipsoidal assembly. The equivalent diameter can be calculated aswhere Vcis the volume of the clump. The interactions between clump-clump and clumpwall are assumed to follow a generalized linear contact law (Yan,2009; Abedi and Mirghasemi, 2011; Stahl and Konietzky, 2011;Gu et al., 2014; Yang et al., 2016). The friction coefficient between clump-clump during triaxial compression is set to 0.5,whereas for the clump-wall,it is zero.According to Goldenberg and Goldhirsch(2005),the value of stiffness ratio(i.e. the ratio of shear to normal stiffness, Kn/Ks) for granular material is in the range of 1.0

Fig. 3. Effects of varying β and ξ for ellipsoid clump templates.

Fig. 4. The particle size distribution of assemblies.

4. Results and discussion

4.1. Macroscopic parameters

Previously, DEM assembly containing a minimum of 400 particles was shown to be sufficient for the representative elemental volume (REV) (Guo and Zhao, 2013, 2014; Nguyen et al., 2013, 2014). The stress state of granular materials during shearing can be quantified using a Cauchy stress tensor. In a simple form, it can be described by a second-order tensor, a volumetric average of the dyad comprising the contact force and branch vectors within the REV (see Cundall and Strack, 1979;Christoffersen et al.,1981):

where VRis the volume of the REV; andandrepresent the contact force and branch vectors at each contact point c, respectively. The mean stress p and deviatoric stress q are given respectively by p = σii/3 (summation convention used) and(summation convention used), whereis the deviatoric part of σij, equivalent to σij- pδij, in which δijis the component of the Kronecker delta tensor. However, the quantities related to the stress tensors in this paper,such as principal values of the stress tensor, are calculated from the wall boundaries directly since the triaxial shear conditions imposed will ensure no rotations of the principal stresses.

Fig. 5. Initial configurations of ellipsoidal assemblies for triaxial simulations.

The axial strain(ε1)and the volumetric strain(εv)are calculated from the displacements of the rigid walls (boundaries) as follows:

where H and V represent the height and the volume of the specimen at the current time, respectively; and H0and V0are the corresponding initial values at the onset of triaxial shearing.Similarly,ε2and ε3, strains in x and y directions, respectively, can be calculated in the same manner.

4.2. Evolution of stress-strain

Figs.6 and 7 show the stress ratio(deviatoric stress normalized by mean stress, q/p) and volumetric strains versus axial strain for the specimens with different values of ξ and β, respectively. All specimens were found to exhibit similar behavior, i.e. post-peak strain-softening, which is associated with dilatancy in terms of volumetric strain, as shown in Figs. 6b and 7b. This behavior is expected as it is typical for dense non-cohesive granular materials.Moreover, the normalized deviatoric stress and volumetric strains flatten off and tend to reach a steady value after an axial strain of 40%.This indicates that at this point,all specimens have reached a critical state (Schofield and Wroth, 1968; Gong et al., 2012a).Furthermore, the peak shear strength is found to depend on the morphological descriptors,i.e. the lower the value of β, the higher the shear strength, and the higher the value of ξ, the higher the shear strength. As compared to the rough surface of a particle where multiple contacts exist around the particle,a limited number of contacts may be present in-between smoother particles. It is expected that multiple contacts between rough surfaces of the particles would result in stronger interlocking,and as such,lead to an heightened shear strength within the assembly (Ludewig and Vandewalle, 2012; Shamsi and Mirghasemi, 2012). Further discussion on this in terms of friction angles is provided in the following section.

The evolutions of void ratio(e)versus axial strain for specimens with varied values of β and ξ are shown in Fig. 8. These curves follow a similar trend to the volumetric plots presented in Figs.6b and 7b. From Fig. 8, it can be seen that the particle morphology affects the initial void ratios of the granular assemblies. The initial void ratios increase with an increase of β and decrease with an increase of ξ. These variations could be ascribed to the transformation of particle geometry and surface textures with corresponding changes in morphological descriptors.

4.3. Friction angles

Fig. 6. Evolutions of (a) stress ratio versus axial strain and (b) volumetric strain versus axial strain for the specimens with ξ of 0.4 and different values of β.

Fig. 7. Evolutions of (a) stress ratio versus axial strain and (b) volumetric strain versus axial strain for the specimens with β of 150° and different values of ξ.

Fig. 8. Evolutions of void ratio against axial strain for the specimens of: (a) ξ = 0.4 with different values of β; and (b) β = 150° with different values of ξ.

The effects of β and ξ on the macroscopic responses of granular assemblies at peak and critical states are explored in this section. The initial state corresponds to the end of isotropic compression or beginning of shearing, i.e. at an axial strain of εa= 0%. The peak state is identified where the stress ratio exhibits a peak against axial strain. The critical state is considered as the one in which the granular specimens are continuously deforming under constant stresses and constant volumes.According to Figs. 6 and 7, the critical state is reached after an axial strain of 40%. In triaxial conditions, the internal friction angle, φ, represents the shear strength of granular materials, as described below:

Fig. 9. Effects of β and ξ on the internal friction angle.

Fig.10. Evolutions of the average coordination number for the specimens of: (a) ξ = 0.4 with different values of β; and (b) β = 150° with different values of ξ.

Fig.11. Probability distributions of average coordination number for the specimens of: (a) ξ = 0.4 with different values of β; and (b) β = 150° with different values of ξ.

where σ1and σ3are the major and minor principal stresses,respectively. The influences of β and ξ on the internal friction angle are shown in Fig. 9. The error bars stand for the standard deviation of the friction angles corresponding to the critical state.In detail, friction angles are dependent on morphological descriptors at both peak and critical states. However, at the critical state,such effects are less significant compared to that at the peak state. It is noted that the friction angles at the peak state decrease with an increase of β and increase with an increase of ξ.Furthermore, the deviation of friction angles at peak and critical states appears to have an increasing trend as ξ increases, which suggests an increase in dilation angle. This observation is consistent with the trend between volumetric strain versus axial strain,as presented in Fig. 7b. Similar behavior was noticed with an increase in asphericity under drained triaxial simulations (Zhao and Zhou, 2017).

Fig.12. The contact force chain networks at initial (left), peak (middle) and critical (right) states for β = 100° and β = 160°, respectively.

Fig.13. The contact force chain networks at initial (left), peak (middle) and critical (right) states for ξ = 0.4 and ξ = 0.8, respectively.

4.4. Microscopic parameters

4.4.1. Structural stability

The structural stability of a granular assembly correlates with the average coordination number(Z),which is one of the essential microscopic parameters (Thornton, 2015). Z is defined as the average number of contacts per particle,i.e.Z=2Nc/Np,where Ncis the number of contacts, and Npis the total number of particles in the assembly.The number 2 in the definition of Z indicates that the number of particles corresponding to each contact is two. In a 3D assembly,considering the friction between the particles,assuming there is no sliding at each contact, the number of degrees of freedom for a particle is six,which means that the total number of degrees of freedom in the assembly is 6Np.Similarly,the number of constraints at each contact point is three,thus the total number of constraints in the system is 3Nc.If the assembly is considered to be statically determinate,Z equals 4(i.e.3Nc=6Np).According to Gong(2008),if the Z value is larger than or equal to 4,the assembly can be treated as redundant or stable; otherwise, it is structurally unstable.In this study,irrespective of ξ and β,the Z values are all above 4 which indicates structural stability during the entire simulation.A more robust parameter in terms of redundancy index describing structural stability for particle assemblies was suggested by Gong et al. (2012b, c) and Gong and Zha (2013). The evolutions of Z for the specimens with different β and ξ values are recorded during shearing and are shown in Fig.10,where the average coordination number experiences a significant exponential decay prior to reaching ε1=15%where it stabilizes during further increase in axial strain.The contacts used for calculating Z for clumped spheres refer to the clump-clump contacts, which means that the contacts between the spheres constituting the clump are not counted.

The probability distributions of Z for specimens with different morphological descriptors are presented in Fig. 11. It should be noted that the whole range of axial strain (i.e. εa= 60%) is considered for the evolution of probability distributions of Z.These curves are well-fitted by a Gaussian distribution function and there appears to be a systematic skewness to the distributions.Although there is no large variation between the specimens with β = 100°and β=130°,the distribution shifts to the right both as β increases and as ξ decreases. It can be noticed that the peaks are within the range of 6-8 of Z. Furthermore, an increase in β or decrease in ξ leads to a decrease in Z.

4.4.2. Contact force chain networks

In a granular assembly,the transmission of contact forces occurs through inter-particle contacts. Intuitively, the distribution of contact forces should depend upon the contact distribution and as a result, the contact forces will not be distributed uniformly within the system. In order to understand the internal physics inside an assembly during shearing, a visualization of the contact force networks at various states is helpful (Kodicherla et al., 2019). Figs.12 and 13 show the contact force chain networks corresponding to the initial, peak and critical states with different values of β and ξ,respectively. In both figures, before shearing the specimen, i.e. at the initial stage or end of isotropic compression(see left columns in Figs. 12 and 13), contact forces are more or less distributed uniformly throughout the specimen,because the specimen is under a state of isotropic compression. From Fig. 12a and b, at the peak(middle figure), the strong contact forces (i.e. red-colored spots)(see Section 4.5) are transmitted through the crest and toe of the specimen. There are limited strong contacts represented by redcolored spots. However, they are deemed to be adequate enough to take up the external loads as pillars. Meanwhile, weak contact force chains contributed to the stability of the strong contact force chains to a much less degree. Furthermore, as β increases or ξ decreases,the strong contact force chains were found to be stronger,corresponding to higher peak stresses. These observations are consistent with the trends observed in Figs. 6 and 7, respectively.For the remaining shear process, the contacts gradually disappeared, leading to sparser contact force networks before the critical state was reached.

4.4.3. Proportions of sliding

In the DEM simulations,the Coulomb’s friction law was used to control the sliding of contacts.To explore the proportions of sliding contacts, the sliding is governed by a sliding index Sc= |ft|/（μfn），where ftis the tangential force at the contact, μ is the internal friction coefficient, and fnis the normal force at the contact (Gong and Liu, 2017). The sliding of a contact is assumed to occur when Scis greater than 0.999.Fig.14 shows the evolution of proportions of sliding contacts for the specimens with different values of β and ξ,respectively.Clearly,independent of β and ξ,the sliding contacts increase sharply to the peak and reach a steady value at a higher axial strain.Similar observations were also noted by Gu et al.(2014)and Gong and Liu (2017) for dense granular assemblies under triaxial test conditions. Moreover, it shows that ξ has a noticeable effect on the sliding contacts compared with β. This indicates that the increase of ξ can improve the interlocking between the particles, as evidenced in the observations of the stress ratio in Fig. 7.

4.5. Correlation between macroscopic and microscopic parameters

Rothenburg and Bathurst (1989) pointed out that geometrical and mechanical anisotropies were the main anisotropy sources for governing the shear strength of a granular material. The former refers to the local orientation of a contact plane,producing a global anisotropy, while the latter corresponds to the external forces and depends mainly on the contact forces induced between particles with respect to the local orientation of the contact plane.These two anisotropies affect the overall behavior of granular materials significantly(Guo and Zhao,2013;Zhao and Guo,2013).A secondorder tensor was used to characterize the fabric tensor,which was used to quantify the contact normal distributions of the overall contact network (Oda,1982; Satake,1982):

Fig.14. Proportions of sliding fractions for the specimens of:(a)ξ=0.4 with different values of β; and (b) β = 150° with different values of ξ.

where nkdenotes the unit contact normal vector at contact k with i,j = 1, 2 and 3. A similar second-order fabric tensor was used to describe the inherent fabric anisotropy in terms of vector magnitude based on image analysis approach (Kodicherla et al., 2018).Following Radjai et al.(1998)and Shi and Guo(2018a,b),the entire contact network can be classified as strong and weak subnetworks.Simply,it can be said that if the contact normal force is greater than the average contact normal force,then the contact is said to be in a strong subnetwork; otherwise, it is in a weak subnetwork. The average normal contact force is defined by

where φm= （φ1+φ2+φ3）/3 = 1/3 represents the mean principal fabric.

The principal deviatoric fabric in the strong contact network,is defined by

Figs.15 and 16 show the overall and strong fabric measures of φd/φmandagainst axial strain for the specimens with different values of β and ξ,respectively.Intuitively,both overall and strong fabric measures behave in a similar fashion but higher magnitudes are noted in the strong fabrics. The increase in β or ξ leads to an increase in the fabric measures corresponding to the same axial strain, which indicates that the anisotropy in terms of contact normal orientations increases.

Fig.15. Evolutions of φd/φm against axial strain for the specimens of: (a) ξ = 0.4 with different values of β; and (b) β = 150° with different values of ξ.

Fig.16. Evolutions of against axial strain for the specimens of: (a) ξ = 0.4 with different values of β; and (b) β = 150° with different values of ξ.

Fig.17. Relationship between strong fabric parameter and stress ratio q/p.

Numerous attempts have been made in the literature to relate microscopic quantities with macroscopic parameters (e.g. Radjai et al., 1998; Alonso-Marroquin et al., 2005; Sazzad and Suzuki,2013). Also, from the literature, some relationships between the stress ratio and the strong fabric measure are proposed (Antony,2000; Antony and Kuhn, 2004; Sazzad and Suzuki, 2013; Liu et al., 2018). Following previous investigations, a similar attempt has been made in order to establish a link between the macroscopic and microscopic quantities for different morphological descriptors,as plotted in Fig. 17. The relationships are shown to be similar,indicating that the particle morphology did not alter the trend of the macro-micro relationship. Moreover, it can be inferred from the figure that regardless of β and ξ, the strong correlation is approximately linear before the peaks, i.e. during hardening.

5. Conclusions

A set of 3D DEM simulation tests under drained triaxial compression is performed considering the effects of two morphological descriptors.Initially,the effects of morphological descriptors are examined at the macroscopic level.In addition,the microscopic parameters in terms of the average coordination number, percentage of sliding fractions and fabric anisotropy measures are also investigated.

(1) The stress responses of all assemblies exhibit similar behavior and show a post-peak strain-softening response.The normalized deviatoric stress and volumetric strains flatten off and tend to reach a steady value after an axial strain of 40%where they are considered to have reached the critical state.

(2) Compared to the friction angles at the critical state, the friction angles at the peak state are significantly affected when the morphological descriptors are varied. The difference between the friction angles at the peak and critical states appears to increase with increasing ξ, which corresponds to an increase in dilation angle.

(3) The average coordination number (Z) shows a significant exponential decay prior to reaching εa= 15% and it eventually stabilizes with further increase in axial strain.In all the specimens,the Z values are above 4 indicating a structurally stable system.

(4) Irrespective of the morphological descriptors, an approximately linear relationship is found between the stress ratio and the strong fabric parameterbefore the peak for each assembly.These relationships represent the correlation between macroscopic and microscopic responses of granular materials induced by the changes in morphological descriptors.

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

The authors would like to express their gratitude for the financial support from Xi’an Jiaotong-Liverpool University (XJTLU)(Grant Nos.RDF 15-01-38 and RDF 18-01-23).Also,the support by the Key Program Special Fund at XJTLU (Grant No. KSF-E-19) and Natural Science Foundation of Jiangsu Province (Grant No.BK20160393) is greatly appreciated.