Applicability of discrete element method with spherical and clumped particles for constitutive study of granular materials

Tongming Qu ,Min Wng ,Yuntin Feng Computtionl Engineering,Swnse University,Swnse,Wles

b T-3 Fluid Dynamics and Solid Mechanics Group,Theoretical Division,Los Alamos National Laboratory,Los Alamos,NM,87545,USA

Keywords:Discrete element method(DEM)Granular materials Constitutive behaviour Deviatoric hardening model Rolling resistance model Irregular particles

ABSTRACT Discrete element method(DEM)has been intensively used to study the constitutive behaviour of granular materials.However,to what extent a real granular material can be reproduced by virtual DEM simulations remains unclear.This study attempts to answer this question by comparing DEMsimulations with typical features of experimental granular materials.Three groups of models with spherical and clumped particles are investigated from four perspectives:(i)deviatoric stress and volumetric behaviour;(ii)critical state behaviour;(iii)stress-dilatancy relationship;and(iv)the evolution of principal stress ratio against axial strain.The results demonstrate that DEM with spherical or clumped particles is capable of qualitatively describing macroscopic deviatoric stress responses,volumetric behaviour,and critical state behaviour observed in experiments for granular materials.On the other hand,some qualitative deviations between experiments and the investigated DEM simulations are also observed,in terms of the stress-dilatancy behaviour and principal stress ratio against axial strain,which are proven to be critical for constitutive modelling.The results demonstrate that DEM with spherical or clumped particles may not necessarily fully capture experimental features of granular materials even from a qualitative perspective.It is thus encouraged to thoroughly validate DEM with experiments when developing constitutive models based on DEM observations.

1.Introduction

Granular materials are prevalent in geotechnical engineering.It is not only important but also challenging to predict how granular materials deform when subjected to various external loads.The mechanical complexity of granular media partially arises from its unique characteristics,such as anisotropy and heterogeneity(Ueda and Iai,2019;Petalas et al.,2020),loading rate,loading path and pressure dependence(Gudehus,1996;Li et al.,1999;Das and Das,2019),and strain-localisation within unstable granular materials(Desrues and Andò,2015;Dafalias,2016;Qu et al.,2019a).Although much attention has been paid to the constitutive modelling of granular materials over the past decades,it remains an ongoing challenge to develop a uni fied theoretical model.

Discrete element method(DEM)is a powerful method to perform virtual mechanical testing of granular materials.It re flects the bulk behaviour of granular materials by simulating the particlescale interactions.Despite some idealised simpli fications in terms of particle shape,contact relations and mechanical features of materials,DEM is increasingly popular in the constitutive modelling of granular materials in recent years.Generally,the role of DEM in constitutive modelling can be classi fied as:(i)DEM provides insightful observations for constitutive behaviour of granular materials(Sitharam et al.,2002;Li and Dafalias,2012;Vidyapati and Subramaniam,2012;Yang et al.,2018;Nguyen et al.,2020a,b;Wang et al.,2020,2021a;Zhao and Kruyt,2020),and then is used to improve existing constitutive models;(ii)DEM serves as a validation for newly-developed constitutive models (Sun and Sundaresan,2011;Vidyapati and Subramaniam,2016;Yuan et al.,2019;Hu et al.,2020;Wang et al.,2020,2021a;Xue et al.,2021);(iii)DEM is intended to directly replace the phenomenological constitutive models(Guo and Zhao,2014;Zhao et al.,2020;Qu et al.,2021a;Wang and Zhang,2021)in the framework of hierarchical multiscale modelling of granular materials;and(iv)DEM generates data specimens to train data-driven constitutive models(Wang and Sun,2019;Karapiperis et al.,2021;Qu et al.,2021b).Although DEMhas been extensively applied to various constitutive modelling problems,whether DEM is truly reliable,or to what extent DEM can reasonably replicate the mechanical behaviour of real materials is still not clear.

On one hand,previous work shows the excellent capability of DEM in reproducing complex granular behaviour.Calvetti et al.(2003)showed that DEM models are particularly effective in reproducing complex behaviours of real sand materials across a wide range of loading paths.Huang et al.(2014a)found that DEM simulations agree well with the NorSand constitutive model in terms of macroscopic deviatoric stress ratio and volumetric deformation.Kawamoto et al.(2018)quantitively replicated experimental triaxial testing across both macroscopic stress-dilatancy responses and mesoscopic evolvement of the shear band via LSDEM.Jiang et al.(2020)demonstrated that DEM simulations and the double yielding surface constitutive model can yield satisfactory agreements in predicting the complex mechanical behaviour of granular soils.Recently,Sibille et al.(2021)demonstrated that using DEM with an improved rolling resistance contact law is capable of yielding satisfactory qualitative and quantitative predictions on the deviatoric stress and volumetric strain of both dense and loose soils.

On the other hand,some potential deviations of DEM in simulating real materials are also reported.Fukumoto et al.(2013)found that the presence of rolling resistance quantitively affects both microstructures and macroscopic behaviours of granular materials.The failure of realistic modelling on particle-scale motions will lead to a deviation of the modelled microstructural evolution within materials and then affects the macroscopic volumetric behaviour.Huang et al.(2014b)reported that for sphere-based DEM simulations with a relatively high sliding friction coef ficient(i.e.0.5),the relation between the critical state void ratio and the con fining pressure con flicts with the experimental observations of real sand.Santos et al.(2020)showed that only DEM models that consider rolling and twisting resistances are capable of reproducing the experimental volume fractions.

In this study,we attempt to explore the ability of DEM in predicting constitutive responses of granular materials by performing a series of DEM simulations with spherical and non-spherical particles.Some critical mechanical behaviours of virtual DEM materials are investigated in detail,including macroscopic deviatoric and volumetric behaviours,the theory of critical state soil mechanics(CSSM)(Scho field and Wroth,1968),stress-dilatancy relation,and the evolution of principal stress ratio against axial strain.We further discuss some potential issues arising from the deviation in simulating granular materials when taking DEMinto constitutive research.

2.DEM model preparation

To investigate the applicability of DEMto the constitutive study of granular materials,three groups of different specimens are prepared,as shown in Fig.1.Group A comprises spheres with only sliding friction.Group B consists of spherical particles with both sliding and rolling resistances,which implicitly introduces a torque acting on two contacting particles as a complement to the interlocking effects of irregular particle shapes.Note that several different rolling resistance models can be found in Ai et al.(2011).In this paper,the model proposed by Wensrich and Katterfeld(2012)is adopted.Group C is made of clumped particles with varied sizes and the shape-induced interlocking effects are directly simulated.

Fig.1.Model con figurations:(a)Spherical particle-based specimen(Groups A and B),and(b)Clumped particle-based specimen(Group C).

Groups A and B share the same particle parameters except for the presence of rolling resistance in Group B.Some key particlescale parameters used in simulations are summarised in Table 1.Similar parameters are commonly used to model geomaterials.Some existing studies(Goldenberg and Goldhirsch,2005;Guo and Zhao,2013)recommend that the ratio of normal stiffness(kn)to shear stiffness(ks)should range from 1 to 1.5 for realistic granular materials.Thus,most simulations in this work adopt a value ofkn/ks=1.5,but other values from 1 to 10 are also explored.A total of 4120 particles are involved for all the simulated sphere-based specimens(Groups A and B).For Group C,the clumped particles are generated by permanently fixing two identical spherical particles together and each has a radius ranging from 1 mm to 2 mm with a constant aspect ratio(AR)of 0.5.A total of 5000 clumped particles are generated for Group C.The particle sizes are uniformly distributed in all specimens.Both linear and nonlinear contact models are considered in this work.The servo-wall boundary conditions are adopted in all the simulations.Note that our previous research has shown that the in fluence of boundary conditions on the macroscopic mechanical behaviour is negligible for a triaxial specimen involving around 2300 particles,although the microstructural aspects(the distribution of contact forces and fabrictensor)will be different(Qu et al.,2019a).In this study,we mainly investigate the macroscopic responses of granular materials and thus the effects of boundary conditions are not discussed here.

Table 1Key particle-scale parameters for DEM models.

Table 2Initial void ratios for different specimens after isotropic compression.

Fig.2.Deviatoric stress and volumetric responses for Group A(μ=0.5):(a)Deviatoric stress against axial strain;and(b)Volumetric strain against axial strain.

In the model preparation phase,all particles are initially generated in a large space and then isotropically compressed to a prescribed con fining pressure state via the servo-wall algorithm.In this study,the constitutive behaviours of granular materials are studied under five different con fining stresses.Speci fically,we restrict most investigations to relatively small con fining stresses(200 kPa,400 kPa and 600 kPa),and higher con fining stresses(1000 kPa and 2000 kPa)are simply used to describe the critical state behaviour,thus particle breakage is not considered in our DEM models.Different initial contact friction coef ficients between particles are used to obtain isotropically dense or loose specimens separately.Similar specimen preparation methods can be found in Gong et al.(2019a,b)and Nie et al.(2020).After the specimens have been made,a low friction coef ficient of 0.2 and relatively high friction coef ficient of 0.5 are used to generate two types of specimens with different frictional strengths.Then the specimens are ready for axial loading.Since the mechanical responses of granular materials depend on the state parameters(i.e.the con fining stress and initial void ratio)after isotropic compression(Been and Jefferies,1985),the initial state parameters of each simulated specimen are provided in Table 2.

3.Numerical model results

3.1.Deviatoric stress and volumetric responses

As a specimen deforms incrementally during triaxial testing,the logarithmic strain is adopted to calculate axial strainε33and volumetric strainεv.Also,compressive strains(including volume)and stresses are assumed positive:

whereH,H0,VandV0are the current axial length,the initial axial length,the current volume,and the initial volume of the specimen,respectively.Note that engineering strain,instead of logarithmic strain,is also commonly adopted in laboratory element testing,but it is found that different strain de finitions do not affect the upcoming investigations.A detailed comparative study on the influences of strain formulation can be found in Appendix A.

Figs.2-4 show the macroscopic stress and volumetric responses of dense and loose granular materials.For dense specimens,granular materials show an apparent stress softening phenomenon and dilative deformation after the initial compaction,while the loose specimens show signi ficant stress hardening and contractive deformation during the axial loading.Besides,the volumetric strain is approximately constant when the loading strain is suf ficiently large.

Fig.3.Deviatoric stress and volumetric responses for Group B(μ=0.5,andμr=0.5):(a)Deviatoric stress against axial strain;and(b)Volumetric strain against axial strain.

From a qualitative perspective,all the macroscopic stress-strain features and volumetric behaviour are consistent with the basic characteristics of real granular materials.From a quantitative perspective,the peak and residual deviatoric stresses measured from Group A are signi ficantly lower than those in Groups B and C.It is understandable that spherical particles without rolling resistance lack interlocking among particles,which is one of the particle-scale origins contributing to the macroscopic strength of granular materials.In this case,spherical particles without considering rolling resistance may be insuf ficient to quantitively represent real granular materials with high shear strength.Some existing work uses an arti ficially increasing sliding friction coef ficient to enhance the shear resistance of sphere-based specimens(Qu et al.,2019b).Mollon et al.(2020)showed that this approach is possible,to a certain extent,via the observations of shear localisation simulations.However,whether such a non-physical representation may cause other problems still deserves further investigation.

3.2.Critical state behaviour

The constitutive behaviour of granular materials is highly dependent on their states(void ratio and stress level).A typical feature for granular materials is that a critical state line(CSL)is found when materials are sheared to their critical states(Been and Jefferies,1985;Huang et al.,2014c,d).For sand materials,the CSL can be depicted by(Li et al.,1999):

whereeandp′are the void ratio and the corresponding mean effective stress,respectively;pais the atmospheric pressure(taken as 100 kPa);elim,Λ,andξare the curve-fitted critical state parameters.

Fig.4.Deviatoric stress and volumetric responses for Group C(μ=0.5):(a)Deviatoric stress against axial strain;and(b)Volumetric strain against axial strain.

To gain a thorough understanding of the critical state behaviour of virtual DEM specimens,two more groups of con fining stresses with 1000 kPa and 2000 kPa for both dense and loose granular materials are added to this part.Two different friction coef ficients are used to represent granular materials with different frictional strengths.The initial states,the critical states and the evolution process of state parameters for all the specimens are shown in thee-log10p′space in Figs.5-7.

All the specimens with both low and relatively high friction coef ficients show a very similar trend.Although these specimens have different initial void ratios,the critical states of the two specimens in the same group under the identical con fining pressure have similar mean effective stresses and critical void ratios.With an increase in the mean effective stress,the critical void ratio decreases gradually.A total of 10 critical state points for each specimen can be excellently described by the critical state equation(Eq.(3)).The squared correlation coef ficients exceed 0.99 for most of the specimens.

Fig.5.Evolution of the critical state and CSL in the e-log10 p′space for Group A:(a)μ=0.2;and(b)μ=0.5.

Huang et al.(2014b)reported that for sphere-based DEM simulations with a relatively high sliding friction coef ficient(i.e.0.5),the basic trend between the critical state void ratio and the con fining pressure may con flict with the experimental observations of real sand.However,such an abnormality is not found in our study.Similarly,Barnett et al.(2020)also did not observe unusual results when using a large friction coef ficient.The reasons may be attributed to different ways to calculate the void ratio,i.e.the‘floating’particles which do not transmit contact forces are eliminated from the void ratio calculation in Huang’s work,while we view these‘floating’particles as solid constituents and the volume of‘floating’particles should be considered in the calculation of void ratio.

The results in Figs.5-7 show that the critical state theory is universal for all the virtual DEMsimulations.It is found that,under the same con fining pressure,the critical void ratios of spheres with rolling resistance and clumped non-spherical particles are higher than that of spheres with only sliding friction.Furthermore,the critical void ratio for the specimen with a large friction coef ficient is higher than the counterpart with a smaller friction coef ficient in each group under the identical con fining stress.The results may support that granular materials with a higher shear strength tend to have a larger critical void ratio.On the other hand,except for the differences in quantity,no signi ficant qualitative differences are observed in critical states for the investigated specimens.

3.3.Stress-dilatancy relation

Stress-dilatancy behaviour is another important component for the constitutive modelling of granular materials(Zhang and Salgado,2010;Yin and Chang,2013).Ignoring the contributions of elastic strain for specimens under large strain conditions,Rowe(1962)’s stress-dilatancy equation can be expressed as(Wan and Guo,1998;Sun et al.,2019):

Fig.6.Evolution of the critical state and CSL in the e-log10 p′space for Group B:(a)μ=0.2,andμr=0.2;and(b)μ=0.5,andμr=0.5.

Figs.8-10 demonstrate the stress-dilatancy relations of the simulated specimens.An approximately linear relation betweenand 1-dεv／dε1can be observed in all the investigated loose specimens with both low and high friction coef ficients.This is qualitatively consistent with the understanding that Eq.(4)can be applied to the entire strain history for loose sands(Been and Jefferies,2004)and the experimental results reported in Wan and Guo(1998)and Sun et al.(2019).

Fig.7.Evolution of the critical state and CSL in the e-log10 p′space for Group C:(a)μ=0.2;and(b)μ=0.5.

On the other hand,for the modelling of dense sands,it seems that sphere-based simulations with a relatively low interparticle friction cannot reproduce the“loop”between the pre-and postpeak regimes,as reported in Wan and Guo(1998).In particular,although the“loop”can be found in dense specimens with a large friction coef ficient,the curvature direction is opposite to the one observed from experimental results(Wan and Guo,1998).This abnormal phenomenon may arise from the unreasonable volumetric strain rate that happens in the investigated simulations,although the underlying mechanism responsible for the deviation is not clear yet.

Furthermore,it seems that the stress-dilatancy relations for dense specimens of Groups A and B exhibit slightly pressuredependency.However,the empirical stress-dilatancy relation described by Eq.(4)shows that the relation betweenand 1-dεv／dε1is irrelevant to pressure levels.The experimental results reported in Wan and Guo(1998)and Sun et al.(2019)also demonstrate that the effects of con fining pressure on the stressdilatancy relations are limited.

The King grew as white as a sheet when he heard the Prince s story, and said, Woe19 is me, my son! The time has come when we must part, and with a heavy heart he told the Prince what had happened at the time of his birth

To examine the role of nonlinear contact models in the simulated stress-dilatancy behaviour,we further run the above si mulations using the nonlinear Hertz-Mindilin contact model(Qu et al.,2020).All the particle-scale parameters and procedures of specimen preparation are kept the same as those in Groups A,B and C.It seems that the deviations in stressdilatancy relati on between DEM si mulations and experiments are independent of the type of contact models.For instance,the stress-dilatancy relations for the simulations of Group A with the Hertz-Mindili n contact model are shown in Fig.11.The results show that the reported issues found in Groups A,B and C are irrespective of linear or nonlinear contact models adopted in DEM.

Fig.8.Stress-dilatancy relation for Group A:(a)μ=0.2;and(b)μ=0.5.

3.4.The principal stress ratio against deviatoric strain

The evolution of the principal stress ratio(σ1／σ3)also provides an important perspective to understand the shear behaviour of granular materials.Figs.12-14 show the principal stress ratio against the axial strain for the DEM simulations of Groups A,B and C with the linear contact model.Fig.15 presents the results of DEM simulations for Group A with the Hertz-Mindlin contact model.Similar to the numerical modelling results in Section 3.1,the basic stress hardening or softening behaviour can be well reproduced for both dense and loose materials in the simulations.

Fig.9.Stress-dilatancy relation for Group B:(a)μ=0.2,andμr=0.2;and(b)μ=0.5,andμr=0.5.

However,it is found that the evolutional trend of the principal stress ratio curves under different con fining stresses is slightly different from each other in DEM simulations and thus signi ficant intersection points among them can be found.For dense specimens,the intersection points tend to occur before the peak strength in Group A,while these points are normally observed in the post-peak regime of principal stress ratios in Group B.At the early stage ahead of the intersection points,the relations of principal stress ratios among different con fining pressures are in agreement with experiments,i.e.the principal stress ratio with a lower con fining pressure is higher than that with a higher con fining pressure.This relation is reversed after these intersection points.These observations are inconsistent with most experimental findings of granular materials under a relatively low con fining pressure(＜1 MPa),such as rock fill materials(Xiao et al.,2014;Sun et al.,2019)and river sand(Lade,1977;Kolymbas and Wu,1990;Wan and Guo,1998;Alshibli and Cil,2018).A representative experimental example of this relation for granular materials can be found in Fig.16,which shows that the principal stress ratios under different stress conditions have a similar evolutional trend for experiments and no apparent intersection points are found among different curves.

This abnormal phenomenon may originate from some deviations when the virtual DEM describes the microstructural evolvement of physical granular materials and is potentially related to the improper re flection of particle-scale motion in the current DEM simulations.Note that some phenomenological constitutive models,such as the deviatoric hardening model,use the principal stress ratios to fit the hardening parameters.If DEM cannot properly reproduce similar experimental observations,the phenomenological model fitted from the DEM results will introduce inevitable deviations.The signi ficance of a high-fidelity re flection of the principal stress ratios for constitutive modelling of granular materials is further exempli fied in Section 4.

Fig.10.Stress-dilatancy relation for Group C:(a)μ=0.2;and(b)μ=0.5.

4.Potential issues of DEM for constitutive models

The previous findings show that DEMis capable of qualitatively reproducing many typical mechanical behaviours,such as the critical state of granular materials.However,we also find that some other phenomena,such as stress-dilatancy behaviours and the relation between the principal stress ratio and the deviatoric strain,are inconsistent with experimental observations.It is noted that we also performed extra DEMsimulations withkn/ksranging from 1 to 10,a range used in most existing DEMapplications,but the abovementioned inconsistencies still exist.

In the condition when DEM is used to provide numerical observations to improve a phenomenological constitutive model,the inconsistency between DEM simulations and experimental results in the mechanical behaviour of real granular materials will be reflected in speci fic constitutive models.For instance,it may yield different plastic flow rules or hardening laws.

A phenomenological constitutive model normally requires a yield function,a hardening law and a flow rule to formulate an incremental stress-strain analytical relation.Particularly,thehardening law governs the change in material strength with the plastic deformation and is thus essential for a constitutive equation.In this section,we will illustrate some potential issues using the classical deviatoric hardening model,in which the deviatoric plastic strainεpqis selected as the hardening parameter(Poorooshasb and Pietruszczak,1985;Schweiger et al.,2009;Pietruszczak,2010;Wang et al.,2021b).In the deviatoric hardening model,the yield functionfand the plastic potential functionΦin thep-qstress space are respectively de fined as

where ηcis the slope of the zero-dilatancy line.A hyperbolic hardening law is adopted and described by

whereMis a material constant associated with the hardening law,andηfis the slope of the failure line.

whereKandGare the small-strain bulk and shear moduli,respectively;andHeandHpare the elastic and plastic hardening moduli,respectively,as de fined below:

Fig.13.Principal stress ratio against axial strain for Group B:(a)μ=0.2,andμr=0.2;and(b)μ=0.5,andμr=0.5.

The required parameters for the deviatoric hardening model are determined via DEM-based triaxial testing in Group A.In the deviatoric hardening model,the hardening parameterMis usually fitted from therelation,as shown in Fig.17.It should be highlighted that therelation is an alternative form of the relation between the principal stress ratio and the deviatoric strain discussed in Section 3.4.The former can be transformed from the latter.

The basic mechanical parameters of the dense specimen are listed below:K=115.2 MPa,G=108.9 MPa,ηf=1.13,ηc=1.045,andM=0.001.

These parameters are determined by performing DEM simulations.Kis measured by performing small-strain isotropic compression tests,and calculated byK=dp/dεv,with dpand dεvbeing increments of the mean effective stress and volumetric strain,respectively.Gis determined byG=dq/dεqvia small-strain triaxial tests with constant mean stress,where dqand dεqare increments of deviatoric stress and strain,respectively.ηcis determined by fitting the ratio ofq/punder different con fining pressures at the zero-dilatancy point where the volumetric deformation transits from compaction to dilatancy.ηfis obtained from thep-qevolution under different con fining pressures,as shown in Fig.18.

Fig.19 compares the mean stresspvs.volumetric strainεvrelations obtained from the DEM simulations and the deviatoric hardening model.The relations of deviatoric stress-strain curves of the DEM simulations and the deviatoric hardening model under different con fining pressures are compared in Fig.20.

Fig.14.Principal stress ratio against axial strain for Group C:(a)μ=0.2;and(b)μ=0.5.

We can observe that a satisfactory consistency between the DEM simulations and the deviatoric hardening model is observed from the results shown in Figs.19 and 20.However,if we compare theεq-εvrelation obtained from the DEM simulations with the classical deviatoric hardening model,a signi ficant deviation can be found in Fig.21.In fact,the issue presented in Fig.21 can also be observed from experimental results.This issue is a drawback of the classical deviatoric hardening model and arises from the original hardening law described by Eq.(5).A hyperbolic hardening law with a pressure-dependent hardening parameterMcan solve the problem of the classical deviatoric hardening model.

5.Conclusions

Fig.15.Principal stress ratio against axial strain for DEM simulations of Group A with the Hertz-Mindilin contact model:(a)μ=0.2;and(b)μ=0.5.

Three different types of DEM specimens(spheres with only sliding friction,spheres with both sliding and rolling resistances,and clumped particles)are prepared for triaxial simulations of dense and loose granular materials.The microscopic parameters in this work are selected from a commonly used parameter space in the existing literature.Four important aspects pertaining to the constitutive modelling of granular materials are examined.The results show that basic mechanical behaviour and critical state phenomenon of real granular materials can be properly reproduced by DEM,while some inconsistencies exist.Particularly,these inconsistencies may become critical when using DEM to improve phenomenological constitutive models.The detailed findings are:

(1)Virtual DEM simulations with either spheres or irregular particles represented by clumped particles are capable of properly reproducing stress hardening/softening,volumetric shrinkage/expansion,and critical state behaviour of granular materials.

(2)The stress-dilatancy relations obtained from DEM simulations with spheres and irregular particles for dense specimens may not be a reasonable re flection of real granular materials,and this observation is also found for the simulations with the nonlinear contact model.

(3)The evolution of the principal stress ratio against axial strain under different con fining stresses seems inconsistent with experimental results.Moreover,the deviation seems to be independent of particle shape and contact models.

Fig.16.A representative experimental example of principal stress ratio against axial strain for granular materials(Kolymbas and Wu,1990):(a)Dense sand;and(b)Loose sand.

Fig.17.Determination of material constant M for a dense DEM specimen.

Fig.18.p-q relation(ηf=1.13).

Fig.19.p-εv relation.

Fig.21.εq-εv relation.

Fig.22.Pressure-dependent hardening(experiment results)(Kolymbas and Wu,1990).

These results show that although DEMis capable of reproducing the basic mechanical behaviour of granular materials from the simulation of particle-scale interactions,it may not be necessary to fully capture the experimental features of granular materials,even from a qualitative perspective.The causes of the reported issues and associated physics will be further explored.In this work,only clumped particles with two spheres fixed together are used to represent non-spherical particles.Particles with truly arbitrary shapes will be used to investigate the simulated constitutive behaviour of granular materials in our next work.Also,this study explores the performance of DEM simulations mainly from a qualitative perspective.More quantitative comparisons with varied particle-scale parameters will be examined as well.

Fig.20.εq-q relation.

Declarationofcompetinginterest

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

AppendixA.Supplementarydata

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.09.015.