of 12
Received: 3 July 2023
Revised: 3 October 2023
Accepted: 2 November 2023
DOI: 10.1002/nme.7400
RESEARCH ARTICLE
A discretization-convergent
level-set-discrete-element-method using a continuum-based
contact formulation
Shai Feldfogel
1
Konstantinos Karapiperis
2
Jose Andrade
3
David S. Kammer
1
1
Institute for Building Materials, ETH,
Zurich, Switzerland
2
Department of Mechanical and Process
Engineering, ETH, Zurich, Switzerland
3
Department of Mechanical and Civil
Engineering, Caltech, Pasadena,
California, USA
Correspondence
David S. Kammer, Institute for Building
Materials, ETH, Zurich, Switzerland.
Email:
dkammer@ethz.ch
Funding information
Staatssekretariat für Bildung, Forschung
und Innovation; ESKAS, Grant/Award
Number: 2021.0165
Abstract
The level-set-discrete-element-method (LS-DEM) was developed to overcome
the shape limitation of traditional discrete element method. LS-DEM’s shape
generality relies on a node-based surface discretization of grain boundary, and it
has been used to shed new light of a variety of granular mechanics applications
with realistically shaped grains and structural assemblies made of unbonded
building blocks. Due to the node-based discretization of grain boundary, the
original LS-DEM is discretization-sensitive and it suffers from divergence of the
response with discretization refinement, particularly for highly compressible
problems. Previous studies have identified and addressed this issue in differ-
ent ways, each with its own advantages and shortcomings. Here, we propose
a methodologically-rigorous and computationally-efficient adapted formulation
which solves LS-DEM’s discretization-sensitivity issue. It adopts the classical
contact description of continuum mechanics, wherein the contact interactions
are traction-based. We demonstrate the convergence of the adapted LS-DEM
in several highly compressible cases studies, show that it is key to correctly
capturing the mechanical response, and compare it to alternative formulations.
KEYWORDS
contact, convergence, DEM, granular material, level set, level-set-DEM, mesh-sensitivity,
topologically interlocked structures
1
INTRODUCTION
Modeling granular assemblies with realistically-shaped non-spherical grains has been a long-standing challenge in
the discrete element method (DEM) community.
1-5
The level-set-discrete-element-method (LS-DEM) was developed to
address this very challenge.
6,7
LS-DEM’s ability to represent arbitrarily shaped grains and resolve contact between them
relies on a node-based discretization of the grain boundary.
The shape versatility and the other mathematical properties of the level-set geometrical representation used in
LS-DEM(e.g.,theabilitytoneatlysplitgrains)havebeenusedtoshednewlightonadiverserangeofproblemsingranular
mechanics. This includes triaxial compression test of Martian-like sand,
7
constitutive elasto-plastic response and shape
effects upon loading and unloading of Hostun sand,
8
the effects of reduced gravity on the strength of sand,
9
the effects of
This is an open access article under the terms of the
Creative Commons Attribution
License, which permits use, distribution and reproduction in any medium, provided the
original work is properly cited.
© 2023 The Authors.
International Journal for Numerical Methods in Engineering
published by John Wiley & Sons Ltd.
Int J Numer Methods Eng
. 2024;125:e7400.
wileyonlinelibrary.com/journal/nme
1of12
https://doi.org/10.1002/nme.7400
2of12
FELDFOGEL et al.
grain breakage on crushable sand and rock assemblies,
10,11
and the formation and evolution of stress transmission paths
in entangled granular assemblies.
12
The range of application of LS-DEM has recently been extended to structures made of un-bonded building blocks,
including multi-block towers
13-15
and topologically interlocked structures (TIS).
16-18
Using calibration of the normal con-
tact stiffness
k
n
, it has been shown to correctly capture and predict the structural response and the failure mechanism
under dynamic
15
and static
16
loads. Based on a closed-form expression relating
k
n
and the elastic modulus
E
, it has been
used to obtain new insights on the effects of
E
and the friction coefficients on slab-like TIS,
17
and on their deflection
capacity.
18
One of the main shortcomings of the original LS-DEM formulation, first identified and addressed in Reference
19
and later acknowledged in References
20-23
, is that the response is sensitive to the surface discretization of grain bound-
ary. Specifically, the grains become stiffer as the number of surface nodes increases,
19,22
and, as a result, the mechanical
response diverges with discretization refinement. LS-DEM’s discretization sensitivity is manifested mostly in highly com-
pressible cases, that is when the forces or pressures are large enough to induce particle deformation beyond the rigid-body
limit. In contrast, when the response is governed by rearrangements and rigid-body kinematics (e.g., gravity-driven parti-
cle flow), the discretization sensitivity is likely harmless because the stiffness of the grains has very little influence on the
response. One can distinguish between the two regimes through the value of the compressibility number
=
p
d
k
n
,
24
where
p
is the pressure, and
d
is a characteristic dimension of the particle, and
k
n
is the normal contact stiffness.
24
The first strategy proposed in the literature to address LS-DEM’s discretization sensitivity is a single-node approach.
19
It is based on replacing the multiple contact interactions at a contact region between two grains
*
with a single contact
interaction in which a single representative penetration—either the deepest penetration or the average one—is used in
the contact law. Taking the deepest penetration was found to be the more robust approach
19
and it has since been the
more commonly-used one.
19-22
The main problems with the single-node approach, as illustrated and discussed in detail
in Reference
23
, are that it is possible to have the same contact force for two very different penetration states, that the
contact force is prone to discontinuities (even within a single time step), and that it does not conserve elastic energy.
Another problem is that, since only one penetration is used in a contact between two grains, it fails when the grains have
more than one contact region, as is typical with non-convex grains.
The recently-developed volume-interacting level-set discrete-element-method(VLS-DEM)
23
has been presented as an
approach that solves LS-DEM’s discretization-related issues, mainly the non-conservation of elastic energy. It abandons
LS-DEM’s surface discretization altogether, basing instead the contact law on the volume of penetration, that is, the inter-
section volume of penetrating grains (i.e., level-sets). Aside from conserving elastic energy, it is likely that this approach
would indeed avoid LS-DEM’s discretization-sensitivity issue, even though this has not been extensively looked into in
Reference
23
. The two main shortcomings of the VLS-DEM approach are that it is 2–4 times slower than LS-DEM
23
(due
to the computationally expansive search of the intersection volume) and that the contact law deviates from the classical
continuumcontactdescription.Thisdeviationleadstosomeinternalcontradictionsintheformulationofthevolume-type
contact law and it makes it necessary to rely on assumptions regarding the volumetric normal stiffness and the normal
vector, assumptions whose validity has not been tested.
The objective of this study is to present a contact formulation for LS-DEM that solves the discretization-sensitivity in
a straight-forward, methodologically-rigorous, and computationally-efficient manner. This formulation is based on the
classical continuum contact description and its application implies minimal, yet essential, modifications to the original
LS-DEM, making it not only effective, but also elegant.
Next, in Section
2
, we examine LS-DEM’s discretization-sensitivity on a simple two-sphere benchmark and pinpoint
thereasonforitinLS-DEM’soriginalcontactformulation.Wethenpresenttheproposedadaptedcontactformulationand
show that is converges in the benchmark example. In Section
3
, we examine the convergence properties of the adapted
formulation for several highly-compressible cases with multiple grains, and discuss their significance.
2
METHODOLOGY
2.1
Benchmark divergence example
We illustrate the divergence of the original LS-DEM formulation in the two-sphere compression case depicted in
Figure
1A
, with two identical grains. Both grains are initially fixed in space such that Grain 0’s rightmost point is just
10970207, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7400 by California Inst of Technology, Wiley Online Library on [12/07/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
FELDFOGEL et al.
3of12
(A)(B)
(D)
(C)
FIGURE 1
Benchmark example for divergence in the original LS-DEM with surface discretization refinement under
displacement-controlled loading. (A) Configuration of the two-sphere central-compression benchmark (B) two different triangle-mesh-based
surface discretizations (C) force-displacement (
P
) curves with varying number of surface nodes (Nnodes) and a fixed
k
n
=
1GN/mm(D)
P
curves with varying
k
n
’s and a fixed Nnodes
8000.
touching Grain 1’s leftmost one. Then, Grain 1 is displaced leftwards at a constant rate along the line connecting their
centers, such that its displacement
is equal to the penetration between the grains. The response parameter of interest is
the contact force
P
.
To examine the effect of surface discretization (SD) refinement on
P
, we considered four progressively refined surface
discretizations
, see examples in Figure
1B
, and examined the corresponding
P
curves. We found the
P
curves
diverge with SD refinement. We also found that SD refinement is equivalent to increasing the normal contact stiffness
parameter
k
n
,seeFigure
1D
, an observation we explain later in Section
2.2
. The divergent behavior observed in Figure
1C
underscorestheneedtoadaptLS-DEMsothatitbecomesdiscretization-convergent,achallengethatdefinestheobjective
of this study.
2.2
Why exactly does the original LS-DEM diverge?
The root cause of LS-DEM’s divergent behavior is that the inter-grain penetration stiffness, which controls the global
strength and stiffness, scales linearly with the number of surface nodes. In ordinary DEM, the grains are spherical and
each contact comprises a single penetration point
a
,seeFigure
2A
. In LS-DEM, due to the surface discretization, there
are many contact nodes between two contacting grains, see Figure
2B
. At contact node
a
, the normal force acting on grain
i
by grain
j
is, similarly to ordinary DEM:
F
i
n
,
a
=
k
n
d
j
,
i
a
̂
n
j
,
i
a
,
(1)
where
k
n
is the normal penalty parameter with dimensions (force/penetration), and where
d
j
,
i
a
and
̂
n
j
,
i
a
denote the
penetration depth and the unit normal vector to grain
j
at penetration point
a
, respectively.
The resultant normal force
F
i
n
acting on grain
i
by grain
j
is the vector sum of nodal forces over the
N
contact nodes:
F
i
n
=
N
a
=
1
F
i
n
,
a
=
N
a
=
1
k
n
d
j
,
i
a
̂
n
j
,
i
a
,
(2)
where
N
is the total number of surface nodes of grain
i
in contact with grain
j
.
10970207, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7400 by California Inst of Technology, Wiley Online Library on [12/07/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
4of12
FELDFOGEL et al.
(A)
(B)
(C)
(D)
FIGURE 2
Illustration of differences between the contact descriptions in DEM, original LS-DEM, and adapted LS-DEM. (A) DEM
contact approach: Penetration of grains
i
and
j
is defined by the position of the single penetration point a. (B) The original LS-DEM contact
approach: Surface discretization nodes are indicated by black dots. Penetration
d
j
,
i
a
of grain
j
into grain
i
at discretization node a is shown as a
dashed line. The resultant normal force
F
i
n
is obtained as the sum of nodal forces
F
i
n
,
a
. They are indicated by full lines and a dashed one,
respectively. (C) In a continuum-based contact description, the resultant normal force
F
i
,
c
n
is obtained by integration of the normal tractions
f
i
,
c
n
over the contact area s. (D) The adapted LS-DEM formulation is based on the continuum description. Each node has a tributary area
A
a
,
and the resultant force
F
i
,
n
is obtained as a sum of nodal forces
F
i
,
n
,
a
, analogously to the continuum-based integral. The penetrations are
grossly exaggerated for illustrative purposes.
By defining the average penetration
d
j
,
i
a
, we find that the total grain penetration stiffness, which is defined as
k
grain
|
F
i
n
|
d
j
,
i
a
=
N
a
=
1
|
k
n
d
j
,
i
a
̂
n
j
,
i
a
|
N
a
=
1
d
j
,
i
a
N
=
N
k
n
,
(3)
is linear in both
N
and
k
n
.
The distinction between the grain penetration stiffness
k
grain
, which refers to interaction between grains as a whole,
and the normal penalty parameter
k
n
, which refers to contact interactions at the node level, does not exist in ordinary
DEM. There,
N
=
1and
k
grain
and
k
n
are therefore one and the same.
10970207, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7400 by California Inst of Technology, Wiley Online Library on [12/07/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
FELDFOGEL et al.
5of12
The linearity of
k
grain
with
k
n
in Equation (
3
) is reflected in Figure
1D
, and is common in DEM. It is useful because
it allows to relate
k
n
to the Young’s modulus, either through calibration, for example, References
7
and
25
, or through
a closed-form expression.
16,17
In contrast, the linearity of
k
grain
with
N
is unique to the original LS-DEM formulation
and it is problematic because, as the surface discretization is refined,
N
and
k
grain
increase proportionally, leading to the
divergence observed in Figure
1
.
2.3
How we eliminate LS-DEM’s divergence
The key to fixing the divergence of the original LS-DEM formulation is to eliminate the discretization dependence of the
penetrationstiffnessexpressedinEquation(
3
).Todoso,weadoptthecontinuum-basedapproachillustratedinFigure
2C
,
wherein the contact forces are integrals of surface tractions over the contact area
s
, rather than sums of nodal forces.
The dimensions of the continuum-based penalty parameter
k
n
are those of an elastic foundation (traction/penetration),
differently from
k
n
’s (force/penetration). This approach provides the resultant normal force
F
i
,
c
n
and the commensurate
penetrationstiffness
k
grain
withthenecessarylinkagetoandboundingbythefinitecontactarea
s
.Theoriginalformulation
expressed in Equation (
2
) lacks this bounding property.
The continuum-based resultant normal force
F
i
,
c
n
reads:
F
i
,
c
n
=
s
d
F
i
,
c
n
,
(4)
where the superscript
,
c
denotes the continuum-based approach,
s
denotes the contact area between the grains, and the
elemental
d
F
i
,
c
n
reads:
d
F
i
,
c
n
=
f
i
,
c
n
ds
=
k
n
d
j
,
i
̂
n
j
,
i
ds
,
(5)
where
f
i
,
c
n
are the normal tractions,
d
j
,
i
the normal penetrations, and
̂
n
j
,
i
the surface normals.
WeevaluatetheintegralinEquation(
4
)numericallyasasumofnodalcontributionsatcontactingnodes,asillustrated
in Figure
2D
. The resultant normal force
F
i
,
n
in the adapted LS-DEM therefore reads:
F
i
,
n
=
N
a
=
1
F
i
,
n
,
a
,
(6)
where the superscript
,
denotes the adapted formulation. The nodal normal force
F
i
,
n
,
a
, which is the finite analog of
d
F
i
,
c
n
from Equation (
5
), reads:
F
i
,
n
,
a
=
k
n
d
j
,
i
a
̂
n
j
,
i
a
A
a
,
(7)
where
A
a
is the tributary area of node
a
,seeFigure
2D
.
Assuming that
A
a
is taken equal to the total surface area of the grain divided by the total number of surface nodes
§
,
the penetration stiffness
k
grain
,
in the adapted model then reads:
k
grain
,
N
a
=
1
|
F
i
,
n
,
a
|
d
j
,
i
a
=
N
a
=
1
k
n
A
a
=
k
n
N
A
a
.
(8)
The “corrective” factor
A
a
and the elastic-foundation nature of
k
n
are the two keys to resolving the divergence issue.
Increase in
N
with SD refinement is off-set (“corrected”) by a proportional decrease in
A
a
, so that the sensitivity of
k
grain
,
on the SD becomes minor, and the global response converges upon SD refinement, as will be shown in the next section.
The tangential contact formulation requires a treatment similar to that of the normal contact formulation, for similar
reasons.ComparingEquations(
1
)and(
7
),themodificationofthenormalnodalforceintheadaptedformulationamounts
to replacing
k
n
with
k
n
A
a
. The modification of the nodal tangential force in the adapted formulation is analogous to
this replacement. In the original LS-DEM (see eq. (15) in Reference
7
), the elastic tangential force at node
a
,whichis
the basis for the commonly-used Coulomb friction law, is based on the product
k
s
Δ
S
a
,where
k
s
(force/displacement)
10970207, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7400 by California Inst of Technology, Wiley Online Library on [12/07/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
6of12
FELDFOGEL et al.
is the tangential stiffness analogous to
k
n
and where
Δ
S
a
is the (incremental) tangential slip vector analogous to the
penetration vector
d
j
,
i
a
̂
n
j
,
i
a
. Hence, the modification of the nodal shear force amounts to replacing
k
s
with
k
s
A
a
,where
the distributed tangential stiffness
k
s
has dimensions (traction/penetration).
The major advantage of the adapted contact formulation compared to the volume-based approach
23
and the
single-node approach
19
is that it is derived directly from the well-established continuum contact description based on
surface tractions. As such, it account for the spatial variation of contact tractions across the contact area, preserves the
classical significance and dimensions of the penetration stiffness (traction/penetration), and obtains the normal and tan-
gential vectors directly from the geometrical description, without any associated assumptions. All these elements make
the proposed approach not only methodologically sound, but also numerically smooth, stable, and economical.
2.4
Adapted LS-DEM formulation
With the modifications to the contact formulation described above, the adapted LS-DEM’s formulation is neatly obtained
fromtheoriginalonebyreplacing
k
n
and
k
s
inEquations(5)and(9)fromReference
7
with
k
n
A
a
and
k
s
A
a
,respectively.
Regarding the computational aspect, it is pointed out that the additional computational cost in the adapted formulation
is effectively zero, consisting merely of calculating
A
a
once for each node in the pre-processing stage and appending it to
k
n
and
k
s
whenever the node is in contact. This speaks to the computational efficiency of our approach.
3
RESULTS AND DISCUSSION
The convergence properties of the adapted LS-DEM in different deformation-governed cases are demonstrated and
discussed next.
3.1
Two-sphere central-compression
We begin exploring the convergence of the adapted LS-DEM formulation in the illustrative two-sphere benchmark for
which the original LS-DEM diverged, see Figure
1
. The force-displacement with progressively refined SD’s in Figure
3A
show that the adapted LS-DEM converges upon SD refinement. Further, for a given surface discretization, the response
scaleslinearlywith
k
n
,seeFigure
3
,reflectingthelinearityof
k
grain
,
with
k
n
inaccordancewithEquation(
8
).Theseresults
indicate that the cause of divergence in the original LS-DEM formulation has been correctly identified and eliminated as
described in Section
2
.
(A)(B)
FIGURE 3
Force-displacement curves for the two-sphere central-compression benchmark case with the adapted LS-DEM
formulation: (A)
P
curves for a fixed
k
n
=
1 GPa/mm and five progressively refined surface discretizations. (B)
P
curves for a fixed
Nnodes
8000 and varying
k
n
.
10970207, 2024, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7400 by California Inst of Technology, Wiley Online Library on [12/07/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License