1 INTRODUCTION

In a plane strain compression on dense sands, homo-

geneous deformation first takes place. Near peak
stress, the deformation suddenly localizes into thin
zones (called shear bands), and the stress drops to a
residual stress state. In soil mechanics, this mode of
failure is very important as well as in engineering
design problems, where one is interested in ultimate
and residual bearing capacities of the analyzed struc-
tures.
To study shear bands, Hill (1962) Mandel (1964)
and Rudnicki & Rice (1975) analyzed the occurrence
and inclination of shear bands. But, their constitu-
tive equations are function only on the first gradient
of displacement then the shear band thickness is un-
determined. Mhlhaus & Vardoulakis used the
Cosserat theory to obtain the shear band thickness as

1
Ingeniera de Diseo INGETEC Ingeniera y Diseo S.A. Co-
rreo electrnico: anagomez@ingetec.com.co
2
Director del Departamento de Ingeniera Civil y Ambiental,
Universidad de los Andes. Correo electrnico: aliz-
cano@uniandes.edu.co
a material length, but Huang & Bauer (2001) and
Tejchman & Gudehus (2001) showed that the thick-
ness of shear bands is not a material constant.
Huang et al.(2002) find that the thickness of shear
bands is proportional to the mean grain diameter and
a friction parameter related with inter-granular fric-
tion and these one is influenced by the density and
stress rate. In all cases, the analyses are developed
in plane strain conditions with constitutive equations
in 2D and supported with quantities of stress and
strain at the boundary of the biaxial specimens but
they did not consider the micro-structure behavior
and this relation with the localization of deformation
process. Several authors has develop studies related
with the micro-structure of shear bands: Bardet &
Proubet (1992) with numerical simulations of granu-
lar assemblies suggest that the particle rotations at
shear bands are induced by macrorotations, Iwashita
& Oda (1998) and Oda & Kazama (1998) showed
with numerical simulation that the microdeformation
mechanism of dilatancy in a shear band is the com-
bination of buckling and rolling of columns, but do
not mention about the implementation algorithm to
boundary conditions. To develop parametric studies
Shear band formation in granular materials: a micromechanical
approach
Ana Mara Gmez Gmez
1

Arcesio Lizcano Pelez
2








RESUMEN: Las bandas de corte como modo de falla de suelos granulares han sido ampliamente estudiadas
sin explicar fsicamente su mecanismo de formacin ni reproducir exactamente un patrn de falla definido.
Se desarrollaron varios algoritmos usando un mtodo de elementos discretos, en los que se busc reproducir
la formacin de una banda de corte nica en una simulacin de un ensayo biaxial sobre un ensamblaje de dis-
cos rgidos, para obtener informacin de esfuerzos en la zona de corte.
ABSTRACT: Shear bands as mode of failure in granular soils has been largely studied without to physically
explain the formation mechanism of that, neither to reproduce exactly a defined failure pattern. Several algo-
rithms were performed with a distinct element method. In this ones, are looking for to reproduce only shear
band in a biaxial simulation on rigid disk assembly, to get information at shear zone.
of the influence of initial density, load velocity and
stiffness at contacts it is necessary first to have an
algorithm that permit to reproduce realistic biaxial
test with particle assemblies. In this paper, two al-
gorithms to reproduce flexible boundaries (like
membranes) are presented. The goal of the work is
to develop a deeper understanding of the microme-
chanics of particle assemblies.

2 DISCRETE MODEL

The discrete model employed to implement the soft
boundary algorithm was the commercially available
code PFC-2D v3.0 (Itasca, 2002). The mathematical
formulation is based upon the solution of Newtons
equations of motion for each particle in a two-
dimensional assembly.
The assumptions of method are: contacts occur at
a point, particles may be represented as rigid disks or
spheres, and particle overlap is allowed at contacts
but these overlaps are small relative to particle size.
The solution include: law of motion, force-
displacement law, as part of calculation cycle, and
contact constitutive model and boundary conditions
that defines the simulation conditions and last ones
are user defined.
2.1 Force-displacement law
Solution of the problem begins with determining the
unit normal vectors of the contacts in the assembly
as:

[ ] [ ]
d
x x
n
A
i
B
i
i

= (1)
where: n
i
is the normal vector; x
i
[A]
and x
i
[B]
are the lo-
cations of particle centers for particles A and B, re-
spectively, and d is the distance between particle
centers. The overlap at contacts, U
n
, is calculated as:

[ ] [ ]
[ ]

+
=
d R
d R R
U
B
B A
n
(2)
In first case is defines the overlap between two ball
and the second one the overlap between ball and
wall;
[ ]
R

is the radius of ball and d is the dis-
tance between two centers of contiguous particles.
The location of the contact point is:

[ ]
[ ] [ ]
( )
[ ] [ ]
( )

+
+
=
i
n B B
i
i
n A A
i C
i
n U R x
n U R x
x
2
1
2
1
(3)
When contact location is known, it is possible to
calculate velocities and invoke the constitutive rela-
tion for contact theory to calculate normal and shear
forces on each particle.
2.2 Motion law
Each new calculation cycle, new particle forces, F
i
,
and moments, M
3
, are used to calculate the new ac-
celeration and angular momentum, and then the
Newtons equations of movement are applied:

) (
i i i
g x m F = & & (4)

i i
H M
&
= (5)

where m is particle velocity, g
i
is acceleration due to
gravity, H
i
is the angular momentum of the particle
and is 0.4 for spheres and 0.5 for disks. The equa-
tions of motion 4 and 5 are discretized with central
finite differences:

( ) ( )
( )
( ) ( )
( )
2 /
3
2 /
3 3
2 / 2 /
1
1
t t t t
t t
i
t t
i i
t
x x
t
x
+
+

=
&
& & & &
(6)

( ) ( )
( ) ( )
t
I
M
t g
m
F
x x
t
t t t t
i
t
i t t
i
t t
i

|
|

\
|
+ =

|
|

\
|
+ + =
+
+
) (
3 2 /
3
2 /
3
) (
2 / 2 /

& &
(7)

( ) ( )
t x x x
t t
i
t
i
t t
i
+ =
+ + 2 / ) (
& (8)

the equations 6 to 8 are respectively, translational
and rotational accelerations, translational and rota-
tional velocities, and displacement.
The critical value of the time step for a stable so-
lution to a DEM problem depends on the minimum
eigenperiod of the entire system. Calculating this
value is computationally very expensive, so PFC-2D
uses a simplified spring and mass logic to calculate
the critical time step (Itasca, 2002):

=
rot
tras
crit
k
I
k
m
t (9)

where m and I are the mass and inertia of particle,
and k
tras
and k
rot
are the translational and rotational
stiffnesses. The minimum of all critical timesteps
computed for all degrees-of-freedom of all particles
is the final critical timestep.
2.3 Contact model
The contact model relates relative displacements
with contact forces via Eqs. 10 and 11 below:

i
n n n
i
n U K F = (10)
s
i
s s
i
U k F = (11)
The normal stiffness, K
n
, is secant stiffness and it re-
lates the total normal force to the total normal dis-
placement; the shear stiffness, k
s
, is a tangent stiff-
ness since it relates the increment of shear
displacement to the increment of shear force. The
contact normal and shear stiffnesses are given by:

[ ] [ ]
[ ] [ ] B
n
A
n
B
n
A
n n
k k
k k
K
+
= (12)

[ ] [ ]
[ ] [ ] B
s
A
s
B
s
A
s s
k k
k k
k
+
= (13)

In PFC-2D program, this model as referred as lin-
ear contact model.

2.4 Boundary conditions
In PFC-2D it is possible to define boundary condi-
tions by two kinds of control entities: velocity-
controlled walls or velocity or force-controlled parti-
cles.
3 PROGRAM ENVIRONMENT

PFC-2D has a series of commands that permit to
create and deal with the model entities as particles
and walls, and a programming language called fish,
with their own language rules and intrinsic func-
tions, and that enables the user to define new func-
tions and to generate his own control boundary con-
ditions to program simulation environments. Itasca
as support codes of fish has some test environments
and models called Augmented fishtank, to develop
a biaxial test, but it consist in a biaxial compression
test that do not allow the shear band formation be-
cause is not possible the lateral deformation of the
specimen. Nevertheless, several algorithms of
fishtank has been used to develop the biaxial simula-
tion with lateral soft boundary.
3.1 Fishtank biaxial environment
The fiishtank biaxial test consists in a set of proce-
dures to: particle generation, set particle properties,
seats the specimen (applying initial stresses), exe-
cute the load of the biaxial test with walls, and final-
ly to get responses from the assembly.
The particle generation is made into a box of four
walls with a radius expansion algorithm that allow
defining the initial porosity of the specimen. First,
with uniform distribution of particle radii and mini-
mum and maximum radii defined, the algorithm
computes the number of particles required inside the
box based on the mean particle radius and the de-
sired porosity, then the particles are generated and
place randomly in the given sample area. The actual
sample porosity is then calculated and the radii of all
particles are increased by the following multiplica-
tion factor:

measured
desired
n
n
m

=
1
1
(14)
where n
desired
and n
measured
are respectively, the de-
sired and measured porosity. Expanding the radii by
this factor will result in overlapping particles and,
thus, large unbalanced forces. The particles are al-
lowed to equilibrate on their own (without applied
stresses). The particle assembly is created using ze-
ro friction and zero shear stiffness between particles.
In this manner, no shear forces are generated be-
tween the particles in this initial state. After equilib-
rium is reached, particles are assigned friction values
and shear stiffnesses. The sample is then allowed to
re-equilibrate. Finally, Figure 1 show the particle
generation process.


(a)

(b)
Figure 1. Particle generation process. (a) Particle assembly af-
ter generating particles at half their final size (b) Particle as-
sembly after equilibrium has reached. (Itasca, 2002)

In the load process, the top and bottom walls act
as loading platens, and the velocities walls are con-
trolled by a servo-mechanism that maintains the re-
quested stresses. Accordingly with fish manual
(Itasca 2002), the biaxial test begins by applying
confining and vertical stresses, next the dimensions
of the specimen at this stage are taken as reference
dimensions to be used in the computation of stresses
and strains during subsequent load phase, the speci-
men is loaded by moving the platens towards one
another at the velocity v
p
. During load phase, the
deviatoric stress,
d
=
y
-
x
, and axial strain are
monitored. The test ends when axial strain reach the
specified limit (
y
<(
y
)
lim
). Figure 2 show a sketch
of load process at the biaxial test.

Figure 2. Sketch of biaxial-test environment (Itasca, 2002).
3.2 Control parameters
The input parameters that control the fish biaxial test
are referred to define the material assembly and to
control the boundary conditions and the end of the
test. These parameters are listed in table 1.

Table 1. Parameters that control the fish biaxial simulation
Parameter Description Unit
Ball density Kg/cm
3

n Gross areal porosity
Ec Ball-ball contact modulus Pa
k
n
/k
s
Ball stiffness ratio
Ball friction coefficient
H Specimen height m
w Specimen width m

x
Lateral wall stiffness-reduction fac-
tor

y
Vertical wall stiffness-reduction fac-
tor

x
t
Target confining stress Pa

y
t
Target vertical stress Pa
Wall-servo tolerance
v
p
Final platen velocity m/s
N
p
Total cycles of platen acceleration
S
p
Accelerate platens in this many stag-
es

(
y
)
lim
Limiting vertical strain

In table above, the wall stiffness-reduction factors
are employed to set the wall normal stiffness as the
factor for the average particle normal stiffness of the
specimen. The wall-servo tolerance is used to con-
trol the wall velocity, such that corresponding veloc-
ity will be zero when |(
t
)/
t
| . Finally, the N
p

and S
p
parameters are used because if the platen ve-
locity v
p
is applied in a single step, the large acceler-
ation will produce inertial forces within the speci-
men that may cause damage. To eliminate these
inertial effects, the platen acceleration is controlled
with the values of S
p
and N
p
. The platen velocity is
adjusted to reach final value of v
p
in a sequence of S
p

stages over total of N
p
cycles.
4 SOFT BOUNDARY ALGORITHM
Two kinds of particle control were implemented in
algorithms to obtain soft boundary control with par-
ticles: velocity-controlled and force-controlled.
4.1 Velocity-controlled boundary
In a first approximation, after initial stress applica-
tion with walls, behind of each wall was constructed
a column of mono-sized particles. Particles in these
columns were controlled with velocities of a similar
procedure of wall servo-control. There were several
difficulties: first, the indentation of boundary mono-
sized particles caused a not realistic interaction be-
assembly
tween soft boundary and specimen (indentation of
boundary particles into the specimen); second, it is
not possible to monitor the confinement pressure for
the indentation of particles and because when one to
fix velocity to boundary particle, also apply a force
of undetermined magnitude; finally, the boundary
particles increase the complexity of model and the
calculation time. Although, shear band not develops,
the displacement field suggests localization (figure
3).

Table 2. Parameters of biaxial simulation with velocity-
controlled soft boundary and complementary information
Parameter Description Value
Ball density 2630Kg/cm
3

n Gross areal porosity 0.14
Ec Ball-ball contact modulus 88x10
9
Pa
k
n
/k
s
Ball stiffness ratio 1.0
Ball friction coefficient 0.5
R
min
Min. radius of ball assembly 7.5x10
-4
m
R
max
/R
min
Ratio max/min radius 1.66
H Specimen height 64.01x10
-3
m
w Specimen width 32x10
-3
m

x
Lateral wall stiffness-reduction
factor
0.01

y
Vertical wall stiffness-
reduction factor
1.0

x
t
Target confining stress 10
6
Pa

y
t
Target vertical stress 10
6
Pa
Wall-servo tolerance 0.01
v
p
Final platen velocity 5.0x10
-2
m/s
N
p
Total cycles of platen accelera-
tion
400
S
p
Accelerate platens in this many
stages
10
(
y
)
lim
Limiting vertical strain -1.5x10
-2

Complementary information

Number of assembly balls 563

Number of boundary balls 260

Radius of boundary balls 2.5x10
-4
m

RAM memory 128Mb

Processor velocity 500 MHz

Time of process 2.5 days

The parameters used in this simulation and other
important information are listed in table 2 above.


(a)

(b)
Figure 3. Biaxial simulation with soft boundary velocity-
controlled. (a) Final displacement field (b) Chain forces in as-
sembly and boundary particle indentation.
4.2 Stress-controlled boundary
The second approximation of soft boundary consists
in to identify and listing the boundary particles of the
specimen, to get his contact forces with the lateral
walls and boundary particle unbalanced force (hori-
zontal component). To this purpose, it is necessary
to generate several arrays to contain and to hand this
information, and several control and monitoring pro-
cedures that lets to maintain the confinement pres-
sure and to prevent damage or macro unbalanced
condition of the specimen.
The procedure includes three steps: to initialize
the membranes (boundary particles), to homogenize
boundary forces to obtain a uniform stress condition
and to control membrane; to monitor variables and
prevents the lost of confinement pressure of the
specimen in load phase.
The calculation time is lesser than the calculation
time of velocity-controlled (1 day in the same pro-
cessor and 6 hours in a 2.8GHz and 1Gb processor)
4.2.1 Initializing the membrane
This procedure is applied after the initial stress
application with walls. This one consists in to get the
ball identification of the balls in contact with lateral
walls by the list of contacts of the each lateral wall,
next, to obtain the corresponding contact forces and
particles unbalanced forces. At this stage I was ob-
serve that observe that the lateral forces on the
boundary of the specimen are unbalanced (in very
small magnitude but sufficient to cause lateral dis-
placement of the specimen). Then, before to apply
the forces on the balls that replace the confinement
of the lateral walls, it is necessary to equilibrate the
magnitudes of the confining forces on the boundary
balls to the total force on the both membranes are
equal.
4.2.2 Homogenization of boundary forces
In the final of last process, the magnitude of the
applied forces on the boundary particles in the same
membrane is very different particle to particle. To
control the relative vertical displacement of the
boundary particles and as the same time to maintain
constant the confinement pressure it is necessary to
modify magnitudes of boundary forces to obtain a
uniform distribution of forces (uniform stress) over
the membrane, but avoided inertial effect that dam-
age the specimen.
On the other hand, to maintain the general equi-
librium condition in specimen, the difference of the
sum of forces over membranes must remain null to
warranty the system equilibrium.
The procedure consists in: to get the relative
height of each particle over each membrane. This
height is the ratio between the external height of a
boundary particle and the specimen height
(h
k
r
=h
k
e
/h) as show in figure 4 where the k-particle
belong to left membrane, then to compute the aver-
age and standard deviation of relative heights to each
membrane; next in each calculation cycle, to com-
pute the average and standard deviation of percent-
age of applied forces on the boundary particles, the
corresponding normalized variable and comparing
this one with average of normalized variable. If the
normalized variable is larger than the average then
the percentage of force of respective particle is re-
duced in a proportion of standard deviation of force
percentages, else, a proportion of the sum of reduc-
tions of force percentages is added to the actual force
percentage of the boundary particle. The homogeni-
zation process ends when the standard deviation of
force percentage is smaller than the standard devia-
tion of relative height in both membranes.
To prevent the lost of local confinement of the
specimen, the boundary particle x-velocity is fixed
in zero.
With this process of homogenization, the chain
forces before and after this process are equal.


Figure 4. Definition of particle external height in a left bounda-
ry particle.

4.2.3 Control and monitoring
To control the confinement force during the axial
load phase in biaxial test, if the assumption of uni-
form distribution of confinement pressure is true, the
percentage of force applied to each boundary particle
must be proportional to his relative height. Then, in
each calculation cycle, it is necessary to compute the
relative heights. But it is possible that inertial ef-
fects cause lost of contacts in some boundary (figure
5). To correct these inertial effect, in each calcula-
tion cycle, is assign zero velocity to the boundary
particles but in the cycle the particle can to change
his velocity for his unbalanced force.
In a simulation with the same control parameters
of table 2, the final configuration and displacement
field was the showed in figure 5.
With the model used was not possible to repro-
duce shear band but I think that it is possible loading
laterally the specimen with the boundary particles af-
ter the particle generation and change the stiffness
and friction of the particles, but it is necessary to
check the internal sample equilibrium, for the iner-
tial effect. I noted that always the minimum unbal-
anced force is around 10
-6
or 10
-7
Pa, then with the
PFC2D simulations the internal perfect equilibrium
never occurs.


(a)

(b)
Figure 5. Biaxial simulation with soft boundary stress-
controlled. (a) Final displacement field (b) Chain forces and fi-
nal assembly configuration

4.2.4 Local fields of stress
To compute the evolution of local fields of stress, it
is possible to implement local areas to obtain aver-
age stress tensors inside them. In each control area,
are obtained the average of unbalanced forces of par-
ticles inside them, and this calculus can be develop
several times during the axial loading process. If the
control areas are defined small enough and are dis-
posed over all specimen area, it is possible to obtain
the local evolution of stress inside the specimen and
to fish the shear band and his stress evolution.
5 CONCLUSIONS
Two soft boundary algorithms were implemented.
The velocity-control boundaries do not allow the
calculation of relevant macroscopic variables be-
cause the velocities on the boundary particles cause
undetermined forces. The stress-controlled bounda-
ries, allow the calculation of macroscopic variables
but the correct control of confinement pressure was
elusive. The inertial effects on PFC2D are difficult
to control. I think that it is possible to reproduce a
shear band loading laterally the specimen with the
boundary particles after the particle generation and
change the stiffness and friction of the particles, but
it is necessary to check the internal sample equilibri-
um, for the inertial effect caused by the walls during
radius expansion process.
6 REFERENCES
Bardet, J. P. & Proubet, J., A numerical investigation of
the structure of persistent shear bands in granular me-
dia, Gotechnique 41 (1991) No. 4, 599-613.
Bardet, J. P. & Proubet, J., The structure of shear bands
in idealized granular materials, Appl Mech Rev, vol
45, No. 3, part 2, March 1992. pp S118-S122.
Bardet, J. P. & Proubet, J., Numerical simulations of
shear bands in idealized granular materials, Solid
State Phenomena, Vol 23 & 24, 1992 ,pp.473-482.
Cundall, P.A., Strack, O.D.L., A discrete numerical
model for granular assemblies, Gotechnique 29
(1979) No. 1, 47-65.
Desrues J., (1987), Naissance des bandes de cisaillement
dans les milieux granulaires: exprience et thorie,
Manuel de rhologie des gomatriaux (Ed. Flix
Darve). Presses de lecole nationale des Ponts et
chausses, pp. 279-298.
Han, C. & Vardoulakis, I.G., Plane-strain compression
experiments on water-saturated fine-grained sand,
Gotechnique 41 (1991) No. 1, 49-78.
Hill, R., Acceleration waves in solids, J. Mech. Phys.
Solids 10 (1962), pp. 1-16.
Huang W., Bauer E.(2001), Numerical investigations of
shear localization in a micro-polar hypoplastic mate-
rial, submitted to Int. J. Num. Anal. Meth. Geomech.
Huang W., Nbel K. and Bauer E., Polar extension of a
hypoplastic model for granular materials with shear
localization, Mechanics of Materials, Sep 2002, Vol.
34 Issue 9, p563, 14p
ITASCA Consulting Group, Inc., Particle Flow Code in 2
Dimensions Fish in PFC. Software Manual PFC2D
Version 3.0.
ITASCA Consulting Group, Inc., Particle Flow Code in 2
Dimensions Theory and Background. Software
Manual PFC2D Version 3.0.
Iwashita, K. & Oda, M., Rolling resistance at contacts in
simulation of Shear Band development by DEM,
Journal of Engineering Mechanics, March 1998. pp.
285-292.
Mandel, J. (1964). Propagation des surfaces de disconti-
nuite dans un milieu elasto-plastique, Stress waves in
anelastic solids. Springer-Verlag, Berlin.
Mhlhaus, H., Vardoulakis, I., The thickness of shear
bands in granular materials, Gotechnique 37 (1987),
271-283.
Oda M. and Kazama H., Microstructure of shear bands
and its relation to the mechanism of dilatancy and
failure of dense granular soils, Gotechnique 48
(1998) No. 4, 465-481.
Oda, M. & Iwashita, K., Study on couple stress and shear
band development in granular media based on numer-
ical simulation analyses, International Journal of En-
gineering Science 38 (2000). 1713-1740.
Rudnicki,J. W., Rice, J.R., Conditions for the localiza-
tion of deformation in pressure sensitive dilatants ma-
terial, J. Mech. Phys. Solids 23 (1975), pp. 371-394.
Tejchman, J, Gudehus, G.: Shearing of a narrow granular
strip with polar quantities, J. Num. and Anal. Meth.
Geomech. 25(2001), pp.1-28.
Vardoulakis I. & Sulem J. (1995). Bifurcation Anlisis
in Geomechanics. Blackie Academic & Professional.