Calculation of the incremental stress-strain relation of a polygonal packing
F. Alonso-Marroquin and H. J. Herrmann
ICA1, University of Stuttgart, Pfaffenwaldring 27, 70569 Stuttgart, Germany
共Received 25 March 2002; published 20 August 2002兲
The constitutive relation of the quasistatic deformation on two-dimensional packed samples of polygons is
calculated using molecular dynamics simulations. The stress values at which the system remains stable are
bounded by a failure surface, which shows a power law dependence on the pressure. Below the failure surface,
nonlinear elasticity and plastic deformation are obtained, which are evaluated in the framework of the incre-
mental linear theory. The results show that the stiffness tensor can be directly related to the microcontact
rearrangements. The plasticity obeys a nonassociated flow rule with a plastic limit surface that does not agree
with the failure surface.
DOI: 10.1103/PhysRevE.66.021301 PACS number共s兲: 45.70.Cc, 83.10.Gr, 61.43.Bn, 83.80.Nb
I. INTRODUCTION
The nonlinear and irreversible behavior of soils has been
described by different constitutive theories 关1,2兴. Here the
stress-strain relation is postulated using a certain number of
material parameters which are measured in experimental
tests. These continuous theories have been used for many
geotechnical applications. Excavations, foundations, and
landslides are some few examples of these applications.
Recently, the investigation of soils at the grain scale has
become possible using numerical simulations 关3兴. They evi-
dence that the stress is transmitted through a heterogeneous
network of interparticle contacts 关4兴. The geometric change
of this network during deformation reveals a structural an-
isotropy induced by shearing 关5兴. Although these results pro-
vide valuable insights into the behavior of soils, few issues
are given to derive the continuous relations from the discrete
models.
In this paper the stress-strain relation of a two-
dimensional discrete model is calculated explicitly using nu-
merical simulations. An internal variable is included in this
continuous relation, which is related to the anisotropy of the
contact network. The results show that it is possible to char-
acterize the mechanical behavior of soils at the macroscopic
scale using particle models. In effect, we demonstrate that
simple mechanical laws at the grain level are able to repro-
duce the complex behavior of the deformation of soils.
Usually, disks or spheres are used in the modeling of
granular materials. The simplicity of their geometry allows
one to reduce the computer time of calculations, but they do
not take into account the diversity of the shapes of the grains
in the realistic materials. A more detailed description is pre-
sented here by using randomly generated convex polygons.
As presented by Tillemans and Herrmann 关6兴, the interaction
between the polygons could be handled by letting the poly-
gons interpenetrate each other and calculating the force as a
function of their overlap. This approach has been success-
fully applied to model different processes, such as fragmen-
tation 关7,8兴, damage 关9兴, strain localization, and earthquakes
关6兴.
This paper is organized as follows. A suitable contact
force law is introduced in Sec. II A, which attempts to com-
bine the Hertz contact law with the Coulomb friction crite-
rion. The boundary condition is introduced in Sec. II B by a
flexible membrane that surrounds the sample. The modeling
with such a membrane is very advantageous since it allows
one to implement a stress-controlled loading without any re-
striction in the deformation of the boundary. The strain re-
sponse is calculated in Sec. III for different stress increments
applied on identically generated samples. The results are dis-
cussed in Sec. IV in the framework of the theory of elasto-
plasticity.
II. MODEL
The polygons of this model are generated using a simple
version of the Voronoi tessellation: First, we set a random
point in each cell of a regular square lattice, then each poly-
gon is constructed assigning to each point that part of the
plane that is nearer to it than to any other point. Each poly-
gon is subjected to interparticle contact forces and boundary
forces. They are inserted in Newton’s equation of motion as
we explain below.
A. Contact force
Usually, the interaction between two solid bodies in con-
tact is described by a force applied on the flattened contact
surface between them. Given two polygons in contact, such
surface is obtained from the geometrical construction shown
in Fig. 1. The points C
1
and C
2
result from the intersection
between the edges of the polygons. The contact surface is
taken as the segment that lies between those points. The
vector S
ជ
⫽ C
1
C
2
ជ
defines an intrinsic coordinate system at the
contact (t
ˆ
,n
ˆ
), where t
ˆ
⫽ S
ជ
/
兩
S
ជ
兩
and n
ˆ
is perpendicular to it.
The deformation length is given by
␦
⫽ a/
兩
S
ជ
兩
, where a is the
overlap area between the polygons. ᐉ
ជ
is the branch vector,
which connects the center of mass of the polygon to the point
of application of the contact force, which is supposed to be
the center of mass of the overlap area.
The normal elastic force is taken proportional to the de-
formation length as f
n
e
⫽ k
n
␦
; the tangential force is calcu-
lated from the simplified Coulomb friction law with a single
friction coefficient
s
⫽
d
⫽
. Here
s
is the static and
d
the dynamic friction coefficient. This tangential force is
PHYSICAL REVIEW E 66, 021301 共2002兲
1063-651X/2002/66共2兲/021301共11兲/$20.00 ©2002 The American Physical Society66 021301-1
implemented by an elastic spring f
t
e
⫽⫺k
t
, where
grows
linearly with the tangential displacement of the contact,
whenever
兩
f
t
e
兩
⬍
f
n
e
. We used the straightforward calcula-
tion of
proposed by Brendel 关10兴,
共
t
兲
⫽
冕
0
t
v
t
共
t
⬘
兲
⌰„
兩
f
t
e
共
t
⬘
兲
兩
⫺
f
n
e
共
t
⬘
兲
…dt
⬘
, 共1兲
where ⌰ is the Heaviside function and
v
ជ
is the relative ve-
locity at the contact, which depends on the linear velocity
v
ជ
i
and angular velocity
ជ
i
of the particles in contact according
to
v
ជ
⫽
v
ជ
i
⫺
v
ជ
j
⫺
ជ
i
⫻ ᐉ
ជ
i
.⫹
ជ
j
⫻ ᐉ
ជ
j
. 共2兲
B. Boundary forces
Let us now discuss how to apply the stress on the sample.
One way to do that would be to apply a perpendicular force
on each edge of the polygons belonging to the external con-
tour of the sample. Actually, this does not work because this
force will act on all the fjords of the boundary. It produces an
uncontrollable growth of cracks that with time ends up de-
stroying the sample. Thus, it is necessary to introduce a flex-
ible membrane in order to restrict the boundary points that
are subjected to the external stress.
The algorithm to identify the boundary is rather simple.
The lowest vertex p from all the polygons of the sample is
chosen as the first point of the boundary list b
1
. In Fig. 2 P is
the polygon that contains p, and q苸 P艚Q is the first inter-
section point between the polygons P and Q in counterclock-
wise orientation with respect to p. Starting from p, the verti-
ces of P in counterclockwise orientation are included in the
boundary list until q is reached. Next, q is included in the
boundary list. Then, the vertices of Q between q and the next
intersection point r苸Q艚R in the counterclockwise orienta-
tion are included into the list. The same procedure is applied
until one reaches the lowest vertex p again. This is a very
fast algorithm, because it only makes use of the contact
points between the polygons, which are previously calcu-
lated to obtain the contact force.
The set of points that are in contact with the membrane
are selected using a recursive algorithm. It is initialized with
the vertices of the smallest convex polygon that encloses the
boundary 共see Fig. 3兲. The lowest point of the boundary is
selected as the first vertex of the polygon m
1
⫽ b
1
. The sec-
ond one m
2
is the boundary point b
i
that minimizes the angle
⬔ (b
1
b
i
ជ
) with respect to the horizontal. The third one m
3
is
the boundary point b
i
such that the angle ⬔ (m
2
b
i
ជ
,m
1
m
2
ជ
)is
minimal. The algorithm is recursively applied until the low-
est vertex m
1
is reached again.
The points of the boundary are iteratively included in the
list m
i
using the bending criterion proposed by Åstro”m et al.
关11兴: For each pair of consecutive vertices of the membrane
m
i
⫽ b
i
and m
i⫹ 1
⫽ b
j
we choose that point from the subset
兵
b
k
其
i⭐k⭐ j
that maximizes the bending angle
b
⫽ ⬔ (b
k
b
i
ជ
,b
k
b
j
ជ
). This point is included into the list, when-
ever
b
⭓
th
. Here
th
is a threshold angle for bending. This
algorithm is repeatedly applied until there are no more points
satisfying the bending condition.
The final result gives a set of segments
兵
m
i
m
ជ
i⫹ 1
其
lying on
the boundary of the sample. In order to apply the boundary
forces, those segments are divided into two groups: A-type
segments are those that coincide with an edge of a boundary
polygon; B-type segments connect the vertices of two differ-
ent boundary polygons.
On each segment of the membrane m
i
m
ជ
i⫹ 1
a force f
i
⫽
i
N
i
is applied, where
i
is the local stress and N
i
is the
90° counterclockwise rotation of m
i
m
ជ
i⫹ 1
. This force is trans-
mitted to the polygons in contact with it: if the segment is of
A type, this force is applied in its midpoint; if the segment is
of B type, half of the force is applied at each one of the
vertices connected by this segment.
C. Molecular dynamics simulation
Before we implement the numerical solution of Newton’s
equations it is convenient to make a dimensional analysis of
FIG. 1. Contact surface as defined from the geometry of over-
lap.
FIG. 2. Algorithm used to find the boundary.
F. ALONSO-MARROQUIN AND H. J. HERRMANN PHYSICAL REVIEW E 66, 021301 共2002兲
021301-2
the parameters. In such way we can keep the scale invariance
of the model and reduce the parameters to a minimum of
dimensionless constants. All the polygons are supposed to
have the same density. The mass m
i
of each polygon is mea-
sured in units of the mean mass m
0
of the Voronoi tessella-
tion. The time is measured in fractions of the total loading
time t
0
. The evolution of the position x
ជ
i
and the orientation
i
of the i
th
polygon is governed by the equations of motion,
2
m
i
x
ជ
¨
i
⫹
兺
c
f
i
c
ជ
⫹
兺
c
b
i
b
k
n
f
ជ
i
b
⫽ 0
ជ
,
2
I
i
¨
i
⫹
兺
c
ᐉ
ជ
i
c
⫻ f
i
c
ជ
⫹
兺
c
b
i
b
k
n
ᐉ
ជ
i
b
⫻ f
ជ
i
b
⫽ 0
ជ
. 共3兲
The sums go over all those particles and boundary seg-
ments that are in contact with the i
th
polygon. The interpar-
ticle contact forces f
i
c
ជ
and boundary forces f
i
b
ជ
are given by
f
i
c
ជ
⫽
共
␦
i
c
⫹
␥
m
v
n
c
兲
n
ˆ
i
c
⫹
共
i
c
⫺
␥
m
v
t
c
兲
t
ˆ
i
c
,
f
i
b
ជ
⫽ N
i
b
ជ
⫺
␥
m
i
v
i
ជ
. 共4兲
Here
␦
i
c
and
i
c
denote the deformation length and the
tangential displacement of the contact, which were defined in
Sec. II A;
i
b
is the stress applied on the boundary segment
T
i
b
, defined in Sec. II B. Artificial viscous terms must be
included in Eq. 共4兲 to keep the stability of the numerical
solution and reduce the acoustic waves generated during the
loading process.
v
ជ
c
denotes the relative velocity at the con-
tact 关Eq. 共2兲兴 and m⫽ (1/m
i
⫹ 1/m
j
)
⫺ 1
the effective mass of
the two polygons in contact.
There are four microscopic parameters in the model: the
viscosity
␥
, the ratio ⫽ t
s
/t
o
between the characteristic pe-
riod of oscillation t
s
⫽
冑
k
n
/m
0
and the loading time t
0
, the
friction coefficient
, and the ratio
⫽ k
t
/k
n
between the
tangential k
t
and normal k
n
stiffness of the interparticle con-
FIG. 3. Membrane obtained with threshold bending angle
th
⫽
,3
/4,
/2, and
/4. The first one corresponds to the minimum convex
polygon that encloses the sample.
CALCULATION OF THE INCREMENTAL STRESS-...
PHYSICAL REVIEW E 66, 021301 共2002兲
021301-3
tacts. The viscosity factor
␥
is related to the normal restitu-
tion coefficient 关12兴. It was taken large enough to have a high
dissipation, but not too large to keep the numerical stability
of the method. The ratio was chosen small enough in order
to avoid rate dependence in the strain response, correspond-
ing to the quasistatic approximation. Technically, it is done
by looking for the value of such that a reduction of it by
half makes a change of the strain response less than 5%.
The two parameters
and
determine the constitutive
response of the system. For example, the micromechanical
analysis of the strain response shows that the Young’s modu-
lus and Poisson’s ratio depend on
关13兴. On the other hand,
can be directly related to the friction angle of the material
关14兴. Although the study of the dependence of the constitu-
tive response on those parameters is an important point, such
quantities have been kept fixed in this work.
The boundary conditions yield more dimensional param-
eters. The initial height H
0
and width W
0
of the sample, and
the characteristic length ᐉ
0
of the polygons define two geo-
metrical parameters, which are the shape ratio W
0
/H
0
and
the granularity ᐉ
0
/H
0
of the sample 共see Table I兲.
In order to keep overlaps much smaller than the charac-
teristic area of the polygons, the ratio
i
/k
n
between the
stress applied on the membrane and the stiffness of the con-
tacts is restricted to small values. This was implemented by
fixing the contact stiffness to a value close to the experimen-
tal granular stiffness k
n
⫽ 160 MPa. Then the stress is cho-
sen in such a way that it does not exceed 1% of this value.
III. STRESS-STRAIN CALCULATION
A. Theoretical background
The macroscopic state of the system is characterized by
the stress tensor and the void ratio e. The area fraction of
voids in the sample defines the void ratio. Initially e
0
⫽ 0 due
to the Voronoi tessellation used. The stress controlled test
was restricted to stress states without off-diagonal compo-
nents. The diagonal components, the axial
␦
1
and lateral
␦
3
stress, define the stress vector,
˜
⫽
冋
p
q
册
⫽
1
2
冋
␦
1
⫹
␦
3
␦
1
⫺
␦
3
册
, 共5兲
where p and q are the pressure and the shear stress. The
domain of admissible stresses is bounded by the failure sur-
face. When the system reaches this surface it becomes un-
stable and fails.
Before failure, the constitutive behavior can be obtained
performing small changes in the stress and evaluating the
resultant deformation. An infinitesimal change of the stress
vector d
˜
produces an infinitesimal deformation of the
sample, which is given by a change of height dH and width
dW. This defines the axial strain d
⑀
1
⫽ dH/H and lateral
strain d
⑀
3
⫽ dW/W increments. The volumetric strain d
⑀
v
and the shear strain d
⑀
␥
increments define the incremental
strain vector,
d
⑀
˜
⫽
冋
d
⑀
v
d
⑀
␥
册
⫽
冋
d
⑀
1
⫹ d
⑀
3
d
⑀
1
⫺ d
⑀
3
册
. 共6兲
Each state of the sample is related to a single point in the
stress space, and the quasistatic evolution of the system is
represented by the movement of this point in the stress space.
The constitutive relation is formulated taking the incremental
strain as a function of the incremental stress and the stress
state
d
⑀
˜
⫽ F
共
d
˜
,
˜
兲
. 共7兲
If there is no rate dependence in the constitutive equation,
F(d
˜
) is an homogeneous function of degree 1. In this case,
the application of the Euler identity 关15兴 shows that Eq. 共7兲
can be reduced to
d
⑀
˜
⫽ M
共
ˆ
,
˜
兲
d
˜
. 共8兲
Where
ˆ
is the unitary vector defining a specific direction in
the stress space,
ˆ
⫽
d
˜
兩
d
˜
兩
⬅
冋
cos
sin
册
,
兩
d
˜
兩
⫽
冑
dp
2
⫹ dq
2
. 共9兲
The constitutive relation results from the calculation of
d
⑀
˜
(
), where each value of
is related to a particular mode
of loading. Some special modes are listed in Table II.
The relation 共8兲 has been proposed by Darve 关15兴 and it
contains all the possible constitutive equations. In order to
interpret our particular results, it is convenient to make some
approximations: First, if the load increments are taken small
enough, the tensor M(
) can be supposed to be linear in
TABLE I. Dimensionless variables.
Variable Ratio Default value
Viscosity
␥
0.1
Friction coefficient
0.25
Time ratio ⫽t
s
/t
o
8.0⫻ 10
⫺ 4
Stiffness ratio
⫽ k
t
/k
n
0.33
Granularity ᐉ
0
/H
0
0.1
Shape ratio W
0
/H
0
1.0
Bending angle
th
0.25
TABLE II. Principal modes of loading according to the orienta-
tion of
ˆ
.
0° Isotropic compression dp⬎0 dq⫽ 0
45° Axial loading d
1
⬎ 0 d
3
⫽ 0
90° Pure shear dp⫽0 dq⬎ 0
135° Lateral loading d
1
⫽ 0 d
3
⬎ 0
180° Isotropic expansion dp⬍0 dq⫽ 0
225° Axial stretching d
1
⬍ 0 d
3
⫽ 0
270° Pure shear dp⫽0 dq⬍ 0
315° Lateral stretching d
1
⫽ 0 d
3
⬍ 0
F. ALONSO-MARROQUIN AND H. J. HERRMANN PHYSICAL REVIEW E 66, 021301 共2002兲
021301-4
each stress direction. Then, we assume that the strain can be
separated in an elastic 共recoverable兲 and a plastic 共unrecov-
erable兲 component,
d
⑀
˜
⫽ d
⑀
˜
e
⫹ d
⑀
˜
p
, 共10兲
d
⑀
˜
e
⫽ D
共
˜
兲
d
˜
, 共11兲
d
⑀
˜
p
⫽ J
共
,
˜
兲
d
˜
. 共12兲
Here, D
⫺ 1
defines the stiffness tensor, and J⫽ M⫺ D the
flow rule of plasticity, which results from the calculation of
d
⑀
˜
e
(
) and d
⑀
˜
p
(
).
B. The method
The numerical method presented here was proposed by
Bardet 关16兴. It allows one to find the elastic d
⑀
˜
e
and plastic
d
⑀
˜
p
components of the strain as functions of the stress state
˜
and the stress direction
. Figure 4 shows the three steps
of the procedure.
共1兲 The sample is driven to the stress state
˜
. First, it is
isotropically compressed until it reaches the stress value
␦
1
⫽
␦
3
⫽ p⫺ q. Next, it is subjected to axial loading, in order to
increase the axial stress
␦
1
to p⫹ q 共see Fig. 5兲. When the
stress state
˜
⫽
关
pq
兴
T
is reached, (A
T
being the transpose of
A) the sample is allowed to relax.
共2兲 Loading the sample from
˜
to
˜
⫹ d
˜
the strain incre-
ment d
⑀
˜
is obtained. This procedure is implemented choos-
ing different stress directions according to Eq. 共9兲. Here the
stress modulus is fixed to
兩
d
˜
兩
⫽ 10
⫺ 4
p.
共3兲 The sample is unloaded until the original stress state
˜
is reached. Then one finds a remaining strain d
⑀
˜
p
that corre-
sponds to the plastic component of the incremental strain.
Since the stress increments are taken small enough, the un-
loaded stress-strain path is practically elastic. Thus, the dif-
ference d
⑀
˜
e
⫽ d
⑀
˜
⫺ d
⑀
˜
p
represents the elastic component of
the strain.
One could be concerned about the dependence of the
strain response on the way how the stress state is reached.
We found that there is not remarkable dependence of the
strain response on the stress path, whenever the stress com-
ponents are quasistatic and monotonically increased. Other-
wise, a strong reduction in the plastic component of the
strain is observed. In fact, when the plastic response is cal-
culated after the sample is unloaded, the plasticity is smaller
than that one calculated after a monotonic load. Furthermore,
there is no plastic component in the strain response when
elastic waves are previously generated in the sample. Those
memory effects suggest that the plastic component of the
strain depends on the history of the deformation, and is kept
unchanged only if the sample is subjected to quasistatic and
monotonic loading.
Figure 6 shows the load-unload paths and the correspond-
ing strain response. They were taken from a stress state with
q⫽ 0.5p. The end of the load paths in the stress space map
into a strain envelope response d
⑀
˜
(
) in the strain space.
Likewise, the end of the unload paths map into a plastic
envelope response d
⑀
˜
p
(
). The yield direction
can be
found from this response, as the direction in the stress space
where the plastic response is maximal. The flow rule can be
obtained taking the direction
of the maximal plastic re-
sponse in the strain space. These angles do not agree, which
reveals the necessity to analyze this behavior in the frame-
work of the nonassociated theory of plasticity 共see Sec.
IV C兲.
IV. CONSTITUTIVE RELATION
Figure 7 summarizes the global elastoplastic behavior.
The elastic response, calculated from Eq. 共10兲, has a centered
ellipse as envelope response. This can be related to the mi-
crocontact structure using a local linear relation in each point
of the stress space 共see Sec. IV B兲. The solid line represents
FIG. 4. Procedure to obtain the constitutive behavior: 共1兲 The
sample is driven to the stress state
˜
, with pressure p and shear
stress q. 共2兲 It is loaded from
˜
to
˜
⫹ d
˜
. 共3兲 It is unloaded to the
original stress state
˜
.
FIG. 5. Axial stress
1
⫽ p⫹ q and lateral stress
3
⫽ p⫺ q in a
stress controlled test. They are applied on the boundary of the tes-
sellated sample of polygons.
CALCULATION OF THE INCREMENTAL STRESS-...
PHYSICAL REVIEW E 66, 021301 共2002兲
021301-5
