Particuology 56 (2021) 152–162

Contents lists available at ScienceDirect

Particuology
journal homepage: [Link]/locate/partic

Discrete element modelling of strength and critical state

characteristics of granular materials under axial compression and
axial extension stress path tests
Shiva Prashanth Kumar Kodicherla a , Guobin Gong a,∗ , Lei Fan a , Stephen Wilkinson b ,
Charles K.S. Moy a
a
Department of Civil Engineering, Xi’an Jiaotong – Liverpool University (XJTLU), China
b
Department of Civil Engineering, University of Wollongong, Dubai

a r t i c l e i n f o a b s t r a c t

Article history: The critical state soil mechanics (CSSM) framework has been widely used across a range of problems in
Received 10 March 2020 geomechanics involving complex loading conditions. However, the uniqueness of the critical state has
Received in revised form 2 November 2020 been disputed for many years and it remains a controversial issue. Motivated by previous investigations,
Accepted 10 November 2020
a series of discrete element method (DEM) simulations were performed under both axial compression
Available online 8 December 2020
(AC) and axial extension (AE) stress paths. All samples were isotropically compressed at varying mean
normal effective stresses (confining pressures) and sheared to a large axial strain of approximately 60%.
Keywords:
It is found that there exist unique values of critical void ratios and stress ratios under critical state, which
DEM
Clumped particles
are independent of the samples’ initial packings but dependent on stress paths. And the critical strength
Critical state (stress ratio) for the AC stress path tests is higher than that for the AE stress path. The critical state lines
Densification index (CSLs) are found to path-dependent but unique for each stress path. A unique linear relationship between
State parameter the critical coordination numbers and critical void ratios is identified under the AC and AE stress paths
respectively, but such a relationship depends on the stress paths. It is also found that there exist unique
values of microscopic parameters in terms of deviator fabric under critical state, which are independent
of the samples’ initial packings but dependent on stress paths. All these simulation results lead to the
conclusion of non-uniqueness of CSLs from both macroscopic and microscopic viewpoints.
© 2020 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of
Sciences. Published by Elsevier B.V. All rights reserved.

Introduction is non-unique (Been, Jefferies, & Hachey, 1991; Verdugo & Ishihara,
1996), while others show a unique CSL regardless of the initial
The critical state soil mechanics (CSSM) framework was orig- conditions (Ng, 2009; Salvatore, Modoni, Ando, Albano, & Viggiani,
inally pioneered by Roscoe, Schofield, and Wroth (1958) and 2017; Sitharam & Vinod, 2009). In addition, Yang and Wu (2016)
Schofield and Wroth (1968). The key concept of the CSSM frame- found that CSL is found to be independent of the initial fabrics and
work suggests that the soil is continuously sheared to a large shearing modes.
strain where both the constant volume and constant stress states A critical state void ratio (ec ) alone cannot well represent the
occur. Although the CSSM framework was based on experimental comprehensive state of the soil sample. Been and Jefferies (1985)
observations and theoretical derivations, the conflicting findings introduced a state parameter = e − ec , where e is the current void

regarding some key issues, such as the uniqueness of the critical ratio and ec is the void ratio at the critical state for a given value of p .
state line (CSL) and achieving the critical state, have been discussed A positive value of signifies a looser state while the negative value
for a long time (Chu, 1995; Mooney, Finno, & Viggiani, 1998; Zhao of indicates a denser state of the assembly. Recently, numerous
& Guo, 2013). Some experimental investigations indicated that the constitutive models for soils have been developed which explicitly
CSL in the void ratio (e) and mean normal effective stress (p ) space take into account in their formulation (Jefferies & Shuttle, 2002;
Jefferies, 1993). However, there is a lack of large-strain labora-
tory test data for three-dimensional generalized loading conditions,
which bounds the ability to develop and to apply these models
∗ Corresponding author.
to field applications where the stress state is not axisymmetric
E-mail address: [Link]@[Link] (G. Gong).

[Link]
1674-2001/© 2020 Chinese Society of Particuology and Institute of Process Engineering, Chinese Academy of Sciences. Published by Elsevier B.V. All rights reserved.
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

Fig. 1. Schematic representation of the stress-path models (after Springman et al., 2013).

(Huang, O’Sullivan, Hanley, & Kwok, 2014; Xie, Yang, Barreto, & Overview of DEM and simulation parameters
Jiang, 2017).
Several researchers have worked on the discrete element The advent of DEM by Cundall and Strack (1979) has enabled
method (DEM) for analysing the mechanical behaviour of granular geotechnical researchers to explore the grain-scale interactions
materials under different loading conditions, together with identi- considering various loading paths that may be difficult to achieve
fying their critical-state behaviour (2009, Fu & Dafalias, 2015; Ng, and analyse via laboratory tests. The properties of DEM assem-
2004; Xie et al., 2017; Yang & Wu, 2016). However, the macro- blies are updated at each iteration timestep. The displacements
scopic behaviour of granular materials is not only influenced by the (i.e., translational and rotational) of each individual particle are
loading conditions but also depends on the stress paths (Gutierrez, found by integrating Newton’s second law of motion, whereas the
Wang, & Yoshimine, 2009; Jiang, Li, & Yang, 2013; Rodriguez & Lade, contact forces between the particles are updated using a force-
2013). A more generalised three-dimensional stress state can be displacement contact law (Cundall & Strack, 1979; Hertz, 1882;
achieved by a true triaxial apparatus (TTA) (Green & Bishop, 1969; Itasca, 2003; Mindlin & Deresiewicz, 1953). The simulations pre-
Head, 1986; Ko & Scott, 1967). Moreover, each stress path has its sented in this paper were executed by employing a windows-based
own critical field significance and thus axial compression (AC) and program called particle flow code (PFC3D ), developed by Itasca Con-
axial extension (AE) stress paths were adopted in this work for sulting Group (Itasca Consulting Group, 2019), on a Dell Precision
simplicity (see Fig. 1). T7500 workstation, which has the following specifications: Intel®
Xeon® CPU X5690 @3.47 GHz; 12 cores; 48 GB of RAM; NVIDIA
Quadro 4000.
As pointed out by many researchers, particle shape has an
• Axial compression (AC) – The specimen is loaded axially (a ↑)
important role in affecting the physical and mechanical behaviours
while the confining pressure kept constant (r = 0) (see
of granular materials. Extensive research has been performed by
Fig. 1(a)).
employing wide ranges of non-spherical particles such as sphe-
• Axial extension (AE) – The specimen is unloaded axially (a ↓)
rocylinders (Pournin et al., 2005), ellipsoids (Ng, 2009), clumped
while the confining pressure kept constant (r = 0) (see
spheres (Härtl & Ooi, 2011; Alizadeh Behjani, Hassanpour, Pasha,
Fig. 1(b)).
Ghadiri, & Bayly, 2017; Alizadeh Behjani, Asachi, Ghadiri, Bayly,
& Hassanpour, 2018; Asachi, Alizadeh Behjani, Nourafkan, &
Hassanpour, 2020), superellipsoids (Cleary, 2010; Wellmann, Lillie,
The core idea of this paper is to investigate the uniqueness of & Wriggers, 2008) and polyhedrons (Zhao, Zhou, & Liu, 2015). How-
critical state lines (CSLs) amongst prior observations with respect ever, there is no standard particle shape which can offer a general
to AC and AE stress paths. To achieve this aim, a series of DEM framework to describe arbitrary particle geometry. In this inves-
simulations were performed under drained AC and AE stress path tigation, we made an attempt to mimic the realistic behaviour of
tests, where the samples were sheared to a large axial strain granular materials considering intrinsically chosen sand particle.
to examine the uniqueness of CSLs. In addition, the relation- In general, one can adopt a clump logic to generate irregular par-
ships between and strength, between and dilatancy, and ticle shapes in PFC. A random sand particle in the form of an STL
the correlation between and some microscopic measures were file (clump template) was imported, which was modelled using the
developed. multi-sphere approach given by Taghavi (2000). In clump template,

153
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

each constituent particle is formed by joining and overlapping

individual spheres (called pebbles in PFC) at different coordinates,
which then behave as a rigid body that never break apart, regardless
of the magnitude of forces acting upon it (Li & Yu, 2010; Yang, Yang,
& Wang, 2013). In PFC, modelling the clump template pebbles can
occur via the BubblePack algorithm (Taghavi, 2000). A triangulated
surface must be specified with geometry, and no manual generation
of pebbles is allowed. These pebbles are filtered based on their mor-
phological descriptors such as  which is the ratio of the smallest to
the largest pebble within the clump (0 <  < 1) and ˇ correspond-
ing to an angular measure of smoothness (0◦ < ˇ < 180◦ ). Taking
computational limitations and associated costs into account, the
following parameters are adopted:  = 0.5 and ˇ = 150◦ . The main
effects of the morphological descriptors (i.e.,  and ˇ) on the macro
and microscopic behaviour of granular materials are described by
Kodicherla, Gong, Fan, Wilkinson, and Moy (2020). Fig. 2 illustrates
the distribution
 of the equivalent diameter of the particle, where
Fig. 2. The grain size distribution of the currently used DEM assembly.
deq = 6V/, in which V is the volume of the clump. The assem-
bly was generated using a commonly adopted random generation

Fig. 3. Map of clumps along the different plane of orientations. (a) XYZ, (b) XZY, (c) YXZ, (d) YZX, (e) ZXY and (f) ZYX.

154
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

approach described by Itasca (2003). In this approach, the random Table 1

Input parameters for DEM simulations.
seed can generate the particle, and the location for each individual
particle is identified within the problem domain. In case the par- Parameter description Value
ticle does not fit in, the radius is retained but another location is Notional particle density,  (kg/m3 ) 2700
randomly identified. A map of clumps with the different plane of Inter-particle friction coefficient 0.5
orientations is presented in Fig. 3. Wall-particle friction coefficient 0.0
A cubic sample of dimensions 5 × 5 × 5 cm was used with each Wall stiffness (N/m) 1.0 × 109
Effective modulus, E ∗ (Pa) 1.0 × 108
individual assembly comprised of 2,845 clumps (each clump con-
Normal to shear stiffness ratio, k∗ (≡ kn /ks ) 4/3
sisting of 27 pebbles), equivalent to 76,815 pebbles for the whole Local damping constant 0.7
assembly. This number was selected based on prior investiga-
tions which raised potential conflicting effects on the simulation
results (Ng, 2004; Wang & Gutierrez, 2010). A simple linear
force-displacement contact law was implemented to mimic the
clump-clump and clump-wall interactions (Abedi & Mirghasemi,
2011; Gong & Liu, 2017; Gu, Huang, & Qian, 2014; Kodicherla, Gong,
Yang et al., 2019; Stahl & Konietzky, 2011; Yan, 2009; Yang, Wang,
& Cheng, 2016). The normal contact stiffness of each individual
contact kn is calculated using kn = r 2 E ∗ / (ra + rb ), where E ∗ is the
effective modulus, ra and rb are the radii of constituting clump in
contact, and r represents the minimum value of ra and rb . For gran-
ular materials, the typical range of stiffness ratios (i.e., the ratio
of normal to shear stiffness, kn /ks ) is between 1.0 < kn /ks < 1.5
(Goldenberg & Goldhirsch, 2005). Thus, we adopted the stiffness
ratio of 4/3 (Gong & Liu, 2017; Kodicherla, Gong, Yang et al., 2019,
2020). The input parameters used in the DEM simulations are pre-
sented in Table 1.
During the isotropic compression, an in-built stress-controlled
servo mechanism was triggered in order to achieve the desired Fig. 4. Evaluation of ID .
stress levels. The numerical samples were considered to be equili-
brated until the tolerance level of stress obtained from the walls
are below 0.001. During shearing, the walls were permitted to generation phase, and then e obtained at the end of the isotropic
move under an initial specified loading rate of 0.004 m/s (also with compression phase will define emin . Similarly, for the loosest state,
servo-control), which is small enough to ensure the existence of the interparticle friction coefficient in the successive isotropic com-
quasi-static shear (Zhao & Zhao, 2019). Moreover, the assemblies pression phase shall be gradually increased until e becomes almost
were monitored by the inertia number (I), which is given as (da stable even though the change is imposed in the interparticle fric-
Cruz, Emam, Prochnow, Roux, & Chevoir, 2005): tion coefficient. e obtained at this particular isotropic compression
.  
phase will define emax . At the end, it should be mentioned that the Dr
I = εdi g /p (1) and ID are somewhat different as ID relies on the confining pressure
 
.  p , and consequently, the same e can vary with p .
where ε is the axial strain rate, d is the particle diameter, p is
the mean normal effective stress and g is the particle density. It
should be mentioned that the quasi-static shear can be achieved by Simulation program
a small value of I, thus, I was maintained at a value less than 10–4
throughout the simulations. In total, 24 numerical simulations were performed for 6 differ-
ent series of numerical simulations. In each series, 4 simulations
Evaluation of ICL and ID are performed with different confining pressures (100, 200, 400,
and 800 kPa). Series I–III involve AC stress path tests on three dif-
The state of soil can be well represented by the relative den- ferent states of assemblies, i.e., AC L, AC M and AC D. These states
sity Dr = (emax − e)/(emax − emin ), where emax and emin are the of assemblies typically represent loose, medium and dense states,
maximum and minimum void ratios, respectively. However, the respectively based on their ID . Likewise, the Series IV–VI are pre-
evaluation of emax and emin is challenging in DEM simulations. Thus, pared for the AE stress path. All samples are sheared to a maximum
a substitute measure of Dr , is the densification index (ID ), which is εa of 60%, at which the critical state failure can be clearly identified.
used to capture the state of numerical assembly. The expression for Table 2 shows a summary of the initial characteristics of samples,
ID is defined as (Zhang, Lo, Rahman, & Yan, 2018): where the subscript 0 denotes the onset shear or end of isotropic
eup − e compression.
ID = (2)
eup − elow
Results and discussions
where eup and elow are the corresponding void ratios of upper
and lower consolidation lines, respectively. The upper boundary
Macroscopic evolution
of the isotropic consolidation line (ICL) represents the loosest sam-
ple whereas the lower boundary of the ICL represents the sample
The macroscopic behaviour of granular materials while shear-
exposed to densest state or maximum densification (see Fig. 4).
ing can be evaluated using a Cauchy stress tensor (Christoffersen,
In fact, in DEM, the state of the sample mainly depends on the
Mehrabadi, & Nemat-Nasser, 1981; Cundall & Strack, 1979):
application of interparticle friction coefficient during the specimen
generation or isotropic compression phase. In order to achieve the 1 c c
densest state of the assembly, a zero interparticle friction coef- ij = fi bj (3)
V
ficient should be assigned to the particles during the specimen cεV

155
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

Table 2
Summary of initial characteristics of numerical assemblies.

Test series Simulation ID e0 ID0 0 Zo

AC L 100 0.629 0.109 0.124 5.56

AC L 200 0.608 0.143 0.112 7.09
Series I
AC L 400 0.567 0.204 0.097 9.75
AC L 800 0.527 0.206 0.085 13.5
AC M 100 0.460 0.504 –0.045 8.06
AC M 200 0.451 0.506 –0.045 9.55
Series II
AC M 400 0.435 0.514 –0.035 11.9
AC M 800 0.406 0.505 –0.036 15.8
AC D 100 0.362 0.733 –0.143 10.7
AC D 200 0.353 0.735 –0.143 12.3
Series III
AC D 400 0.337 0.741 –0.133 14.8
AC D 800 0.311 0.741 –0.131 18.9
AE L 100 0.629 0.109 0.109 5.56
AE L 200 0.608 0.143 0.093 7.09
Series IV
AE L 400 0.567 0.204 0.077 9.75
AE L 800 0.527 0.206 0.053 13.5
AE M 100 0.460 0.504 –0.058 8.06
AE M 200 0.451 0.506 –0.061 9.55
Series V
AE M 400 0.435 0.514 –0.070 11.9
AE M 800 0.406 0.505 –0.082 15.8
AE D 100 0.362 0.733 –0.158 10.7
AE D 200 0.353 0.735 –0.161 12.3
Series VI
AE D 400 0.337 0.741 –0.168 14.8
AE D 800 0.311 0.741 –0.179 18.9

where V is the total volume of the assembly, f c and bc are the con-
tact force and branch vector at each contact point c, respectively.
The deviator stress q and mean normal effective stress 
 p are given
   
by p = ii /3 (summation convention used) and q = 3ij ij /2 (ij

is the deviatoric part of ij , equal to ij − p ıij where ıij is the Kro-
necker delta) (summation convention is used). The axial strain ε1
and volumetric strain εv can be obtained from the displacements
of the rigid walls (boundaries), as:

H H 
dh 0
ε1 = − = ln (4)
h H
H0

V V 
dv
εv = − (ε1 + ε2 + ε3 ) = = ln (5)
v V0
V0

where H and V are the height and volume of the specimen at the
current time, respectively, H0 and V0 are their initial values at the
Fig. 5. Drained AC stress path results: (a) εa − , (b) εa − εv , (c) εa − e.
onset of shear. Similarly, the strains ε2 and ε3 in x and y directions
can be respectively obtained in a similar fashion.
The macroscopic evolution of stress ratio () against axial strain Figs. 5(c) and 6(c) present the evolution of void ratio (e) against

(εa ) for the two extreme values of confining pressures (p0 ) and three εa for the two different stress path tests considered in this investi-
different states of assemblies under AC and AE stress path tests are gation, respectively, which is similar to the evolution of volumetric
presented in Figs. 5(a) and 6(a), respectively. It can be seen that strain. e increases with increasing εa in dense samples in general,

independent of the applied p0 , dense samples exhibit post-peak whereas e decreases with increasing εa in loose samples. For both
strain-softening behaviour, whereas loose samples show a strain the cases, e flattens off after εa ∼ 50% and reaches a critical value. A

hardening response. In addition, the medium dense samples exhibit unique value of ec was achieved for a given confining pressure p0 .

initial hardening with minor post-peak strain-softening. Indepen- Moreover, the ec decreases with increasing p0 . It can be found that
dent of stress paths, in all samples, the critical states were attained for a given confining pressure, the value of critical void ratio under
after εa > ∼50%. the AE stress path is higher than that under the AC stress path.
In terms of volumetric responses, it was observed that dense More discussion regarding the ec and the critical state parameters
samples were more dilative while loose samples exhibited extreme are presented in section 5.2.
contractive responses. Initial contraction followed by dilation was Fig. 7 illustrates the comparison of stress path tests results for
found for dense samples in general (see Figs. 5 and 6). These medium dense samples with different confining pressures (100,
behaviours are irrespective of the stress path tests considered dur- 200, 400, and 800 kPa) for Series II and Series V, respectively. For
ing this study. Furthermore, these characteristics well represent both AC and AE stress paths, the q exhibits minor post-peak strain-
the typical behaviour of granular materials documented in the lit- softening behaviour and reaches a critical value after εa > ∼50% (see
erature (Gong & Zha, 2013; Gong, Zha, & Wei, 2012; Guo & Zhao, Fig. 7(a)). The volumetric behaviour for the AC stress path showed
2013). considerable initial contraction (i.e., up to εa ∼ 10%) whereas, for

156
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

Evolution of critical state characteristics

The framework of critical state theory assumes existence of a

unique critical state line (CSL) for a given soil sample in e − q − p
space. The soil sample under shear will be expected to reach crit-
ical state at large strains where the stress state and volume of the
sample both keep steady (while shearing strain will continue to
increase). The CSL can be viewed in a 2D version as the locus of void
ratio – mean effective stress space (i.e. e − log p ) and in q − p
space. As mentioned by Roscoe and Burland (1968) and Schofield
and Wroth (1968), the CSL was originally obtained from the exper-
 
imental observations on remoulded clays and it was found to be
a straight line in the e − log p space. However, the behaviour of
sands is somehow different from clays. After identifying the compli-
cations which arise from the non-linearity, a linearization approach
was proposed by Li and Wang (1998), which is given as:
  ˛
ec = − p /pa (6)

where pa = atmospheric pressure taken as the standard pressure


for normalization (∼ 101.325 kPa); = intercept of CSL at p /pa =
0 axis; = slope of the CSL; and ˛ = material parameter constant
which is taken as 0.7, obtained for a standard Toyoura sand.
Fig. 8 presents the critical state characteristics for AC and AE
stress path tests. From Fig. 8(a), it can be seen that the CSLs show
  paths. The CSL of AE is located
an obvious dependency on the stress
above that of AC in the e − log p space, which is in agreement
with the simulated results performed by Li (2006). In addition,
these CSLs collapses into two different trend lines with different
critical state parameters (see Fig. 8(b)). The AC stress path test
has c = 0.524 and c = 0.013, while the AE path has e = 0.528
and e = 0.01. Both the CSLs fit very well with significant spac-
ing between them and each has a regression value greater than
R2 = 0.99. As is shown in Fig. 8(c), the critical state data points in

the q − p space were linearly fitted and considered to pass through
the origin with slopes of Mc = 1.092 and Me = 0.837, respectively.
Moreover, the critical strength of granular materials in the AC stress
path test is found to be higher than that obtained from the AE stress
path. These findings are in agreement with a study conducted by
Yang and Wu (2016) for undrained loading paths using clumped
particles with an aspect ratio of 0.6.

Relationship between the strength and state parameter

Many investigators (Castro, 1969; Lade, 1972; Norwegian

Geotechnical Institute (NGI) (1982); Jefferies & Been, 2006) have
Fig. 6. Drained AE stress path results: (a) εa − , (b) εa − εv , (c) εa − e. made quantitative comparison between the strength and the state-
dependent response of granular materials. Fig. 9 compares the peak
angle of shearing resistance ( p ) against the initial state parameter
the AE stress path, a well-represented dilation was noticed from ( 0 ). The DEM results of the AC stress path results are generally in
the onset of shear (see Fig. 7(b)). It was found that increase of agreement with those of experimental triaxial compression tests
the confining stresses leads to a lower value of eo , which is true performed on sands. However, the AE stress path simulation data
of both AC and AE stress conditions. Moreover, the critical state are considerably below the experimental values of sands. The p of

stress-ratio for the AC stress path test (Mc = q/p = 1.08) is higher dense samples is an outlier to the right, as it is difficult to achieve

than that obtained for the AE stress path (Me = q/p = 0.83) (see the densest state in experimental tests. The experimental study by
Fig. 7(c)). The ratio of Mc /Me = 1.30 is very close to the typical Cho, Dodds, and Santamarina (2006)) and a DEM study by Maeda,
range of the standard Toyoura sands evaluated from laboratory Fukuma, and Nukudani (2009)) confirmed that the shear strength
experiments (Kulhawy & Mayne, 1990; Lade & Duncan, 1975; of granular materials depends mainly on particle geometry, and
Yang, Li, & Yang, 2008; Yoshimine, Ishihara, & Vargas, 1998). In the numerical study by Huang, O’Sullivan et al. (2014) found that
a numerical study, Ng (2004) found a relative value of stress ratio the shear strength of spherical particles is considerably below the
of 1.12 and 1.09 for two different aspect ratios of ellipsoids (i.e., experimental values, attributing to the particle geometry and inter-
1.2 and 1.5) respectively, which are slightly smaller than those locking.
obtained in this investigation. This may be attributed to the fact Following the approach suggested by Been and Jefferies (1985),
that particle morphological features can affect the critical state the stress dilatancy angle of shearing resistance ( p − c ) against
stress ratios of granular materials (Kodicherla, Gong, Fan, & Moy, the 0 are shown in Fig. 10. The DEM simulations under both the
2019). stress path tests were compared with the experimental datasets

157
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

Fig. 7. Comparison of triaxial stress path results for Series II and Series V: (a) q − εa , (b) εv − εa , (c)  − εa .

documented in the literature (see Fig. 10). Using this approach the the number of contacts and N is the total number of particles in
responses observed in the DEM simulations match well with those the assembly. The numeric number (i.e., 2) explains that each indi-
in the literature (Castro, 1969; Lade, 1972; NGI, 1982; Jefferies & vidual contact is pooled by two contact entities. For multi-sphere
Been, 2006) and as discussed by Huang, O’Sullivan et al. (2014), the clumped particles, the number of contacts refers to the contacts
relationship between p − c and 0 is not shape-dependent. between clump-to-clump, which means that the contacts between
the pebbles constituting the clump are not counted for evaluating
Relationship between dilatancy and state parameter Z. The evolution of Z for Series II and Series V is presented in Fig. 13.
From Figs. 13(a) and (b), it can be seen that for the initial part of
Rowe (1962) mentioned that the similarity between friction shearing, i.e. εa ∼ <10%, the Z increases for the AC stress path tests
angle ( ) and dilation rate (D) may be predicted theoretically while Z decreases for the AE stress path tests, which is due to the fact
and the difference between peak state friction angle ( p ) and that contraction induces new contact formation whereas dilation
critical state friction angle ( c ) is directly related to the dilation induces more contact disruptions. These behaviours are also well
rate. In practice, it is preferable to relate individual behaviours represented by the volumetric responses for both the stress paths
to the fundamental parameters of the state instead of relating considered (see Fig. 7(b)). In addition, for all the samples, Z reaches

derived behaviours to each other. According to Taylor (1948), the a critical value at large εa . Fig. 13(c) shows the relationship of p − Z
strength of soil can be decayed into a dilational and intrinsic compo- at the critical state, which shows that Z increases with increasing

nents. Consequently, the developments of this decomposition led p0 , reflecting a decrease in e. The data presented in Fig. 13(d) indi-
to the stress-dilatancy theory suggested by Rowe (1962). Been and cates that these relationships are more or less linear and unique
Jefferies (1985) introduced that the dilation rate (D = −dεv /dεa ) at but dependent on the stress paths. Furthermore, it should be men-
the peak (Dp ) was related to the 0 . Fig. 11 shows the relationship tioned that Z at critical state for the AC stress path is higher than

between Dp and 0 for both the stress path tests. In addition, for a that of the AE stress path for a given p0 , signifying a higher critical
better comparison, the experimental data from Jefferies and Been stress ratio for AC than for AE in relation to a given sample with a

(2006) is also superimposed in the plot. A unique relationship was given p0 .
found between Dp and 0 , which is path-dependent and the DEM According to Oda (1982), the fabric denotes the spatial arrange-
data lie within the range of the experimental data. ment of particles and associated voids. During shearing, the spatial
Fig. 12 illustrates the relationships between the strength in arrangements of particles and void spaces tend to exhibit an extent
terms of p and p and Dp , respectively. In Figs. 12(a) and (b), a of anisotropy and to evolve to a specifically preferred orientation
linear relationship between p and Dp well represents their corre- (Yimsiri & Soga, 2011). Among the various definitions available for
lation for standard sands (Vaid & Sasitharan, 1991). Based on the the fabric tensor (Li & Li, 2009; Oda, 1982; Satake, 1982), a contact-
DEM simulations results, it is found that the slope and the position normal based proposition introduced by Satake (1982) and Oda
of best-fit lines appear to be path-dependent. Similarly, Ng (2004) (1982) was used in this study, which can be expressed as:
and Huang, Hanley, O’Sullivan, Kwok, and Wadee (2014)) found the ⎡ ⎤
dependency of strength parameters on intermediate stress ratio in 11 12 13
1 
Nc
their numerical simulations in general. ⎢ ⎥
ij = nki nkj = ⎣ 21 22 23 ⎦ (7)
Nc
k=1
Evolution of microscopic parameters 31 32 33

The coordination number, Z (= 2C/N) is a key index to evaluate where nk = the contact unit normal vector of the contact k with
the internal structural characteristics of the assembly, where C is i, j = 1, 2, 3 and Nc is the number of contacts in the assembly.

158
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

Fig. 9. Relationship between p and 0.

Fig. 10. Relationship between p − c and 0.


Fig. 8. Critical state characteristics of AC and AE stress paths: (a) e − p , (b) ec −
0.7 
(p /pa ) , (c) q − p .

The components 11 , 22 and 33 are the fabric in the princi-

pal directions (i = j), and 12 , 13 , 21 , 23 , 31 and 32 are the
fabric components in shear direction (i = / j). The scalar values of
the shear components in the opposite direction are equal, that is, Fig. 11. Relationship between Dp and 0.

12 = 21 , 13 = 31 and 23 = 32 . The principal values 1 , 2 and

3 along the principal directions of fabric can be evaluated using Figs. 14(a) and (b) show the evolution of the deviator fabrics for
the eigenvalues and eigenvectors of the fabric tensor. The deviator the overall and strong contact networks respectively, for medium
fabric, d (= 1 − 3 ), is the difference between major and minor dense samples under the AC and AE stress path tests. For both
principal fabric components, which is widely used to quantify the overall and the strong deviator fabrics, it is found that the
the structural anisotropy of granular assemblies (Thornton, 2000). structural anisotropy increases to the peak and stabilized as εa fur-
A two-dimensional approach by Radjai, Wolf, Jean, and Moreau ther increases and finally reaches a critical value. In addition, the
(1998)) and a three-dimensional method by Thornton and Anthony strong subnetworks generally follow the trend of the stress-strain
(1988) indicated that if the contact normal force is greater than the behaviour as observed in Fig. 6, which supports the initiation of the
average contact normal force, it can be treated as a strong subnet- strong force chain buckling that can be characterized by the appear-
work, otherwise, it is a weak subnetwork. ance of peak structural anisotropy (Huang, Hanley, O’Sullivan, &

159
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

Fig. 12. Relationship between peak strength and dilatancy at peak: (a) p − Dp , and (b) p − Dp .

 
Fig. 13. Evolution of Z for Series II and Series V: (a) εa − Z, (b) εa − Z, (c) p − Z at critical state, (d) p − Z at critical state.

Kwok, 2014; Tordesillas & Muthuswamy, 2009). Figs. 14(c) and state is achieved, which is independent of the samples’ initial
(d) illustrate the relationships between the overall and the strong packings, under the axial compression (AC) and axial extension

deviator fabrics against p , respectively for Series II and Series V at (AE) stress paths respectively.
the critical states. As can be seen, the AC stress path test is more • For a given confining pressure, the value of the critical void ratio
anisotropic than the AE stress path in general. Moreover, the struc- under the AC stress path is smaller than that under the AE stress

tural anisotropy at the critical state decreases with increasing po path. The critical void ratio decreases as the confining stress
(confining pressure), which is significantly modified in the strong increases for both the AC and AE stress paths.
deviator fabric d, strong in comparison to the overall deviator fab- • The critical state lines (CSLs) are found to be path-dependent and
ric d . the critical strength (stress ratio) for the AC stress path tests is
higher than that for the AE stress path tests.
Conclusions • The relationship between peak friction angles and peak dilatancy
is established for the AC and AE stress paths respectively, but the
This study presents the DEM simulation results of strength and relationship is found to be stress-path dependent.
critical state behaviour of granular materials under AC and AE stress • The coordination number at critical state is found to increase with
path tests. Twenty four numerical samples were generated consid- increasing confining pressures for each stress path, reflecting a
ering the different states of assemblies based on densification index decrease in the critical void ratios. It is found that the critical
(ID ). The results were analyzed at the macroscopic and microscopic coordination number for the AC stress path is higher than that for
levels, examining the uniqueness of CSLs and correlations among the AE stress path under a given confining pressure, signifying a
strength, dilation and state parameters. Based on the analyses of higher critical stress ratio for AC than for AE in relation to a given
results and discussions, the following key conclusions are drawn: sample with a given confining pressure.
• A unique linear relationship between the critical coordination
• All the samples reach critical state after an axial strain of about numbers and critical void ratios is identified under the AC and
50%. For a given confining pressure, a unique void ratio at critical

160
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

 
Fig. 14. Evolution of fabrics for Series II and Series V: (a) εa − d, (b) εa − d, strong, (c) p − d at critical state, (d) p − d, strong at critical state.

AE stress paths respectively, but such a relationship depends on Been, K., & Jefferies, M. G. (1985). A state parameter for sands. Geotechnique, 35(2),
the stress paths. 99–112.
Been, K., Jefferies, M. G., & Hachey, J. (1991). The critical state of sands.
• The overall structural anisotropy in terms of deviator fabric at the
Geotechnique, 41(3), 365–381.
critical state is found to be unique for the AC and AE stress paths Castro, G. (1969). Liquefaction of sand. PhD thesis, Division of Engineering and Applied
respectively, but the critical deviator fabric is larger for the AC Physics, Harvard University.
Cho, G., Dodds, J., & Santamarina, J. (2006). Particle shape effects on packing
stress paths. density, stiffness, and strength: natural and crushed sands. Journal of
Geotechnical and Geoenvironmental Engineering, 132(5), 591–602.
Christoffersen, J., Mehrabadi, M., & Nemat-Nasser, S. (1981). A micromechanical
Declaration of interests description of granular material behavior. Journal of Applied Mechanics, 48(2),
339–344.
Chu, J. (1995). An experimental examination of the critical state and other similar
The authors declare that they have no known competing finan-
concepts for granular soils. Canadian Geotechnical Journal, 32(6), 1065–1075.
cial interests or personal relationships that could have appeared to Cleary, P. W. (2010). DEM prediction of industrial and geophysical particle flows.
influence the work reported in this paper. Particuology, 8, 106–118.
Cundall, P. A., & Strack, O. D. (1979). A discrete numerical model for granular
assemblies. Geotechnique, 29(1), 47–65.
da Cruz, F., Emam, S., Prochnow, M., Roux, J. N., & Chevoir, F. (2005). Rheophysics of
Acknowledgements
dense granular materials: discrete simulation of plane shear flows. Physical
Review: E, 72(2), Article 021309.
The financial supports provided by Xi’an Jiaotong – Liverpool Fu, P., & Dafalias, Y. F. (2015). Relationship between void-and contact
University (Grant Nos. RDF 14-02-44, RDF 15-01-38, RDF 18-01- normal-based fabric tensors for 2D idealized granular materials. International
Journal of Solids and Structures, 63, 68–81.
23 and PGRS1906002) and the Key Program Special Fund at XJTLU Goldenberg, C., & Goldhirsch, I. (2005). Friction enhances elasticity in granular
(Grant No. KSF-E-19) and Natural Science Foundation of Jiangsu solids. Nature, 435(7039), 188–191.
Province (Grant No. BK20160393) are gratefully acknowledged. In Gong, G., & Zha, X. (2013). DEM simulation of undrained behaviour with
preshearing history for saturated granular media. Modelling and Simulation in
addition, Prof. Kristian Krabbenhoft from the University of Liver- Materials Science and Engineering, 21(2), Article 025001.
pool, UK and Prof. Zhongxuan Yang from Zhejiang University, China, Gong, G., Zha, X., & Wei, J. (2012). Comparison of granular material behaviour
provided valuable comments on the relevant work, which is appre- under drained triaxial and plane strain conditions using 3D DEM simulations.
Acta Mechanica Solida Sinica, 25(2), 186–195.
ciated. Gong, J., & Liu, J. (2017). Effect of aspect ratio on triaxial compression of
multi-sphere ellipsoid assemblies simulated using a discrete element method.
Particuology, 32, 49–62.
References Green, G. E., & Bishop, A. W. (1969). A note on the drained strength of sand under
generalized strain conditions. Geotechnique, 19, 144–149.
Abedi, S., & Mirghasemi, A. A. (2011). Particle shape consideration in numerical Gu, X. Q., Huang, M. S., & Qian, J. G. (2014). DEM investigation on the evolution of
simulation of assemblies of irregularly shaped particles. Particuology, 9(4), microstructure in granular soils under shearing. Granular Matter, 16(1),
387–397. 91–106.
Alizadeh Behjani, M., Asachi, M., Ghadiri, M., Bayly, A., & Hassanpour, A. (2018). A Guo, N., & Zhao, J. (2013). The signature of shear-induced anisotropy in granular
methodology for calibration of DEM input parameters in simulation of media. Computers and Geotechnics., 47, 1–15.
segregation of powder mixtures, a special focus on adhesion. Powder Gutierrez, M., Wang, J., & Yoshimine, M. (2009). Modeling of the simple shear
Technology, 339, 789–800. deformation of sand: effects of principal stress rotation. Acta Geotechnica, 4,
Alizadeh Behjani, M., Hassanpour, A., Pasha, M., Ghadiri, M., & Bayly, A. (2017). The 193–201.
effect of particle shape on predicted segregation in binary powder mixtures. Härtl, J., & Ooi, J. Y. (2011). Numerical investigation of particle shape and particle
Powder Technology, 319, 313–322. friction on limiting bulk friction in direct shear tests and comparison with
Asachi, M., Alizadeh Behjani, M., Nourafkan, E., & Hassanpour, A. (2020). Tailoring experiments. Powder Technology, 212, 231–239.
particle shape for enhancing the homogeneity of powder mixtures: Head, K. H. (1986). Manual of soil laboratory testing. London: Pentech Press.
Experimental study and DEM modelling. Particuology (In press). Hertz, H. (1882). Über die berührung fester elastischer Körper. Journal für die Reine
[Link] und Angewandte Mathematik, 29, 156–171.

161
S.P.K. Kodicherla et al. Particuology 56 (2021) 152–162

Huang, X., Hanley, K. J., O’Sullivan, C., & Kwok, C. Y. (2014). Exploring the influence Roscoe, K. H., Schofield, A. N., & Wroth, C. P. (1958). On the yielding of soils.
of interparticle friction on critical state behaviour using DEM. International Geotechnique, 8, 22–53.
Journal for Numerical and Analytical Methods in Geomechanics, 38, 1276–1297. Rowe, P. W. (1962). The stress–dilatancy relation for static equilibrium of an
Huang, X., Hanley, K. J., O’Sullivan, C., Kwok, C. Y., & Wadee, M. A. (2014). DEM assembly of particles in contact. Proceedings of Royal Society A: Mathematical,
analysis of the influence of the intermediate stress ratio on the critical-state Physical and Engineering Sciences, 269(1339), 500–527.
behaviour of granular materials. Granular matter, 16, 641–655. Salvatore, E., Modoni, G., Ando, E., Albano, M., & Viggiani, G. (2017). Determination
Huang, X., O’Sullivan, C., Hanley, K. J., & Kwok, C. Y. (2014). Discrete-element of the critical state of granular materials with triaxial test. Soils and
method analysis of the state parameter. Geotechnique, 64(12), 954–965. Foundations, 57(5), 733–744.
Itasca. (2003). PFC3D: Theory and background. Itasca Consulting Group. Satake, M. (1982). Fabric tensor in granular materials. In In Proceedings of the IUTAM
Itasca Consulting Group. (2019). PFC3D 5.0. 2019. Itasca Consulting Group, Symposium on Deformations and Failure of Granular Materials 1982 (Vermeer PA
Minneapolis. and Luger HJ (eds)). pp. 63–68. the Netherlands: Balkema, Rotterdam.
Jefferies, M. G., & Been, K. (2006). Soil liquefaction: a critical state approach. London, Schofield, A., & Wroth, P. (1968). Critical state soil mechanics. London: McGraw-Hill.
UK: Taylor and Francis. Sitharam, T. G., & Vinod, J. S. (2009). Critical state behaviour of granular materials
Jefferies, M. G., & Shuttle, D. A. (2002). Dilatancy in general Cambridge-type from isotropic and rebounded paths: DEM simulations. Granular matter, 11,
models. Geotechnique, 52(9), 625–638. 33–42.
Jefferies, M. G. (1993). Nor-Sand: a simple critical state model for sand. Springman, S. M., Yamamoto, Y., Buchli, T., Hertrich, M., Maurer, H., Merz, K.,
Geotechnique, 43(1), 91–103. Gartner-Roer, I., & Seward, L. (2013). In C. Margottini, P. Canuti, & K. Sassa
Jiang, M. J., Li, L., & Yang, Q. (2013). Experimental investigation on deformation (Eds.), Rock glacier degradation and instabilities in the European Alps: a
behavior of TJ-1 lunar soil simulant subjected to principal stress rotation. characterisation and monitoring experiment in the Turtmanntal, CH. In Landslide
Advanced Space Research, 52, 136–146. Science and Practice (pp. 5–13). Berlin, Heidelberg: Springer.
Ko, H. Y., & Scott, R. F. (1967). A new soil testing apparatus. Géotechnique, 17, 40–57. Stahl, M., & Konietzky, H. (2011). Discrete element simulation of ballast and gravel
Kodicherla, S. P. K., Gong, G., Yang, Z. X., Krabbenhoft, K., Fan, L., Moy, C. K. S., et al. under special consideration of grain-shape, grain-size and relative density.
(2019). The influence of particle elongations on direct shear behaviour of Granular Matter, 13(4), 417–428.
granular materials using DEM. Granular Matter, 21, 86. Taghavi, R. (2000). Automatic Block Decomposition Using Fuzzy Logic Analysis.
Kodicherla, S. P. K., Gong, G., Fan, L., & Moy, C. K. S. (2019). Effects of particle Proceedings of the 9th International Meshing Roundtable. pp. 187–192. New
morphology on the macroscopic behaviour of ellipsoids: A discrete element Orleans: Sandia National Laboratories.
investigation. Proceedings of the 2nd International Conference on Sustainable Taylor, D. W. (1948). Fundamentals of soil mechanics. New York, NY, USA: Wiley.
Buildings and Structures (ICSBS 2019), Papadikis et al. (Eds), Suzhou, China. pp. Thornton, C., & Anthony, S. J. (1988). Quasi-static deformation of particulate media.
45–50. Philosophical Transactions of the Royal Society A, Mathematical, Physical and
Kodicherla, S. P. K., Gong, G., Fan, L., Wilkinson, S., & Moy, C. K. S. (2020). DEM Engineering Sciences, 356(1747), 2763–2782.
investigations of the effects of particle morphology on granular material Thornton, C. (2000). Numerical simulations of deviator shear deforma- tions of
behaviour using a multi-sphere approach. Journal of Rock mechanics and granular media. Geotechnique, 50(1), 43–53.
Geotechnical Engineering. (In press). Tordesillas, A., & Muthuswamy, M. (2009). On the modeling of confined buckling of
[Link] force chains. Journal of the Mechanics and Physics of Solids, 57, 706–727.
Kulhawy, F. H., & Mayne, P. W. (1990). Manual on estimating soil properties in Vaid, Y. P., & Sasitharan, S. (1991). The strength and dilatancy of sand. Canadian
foundation design. Palo Alto, CA: EPRI. Geotechnical Journal, 29(10), 522–526.
Lade, P. V., & Duncan, J. M. (1975). Elastoplastic stress-strain theory for Verdugo, R., & Ishihara, K. (1996). The steady state of sandy soils. Soils and
cohesionless soil. Journal of Geotechnical Engineering Division, 101(10), Foundations, 36(2), 81–91.
1037–1053. Wang, J., & Gutierrez, M. (2010). Discrete element simulations of direct shear
Lade, P. V. (1972). The stress-strain and strength characteristics of cohesionless soils. specimen scale effects. Geotechnique, 60(5), 395–409.
PhD thesis, University of California at Berkeley. Wellmann, C., Lillie, C., & Wriggers, P. (2008). Comparison of the macroscopic
Li, X., & Li, X. S. (2009). Micro-macro quantification of the internal structure of behavior of granular materials modeled by different constitutive equations on
granular materials. Journal of Engineering Mechanics, 135(7), 641–656. the microscale. Finite Element Analysis and Design, 44, 259–271.
Li, X., & Yu, H. S. (2010). Numerical investigation of granular material behaviour Xie, Y. H., Yang, Z. X., Barreto, D., & Jiang, M. D. (2017). The influence of particle
under rotational shear. Geotechnique, 60(5), 381–394. geometry and the intermediate stress ratio on the shear behavior of granular
Li, X. (2006). Micro-scale investigation on the quasi-static behavior of granular materials. Granular matter, 19, 35.
material. Doctorial Dissertation, Hong Kong University of Science and Technology. Yan, W. M. (2009). Fabric evolution in a numerical direct shear test. Computers and
Li, X. S., & Wang, Y. (1998). Linear representation of steady-state line for sand. Geotechnics, 36(4), 597–603.
Journal of Geotechnical and Geoenvironmental Engineering, 124(12), 1215–1217. Yang, Y., Wang, J. F., & Cheng, Y. M. (2016). Quantified evaluation of particle shape
Maeda, K., Fukuma, M., & Nukudani, E. (2009). Macro and micro critical states of effects from micro-to-macroscales for non-convex grains. Particuology, 25,
granular materials with different grain shapes. AIP Conference Proceedings, 23–35.
1145(1), 829–832. Yang, Z. X., & Wu, Y. (2016). Critical state for anisotropic granular materials: A
Mindlin, R. D., & Deresiewicz, H. (1953). Elastic spheres in contact under varying discrete element perspective. International Journal of Geomechanics, Article
oblique force. Journal of Applied Mechanics, 20, 327–344. 04016054.
Mooney, M. A., Finno, R. J., & Viggiani, M. G. (1998). A unique critical state for sand? Yang, Z. X., Li, X. S., & Yang, J. (2008). Quantifying and modelling fabric anisotropy
Journal of Geotechnical and Geoenvironmental Engineering, 124(11), 1100–1108. of granular soils. Geotechnique, 58(4), 237–248.
Ng, T. T. (2004). Triaxial simulations using DEM with hydrostatic boundaries. Yang, Z. X., Yang, J., & Wang, L. Z. (2013). Micro-scale modelling of anisotropy
Journal of Engineering Mechanics, 130(10), 1188–1194. effects on undrained behaviour of granular soils. Granular matter, 15(5),
Ng, T. T. (2009). Discrete element method simulations of the critical state of a 557–572.
granular material. International Journal of Geomechanics, 9(5), 209–216. Yimsiri, S., & Soga, K. (2011). Effects of soil fabric on behaviors of granular soils:
Norwegian Geotechnical Institute (NGI). (1982). Results of triaxial tests on Hokksund Microscopic modeling. Computers and Geotechnics, 38, 861–874.
sand. Internal Report 52108-52112, Jan. Yoshimine, M., Ishihara, K., & Vargas, W. (1998). Effects of principal stress direction
Oda, M. (1982). Fabric tensor for discontinuous geological materials. Soils and and intermediate principal stress on undrained shear behavior of sand. Soils
Foundations, 22(4), 96–108. and Foundation, 38(3), 179–188.
Pournin, L., Weber, M., Tsukahara, M., Ferrez, J. A., Ramaioli, M., & Liebling, T. M. Zhang, J., Lo, S. C. R., Rahman, M. M., & Yan, J. (2018). Characterizing monotonic
(2005). Three-dimensional distinct element simulation of spherocylinder behavior of pond ash within critical state approach. Journal of Geotechnical and
crystallization. Granular Matter, 7, 119–126. Geoenvironmental Engineering, 144(1), Article 04017100.
Radjai, F., Wolf, D. E., Jean, M., & Moreau, J. J. (1998). Bimodal character of stress Zhao, J., & Guo, N. (2013). Unique critical state characteristics in granular media
transmission in granular packings. Physical Review Letters, 80, 61–64. considering fabric anisotropy. Geotechnique, 63(8), 695–704.
Rodriguez, N. M., & Lade, P. V. (2013). Effects of principal stress directions and Zhao, S., & Zhao, J. (2019). A poly-superellipsoid-based approach on particle
mean normal stress on failure criterion for cross-anisotropic sand. Journal of morphology for DEM modeling of granular media. International Journal for
Engineering Mechanics, 139(11), 1592–1601. Numerical and Analytical Methods in Geomechanics, 43(13), 2147–2169.
Roscoe, K., & Burland, J. B. (1968). On the generalized stress-strain behavior of wet Zhao, S., Zhou, X., & Liu, W. (2015). Discrete element simulations of direct shear
clay. In J. Heyman, & F. A. Leckie (Eds.), Engineering Plasticity (pp. 535–609). tests with particle angularity effect. Granular Matter, 17(6), 793–806.
Cambridge, U.K: Cambridge University Press.

162