Granular Matter (2024) 26:57
https://doi.org/10.1007/s10035-024-01421-7
ORIGINAL REPORT
Cyclic-loading effects in sand: a micromechanical study considering
particle breakage
Jacinto Ulloa
1
·
Ziran Zhou
1
·
John Harmon
1
·
José E. Andrade
1
Received: 21 August 2023 / Accepted: 9 March 2024 / Published online: 19 April 2024
© The Author(s), under exclusive licence to Springer-Verlag GmbH Germany, part of Springer Nature 2024
Abstract
This paper investigates the response of Ottawa sand to cyclic loading using virtual oedometer tests and the level-set discrete
element method. We study both the macroscopic and the micromechanical behavior, shedding light on the grain-scale processes
behind the cyclic response observed in crushable sand, namely stress relaxation under strain control and ratcheting under stress
control. Tests without particle breakage first show that asymmetrical frictional sliding during loading-unloading induces these
cyclic-loading effects. Then, tests considering particle breakage reveal more pronounced stress relaxation and ratcheting,
which decrease in rate over cycles, accompanied by increased frictional sliding and reduced particle contact forces. It is found
that the broken fragments unload the most and promote an enhanced cushioning effect. These micromechanical processes
contribute to a decrease in breakage potential as the cycles progress, implying that cyclically loaded materials may become
more resistant to breakage when compared to the same material loaded monotonically at the same strain level. These new
insights highlight the main contributions of the present work, factoring in real particle shapes from 3D X-ray tomography
and notably contributing to the existing literature on the topic, where most studies rely on idealized particle shapes and rarely
consider crushable grains.
Keywords
Ratcheting
·
Cyclic loading
·
Breakage
·
Crushable sand
·
Grain-size distribution
·
LS-DEM
1 Introduction
Cyclic loading conditions are ubiquitous in both natural and
engineered structural systems, emerging due to the repetitive
action of traffic, heavy machinery, wind, ocean currents, seis-
mic events, and others. Such loading conditions induce pecu-
liar effects in the behavior of materials that strongly modulate
their long-term performance. In the context of geomechan-
ics, the cyclic response of granular soils is characterized
by non-vanishing residual deformations after each loading
cycle, resulting in a gradual plastic strain accumulation. This
phenomenon is called
ratcheting
(in analogy to mechani-
cal ratchet devices) and is also observed in other materials
such as metals [
1
]. Under constant-amplitude cycles of
applied strain, the ratcheting effect naturally translates to
stress relaxation
, where the cyclic stress response gradu-
ally decays as the load cycles progress. These mechanisms
B
José E. Andrade
jandrade@caltech.edu
1
Division of Engineering and Applied Science, California
Institute of Technology, Pasadena, CA 91125, USA
have significant implications for engineering structures that
rely on soil foundations for stability. For instance, prevent-
ing excessive ratcheting is essential to guarantee the integrity
of offshore wind turbines with monopile foundations [
2
]or
structures prone to damage due to differential foundation set-
tlements [
3
].
A particular feature of granular materials is that ratcheting
is observed even when the load amplitude is very small, pre-
cluding the stabilization of the response and the emergence
of an elastic stage [
4
]. This behavior is attributed to com-
plex grain-scale interactions involving particle sliding and
crushing, both of which entail irreversible processes. The
interplay between these processes and the cyclic response of
granular soils is still not fully understood, perhaps due, on
the one hand, to the limited access of experimental studies to
detailedmicromechanicsand,ontheother,tothelimitedabil-
ity of conventional numerical models to describe complex
systems realistically. This knowledge gap possibly explains
why most constitutive models that consider cyclic-loading
effects [
3
,
5
,
6
] do not account for breakage, while continuum
breakage mechanics models [
7
–
9
] do not account for cyclic-
loading effects in the sense that load cycles remain within
123
57
Page 2 of 16
J. Ulloa et al.
a presumed elastic domain. The development of enhanced
versions of these models will likely rely on further insight
from grain-scale processes.
Experimentalworkonthecyclicloadingofgranularmedia
has provided clear evidence of cyclic strain accumulation.
Wichtmann et al. [
10
] studied the effect of different loading
conditions, concluding that ratcheting is relatively insensi-
tive to the loading frequency and strongly influenced by
the number of stress cycles, the average stress, and the
stress amplitude. The physical material properties also mod-
ulate the cyclic response. For instance, an increase in strain
accumulation with increasing initial void ratio has been
reported [
10
–
12
]. Park and Santamarina [
11
] and Cha et
al. [
12
] further highlight that sands asymptotically approach
a terminal void ratio whose value is a function of the
initial density, implying that the material preserves the mem-
ory of its initial fabric. Moreover, a strong influence of
both grain size distribution (GSD) [
10
,
13
] and particle
shape [
14
] has been observed, showing an increase in strain
accumulation with decreasing particle size and increasing
uniformity coefficient. Similar conclusions have been drawn
for ballast samples [
15
]. These observations suggest that
the cyclic response of materials prone to particle breakage
strongly depends on the degree of comminution. In fact,
the evolution of GSD towards finer particles often leads to
looser and more compressible samples, thus promoting strain
accumulation [
16
]. However, after a sufficient number of
cycles, granular materials often reach a stable grading, after
which the rates of breakage and ratcheting decrease substan-
tially [
17
].
Theseexperimentalstudiesprovideinsightintothemacro-
scopic response of granular soils under cyclic loading,
including strain accumulation and related changes in the
mechanical and physical properties of the material. How-
ever, the micromechanical explanation of these phenomena
is merely hypothesized, rendering numerical studies neces-
sary to provide additional insight into the underlying physics.
Alonso-Marroquín and Herrmann [
4
] presented one of the
first studies in this context using the discrete element method
(DEM) for 2D polygons. The authors presented evidence
for the absence of truly elastic behavior in granular pack-
ings under cyclic loading, identifying (i) a short-time regime
with a fast cyclic strain accumulation, particularly during
the first cycles, and (ii) long-time ratcheting regimes char-
acterized by a constant cyclic strain accumulation and often
separated by short time regimes. The ratcheting process was
attributed to an almost periodic ratchet-like stick-slip behav-
ior of sliding contacts. Hu et al. [
18
] further studied the effect
of the packing structure and load amplitude in a simplified
2D scenario and highlighted the role of fabric anisotropy
in promoting strain accumulation. Nguyen et al. [
19
] later
studied the effect of different cyclic loading conditions and
packing densities in a 3D assembly of spheres, observing that
ratcheting is accompanied by a significant frictional energy
dissipation at the sliding contacts, primarily modulated by the
average stress, stress amplitude, and sample density. Sassel
et al. [
20
] also reported a strong influence of load ampli-
tude in an assembly of spheres and observed a steady-state
contact network with persistent but constant cyclic changes.
These numerical studies generally show that cyclic strain
accumulation is accompanied by frictional dissipation while
the average strain energy and the fabric present relatively
small changes from one cycle to another. The root cause of
this behavior seems to be the asymmetrical contact sliding
that takes place during the unloading and reloading stages of
a load cycle, with a higher fraction of sliding contacts at the
peak load.
None of the numerical studies mentioned so far consid-
ered particle breakage, often noting that the applied stress
levels are below the ultimate strength of individual grains.
However, breakage may occur in samples under compression
at stresses below the level required to break individual parti-
cles due to the heterogeneous force distribution in the contact
fabric. The evolving contact fabric also implies that unbro-
ken grains within a load cycle may be subject to larger forces
and eventually break in subsequent cycles, even if the applied
macroscopic stress does not increase [
17
]. Numerical studies
considering breakage are thus required to assess the response
of granular media to cyclic loading. A reduced body of lit-
erature exists on this topic. Some studies have considered
idealized shapes constructed through particle replacement or
particle clumping. For instance, Lobo-Guerrero and Vallejo
[
21
] and Indraratna et al. [
22
] presented 2D studies high-
lighting the strong correlation between particle breakage and
strain accumulation in the cyclic response of railway ballast.
Dahal and Mishra [
23
] confirmed these observations in 3D
assemblies of clumped spheres and further highlighted the
crucial role of particle shape, in particular for breakable par-
ticles, noting substantial differences in the cyclic response
between complex-shaped clumps and ellipsoids. These stud-
ies have mainly focused on the macroscopic response rather
than the micromechanics and are limited by the idealization
of particle shapes. Nie et al. [
24
] presented a very recent
exception using bonded polyhedra to represent realistic par-
ticle morphologies, where the strong correlation between
breakage and strain accumulation was confirmed. In general,
all these works point towards an initial rapid strain accumu-
lation with significant breakage and a slower accumulation
once breakage stabilizes, mainly due to the weakening of
contact forces and the densification of the sample.
The present study is motivated by the gap between
numerics and experiments on granular media under cyclic
loading: while physical experiments provide macroscopic
observations without access to grain-scale kinetics, existing
numerical studies have rarely considered breakage and are
mostly limited to idealized particle shapes. Thus, we conduct
123
Cyclic-loading effects in sand: a micromechanical...
Page 3 of 16
57
a series of numerical cyclic-loading tests using the level set
discrete element method (LS-DEM) [
25
] on a virtual avatar
of a real sample of crushable Ottawa sand. The monotonic
response of the numerical model was verified experimen-
tally in a recent study [
26
], showing a remarkable agreement
with experimental results. In the present work, we show that
some of the observations made by previous works on ide-
alized samples are also present in Ottawa sand under cyclic
loading while further elucidating the underlying microme-
chanics of cyclic particle breakage and its link to ratcheting
under stress control and stress relaxation under strain control.
We first address the purely cyclic-loading effects in a sam-
ple without particle breakage, discussing the role of friction.
Then, we present simulations considering particle breakage
to assess the effect of crushing on the cyclic response. More-
over, we address the fundamentally different question of how
cyclic loading conditions affect the breakage process when
compared to samples loaded monotonically.
2 Virtual testing
2.1 Cyclic oedometer test
Figure
1
a shows the virtual sample of Ottawa sand sub-
jected to one-dimensional compression. The experimental
setup and the preparation of the virtual model were presented
in a recent study [
26
]; a brief summary is provided here for
the reader’s convenience. The physical experiment consists
of an oedometer test, where a sample of oven-dried sand was
first pluviated into a cylinder with a diameter of 4.1 mm and
height of 5 mm. The specimen had an initial void ratio of
0
.
6 and a median diameter of 600
μ
m. Figure
2
shows the
initial GSD curve
F
0
. The physical specimen was subjected
to monotonic strain-control loading at a displacement rate of
0.05 mm/min, taking successive X-ray scans of the compres-
sive breakage process.
In LS-DEM simulations conducted by Harmon et al. [
26
],
a one-to-one virtual avatar of the experimental sample was
generated (Fig.
1
a), following the image processing method-
ology of Vlahini ́
cetal.[
28
]. The initial virtual sample
consists of 404 grains, with a level set grid resolution of 8.6
nm/voxel and surface nodes discretized at 510 points/mm
2
.
The response of this virtual sample was thoroughly studied
in Harmon et al. [
26
] under monotonic loading, showing a
remarkable agreement with the experimental data and further
elucidating the statistics of contact forces. Such a grain-scale
study was possible thanks to the use of breakage-enhanced
LS-DEM, outlined in Section
2.2
below.
In the present study, we employ this computational frame-
work to study the behavior of the oedometer test under cyclic
loading, considering either stress-control or strain-control
simulations (Fig.
1
b and c, respectively). For this purpose, we
first subjected the initial virtual sample to monotonic loading
up to 23% strain. At this point, the initial grains had crushed
into 13720 particles. Figure
2
illustrates the resulting GSD
at this state, denoted as
F
u
. For simplicity, we consider
F
u
as the ultimate GSD. While this curve may not represent the
actual ultimate GSD of the material, generally assumed to be
fractal,
F
u
will serve as a reference state in the subsequent
simulations to track the evolution of the GSD.
2.2 Numerical modeling framework
LS-DEM [
25
] is a variant of discrete element modeling that
allows for arbitrary particle shapes. In this framework, a level
set function
φ
p
(
x
)
is used to describe the morphology of a
particle
p
, measuring the signed distance between a point
x
and the particle’s surface
{
x
:
φ
p
(
x
)
=
0
}
. This rep-
resentation of particle shape allows us to make a seamless
transition from experimental images to high-fidelity virtual
Fig. 1
Graphical representation of the oedometer test, showing (a) the geometry and boundary conditions for (b) stress-controlled and (c) strain-
controlled constant-amplitude cyclic loading
123
57
Page 4 of 16
J. Ulloa et al.
Fig. 2
Initial GSD
F
0
of the specimen in Fig.
1
and assumed ultimate
GSD
F
u
corresponding to an oedometer test up to 23% strain. The
ultimate GSD is employed to generate breakage values [
27
] between 0
and 1
avatars, preserving the shape, position, and orientation of
every particle. Moreover, the contact detection procedure is
very simple, requiring only the evaluation of a particle’s level
set function at the surface points of a potentially contact-
ing neighbor. Rigid particle kinematics are then resolved in
conjunction with simple contact kinetics, for instance, con-
sidering a linear Hookean law capped by Coulomb friction.
The total force
f
c
at a contact point
c
is then decomposed
into normal and shear components, i.e.,
f
c
=
f
c
n
+
f
c
s
, where
f
c
n
=
k
n
δ
c
n
n
c
,
(1)
{
f
c
s
←
Rf
c
s
−
k
s
s
c
,
f
c
s
←
f
c
s
f
c
s
min
{
f
c
s
,μ
f
c
n
}
.
(2)
Here,
n
c
represents the normal at the contact point,
δ
c
n
is
the normal overlap, and
s
c
is the incremental tangential
displacement. The values
k
n
,
k
s
, and
μ
represent the normal
stiffness,shearstiffness,andfrictioncoefficient,respectively.
The rotation matrix from the previous normal direction to the
current is denoted by
R
.
An extension of LS-DEM to breakage was proposed
in Harmon et al. [
29
], where a parent particle defined by
φ
p
(
x
)
is allowed to break into two fragments
φ
p
1
(
x
)
=
max
{
φ
p
(
x
), φ
split
(
x
)
}
and
φ
p
2
(
x
)
=
max
{
φ
p
(
x
),
−
φ
split
(
x
)
}
, where
φ
split
(
x
)
is a splitting surface (Fig.
3
). While in
principle any breakage criterion may be employed to decide
when and how a particle breaks, the reference study of Har-
mon et al. [
26
] considered a stress-based criterion, for which
the average stress in a particle
p
is first computed as
σ
p
=
1
V
p
N
cont
p
∑
α
=
1
sym
(
f
α
p
⊗
l
α
p
),
(3)
Fig. 3
Graphical representation of the breakage process in LS-DEM
where
V
p
is the particle volume,
N
cont
p
is the number of con-
tacts with neighboring grains or boundaries,
f
α
p
is a contact
force, and
l
α
p
is the corresponding branch vector joining the
contact point to the centroid of the grain. Breakage due to ten-
sile fracture is set to occur when the major principal stress
σ
t
p
=
max
{
σ
p
1
,σ
p
2
,σ
p
3
}
reaches a threshold, i.e.,
σ
t
p
≤
σ
cr
p
.
(4)
The critical value
σ
cr
p
is defined according to Weibull’s
theory for brittle fracture due to internal flaws, taking into
account the dependence of strength on the particle diameter
d
p
. Consequently, we consider the relation
σ
cr
p
=
σ
cr
0
[(
d
0
d
p
)
3
ln
(
1
−
P
f
)
]
1
/
m
,
(5)
where
m
is Weibull’s modulus and
P
f
∼
U
(
0
,
1
)
.Therefer-
ence stress
σ
cr
0
corresponds to a reference grain with diameter
d
0
for which
P
f
=
1
−
e
−
1
≈
0
.
63, i.e., a grain with an
approximately 63% failure probability or, equivalently, an
approximately 37% survival probability. Equation (
5
) deter-
mines the strength of the initial particles in a sample, while
the strength of a fragment
q
coming from a parent grain
p
is
computed as
σ
cr
q
=
σ
cr
p
(
d
p
d
q
)
m
/
3
.
(6)
While the method admits arbitrary fracture paths, using
planar surfaces determined by joining the three highest con-
tact forces has been shown to yield accurate results [
26
].
Alternatively, the fracture plane can be determined from the
grain centroid and the direction of the maximum principal
stress.
123