# Study of triaxial loading of segregated granular assemblies through
experiments and DEM simulations

Venakata Rama Manoj Pola  Ratna Kumar Annabattula<EMAIL_ADDRESS>
###### Abstract

A simple position-dependent body force-based confinement for simulating
triaxial tests using the Discrete Element Method is presented. The said method
is used to perform triaxial simulations on mono-disperse and segregated
assemblies of glass spheres. The macroscopic load response obtained in
simulations is validated with experimental load response. A mesh construction
algorithm is presented to check whether the confinement applied in the
triaxial simulations is accurate. The particle displacement data obtained from
triaxial simulations are used to obtain a particle-wise average strain tensor.
This is further used to compare the strain localisation between the mono-
disperse and segregated assemblies. It is observed that, in the segregated
assembly, the interface between the two particle phases acts as a barrier for
strain localisation, and the smaller particles preferentially undergo a higher
degree of shear strain on average.

###### keywords:

granular materials , segregated , discrete element method , triaxial test ,
shear strain , flexible membrane

††journal: arXiv

[inst1]organization=Department of Mechanical Engineering, Indian Institute of
Technology Madras, city=Chennai, postcode=600036, state=Tamil Nadu,
country=India [inst2]organization=Center for Soft and Biological Matter,
Indian Institute of Technology Madras,city=Chennai, postcode=600036,
state=Tamil Nadu, country=India

##  1 Introduction

Discrete Element Method(DEM) [1] is a numerical tool to simulate granular
materials. Most widely used for DEM treats individual particles as distinct
free-moving bodies whose state (positions and velocities) are updated at
regular intervals. To update the state of the particles, the forces on the
particles are calculated at the beginning of each interval. The forces on the
particles are obtained either through contact with other bodies or through a
force field. Various contact force models have been used to predict the force
transmitted through contact between two particles or between a particle and
geometry [2, 3, 4, 5]. After the force is obtained, the net acceleration
vector of the particle is calculated using Newton’s second law of motion. This
acceleration vector is temporally integrated to obtain the updated velocities
and positions of the particles at the subsequent time step. This process is
continued for the desired amount of time.

Triaxial tests have been simulated in DEM using many different approaches. The
main challenge in simulating triaxial tests is accurately modelling the
membrane and confinement pressure. The most common approach is to model the
membrane with bonded particles and construct area elements with the membrane
particles as nodes to calculate the force required for confinement [6, 7, 8,
9, 10]. Another approach is to model the membrane with servo-controlled
geometric elements to maintain the confinement [11]. A variation of this
method is to use a rigid cylinder whose radius can be controlled to maintain
confinement [12, 13]. More advanced approaches include modelling the membrane
with a finite element mesh and coupling DEM with FEA to simulate confinement
[14] and modelling the membrane with deformable discrete elements called
PFacets [15].

In this paper, the membrane is modelled with bonded particles, as has been
done before. However, confinement is applied as a position-dependent body
force on the particles comprising the membrane. There are two main advantages
of this approach. Firstly, it is relatively uncomplicated to incorporate into
DEM, as additional code is not required to track area elements. Secondly, it
is computationally less expensive for the same reason, i.e., additional
computational power is not required to track a mesh. There are some drawbacks
to this approach, too. Firstly, due to the nature of the confinement force, we
can only simulate triaxial tests with a non-rotatable platen. This is because
the inherent assumption in the approach is that the membrane particles move in
an axis-symmetric manner. Such axis-symmetric deformation is only seen when
the platens used for loading are non-rotatable. Also, large lateral
deformations induce erroneous extra confinement onto the specimen. This is
because of another approximation: the configuration of the membrane is such
that any area-element vector on the membrane will be radially outwards from
the membrane. This approximation can only be made when the membrane has not
deformed much.

Segregation is a prevalent phenomenon in granular materials [16]. The sources
of the segregation can be many, including flow, vibration, erosion, etc. In
this study, the cause of initiation or persistence of the segregation is not
examined. Only the effect of segregation on the material’s behaviour when
subjected to triaxial loads is examined. Experimental studies have shown a
correlation between segregation and strength of granular assemblies [17]. In
this paper, we focus on the pattern of deformation. We use DEM to simulate
triaxial loading of segregated specimens and obtain strain from the DEM data
[18, 19]. As will be seen, a counter-intuitive pattern emerges as the
specimens are loaded. We see that irrespective of the size of the larger
particles, the smaller particles preferentially undergo a larger degree of
lateral displacement and, as an implication, a higher degree of shear strain.

##  2 Triaxial test

Triaxial tests on granular materials are performed on cylindrical specimens
held in a thin latex membrane. The membrane is tightly attached to two
platens, one at the bottom and one at the top. The bottom platen is fixed, and
the top platen can move in the axial direction. The entire setup is placed in
a chamber filled with water and subsequently pressurised. A simplified
illustration of a triaxial setup can be seen in Figure 1.

Figure 1: Typical triaxial experimental setup.

Standard $38\text{\,}\mathrm{m}\mathrm{m}$ specimens comprised of soda-lime
glass spheres were used for the test. Initial packing density plays a critical
role in the force-displacement response and rheology of the granular assembly.
In order to maintain uniformity throughout, all the tests were carried out
with 70 grams of the glass beads to obtain an average packing density of 60%.
The sample was loaded by moving the top platen down at a rate of
$5\text{\,}\mathrm{mm}\text{/}\mathrm{min}$, which corresponds to a strain
rate of $0.003\,125\text{\,}{\mathrm{s}}^{-1}$. This strain rate was chosen to
reduce the compute time in the corresponding simulations. Shimadzu AGS-X™
series universal testing machine was used for loading and data acquisition.

The primary purpose of the experiments was to validate the macroscopic load
response of the simulations. Hence, all the testing parameters were identical
between the experiment and simulation. This includes packing fraction, which
was ensured by including the same number of particles in the experiment and
simulation. A total of three cases were experimentally tested, as outlined
below.

  * 1.

Mono-disperse assembly comprising $3\text{\,}\mathrm{m}\mathrm{m}$ particles.

  * 2.

Mono-disperse assembly comprising $5\text{\,}\mathrm{m}\mathrm{m}$ particles.

  * 3.

Segregated assembly comprising $5\text{\,}\mathrm{m}\mathrm{m}$ particles in
the bottom half and $3\text{\,}\mathrm{m}\mathrm{m}$ particles in the top
half.

The rest of the cases were not experimentally tested as the particle sizes lie
within the cases that have been tested. Hence, if the extreme particle sizes
yield a good match between the experiment and simulation, similar behaviour
can be expected for the intermediate particle sizes.

##  3 Triaxial DEM Simulations

###  3.1 Discrete element method

The discrete element method has been used to simulate the triaxial tests
described in section 2. Hertz-Mindlin contact model [2, 3, 4] with linear
rolling friction has been utilised to model contact between the glass beads.
The Hertz-Mindlin contact has four components, as shown below in Equations 1,
2, 3 and 4.

$\mathbf{\mathbf{F_{\text{n}}}}=\frac{4}{3}E^{*}\sqrt{R^{*}}\delta^{\frac{3}{2}}_{\text{n}}$
(1)
$\mathbf{F_{\text{t}}}=-8G^{*}\sqrt{R^{*}\delta_{\text{n}}}\delta_{\text{t}}$
(2) $\mathbf{F^{d}_{\text{n}}}=2*\sqrt{\frac{5}{6}}\left(\frac{\ln
e}{\left(\ln
e\right)^{2}+\pi^{2}}\right)\sqrt{2E^{*}\sqrt{R^{*}\delta_{\text{n}}}m^{*}}\,\mathbf{v^{\text{rel}}_{\text{n}}}$
(3) $\mathbf{F^{d}_{\text{t}}}=2*\sqrt{\frac{5}{6}}\left(\frac{\ln
e}{\left(\ln
e\right)^{2}+\pi^{2}}\right)\sqrt{8G^{*}\sqrt{R^{*}\delta_{\text{n}}}m^{*}}\,\mathbf{v^{\text{rel}}_{\text{t}}}$
(4)

where $\mathbf{F_{\text{n}}}$, $\mathbf{F^{d}_{\text{n}}}$, $\mathbf{F_{t}}$
and $\mathbf{F^{d}_{t}}$ corresponds to force due to normal stiffness, normal
damping, tangential stiffness and tangential damping, respectively.
$E^{*}$,$R^{*}$,$G^{*}$ and $m^{*}$ represent effective Young’s modulus,
effective radius of contact, effective shear modulus and effective mass
respectively and $\delta_{\text{n}}$,$\delta_{\text{t}}$,$v^{rel}_{\text{n}}$
and $v^{rel}_{\text{n}}$ represent the normal overlap, tangential overlap,
normal relative velocity and tangential relative velocity respectively. The
tangential stiffness has a ceiling value dictated by the static friction
coefficient in accordance with the Coulomb law of friction. All the parameters
and material parameters defined in the simulations can be seen in Table 1.

The rolling friction model used in the simulations can be seen in Equation 5,
where $\mathbf{\tau}$ is the torque applied on the particle to resist rolling.
$\mathbf{\tau}$ is applied in the direction opposite to the angular velocity
vector of the particle in consideration. $R_{i}$ is the distance between the
centre of the particle and the point of application of the torque. This point
is the midpoint of the line segment joining the centres of the particles in
contact. $\omega_{i}$ is the unit vector along the direction of angular
velocity.

$\mathbf{\tau}=-\mu_{\text{r}}\mathbf{F_{\text{n}}}R_{i}\omega_{i}$ (5)

The glass particles are first consolidated in a cylindrical container of
radius $38\text{\,}\mathrm{mm}$ and height $80\text{\,}\mathrm{mm}$.
Subsequently, membrane particles (section 3.2) are generated and fixed between
the top and bottom cylindrical platens, similar to the experiment. A
confinement force is then gradually applied (section 3.3) by applying a
central force on the membrane particles. Upon confinement, the top platen is
then moved down at a rate of
$5\text{\,}\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{m}\mathrm{i}\mathrm{n}$ up to
a total displacement of $10\text{\,}\mathrm{m}\mathrm{m}$. A total of five
specimens are simulated. Two mono-disperse specimens comprising of
$3\text{\,}\mathrm{m}\mathrm{m}$ and $5\text{\,}\mathrm{m}\mathrm{m}$ diameter
particles, respectively and three segregated specimens comprising of
$3\text{\,}\mathrm{m}\mathrm{m}$ diameter particles in the bottom half and
$4\text{\,}\mathrm{m}\mathrm{m}$, $4.5\text{\,}\mathrm{m}\mathrm{m}$ and
$5\text{\,}\mathrm{m}\mathrm{m}$ diameter particles in the top half
respectively. The initial and final configuration of the one mono-disperse
specimen and one segregated specimen can be seen in Figure 2. Altair EDEM™has
been used to run the DEM simulations presented in this article. Custom
particle-generator and custom body force models have been incorporated into
the program to simulate the membrane and confinement.

Table 1: Simulation parameters Property | Value | Property | Value
---|---|---|---
Young’s modulus |  | Static coefficient of friction |
glass particles | $70\text{\,}\mathrm{GPa}$ | glass-glass | 0.5
membrane particles | $0.25\text{\,}\mathrm{GPa}$ | glass-membrane | 0.1
loading plate | $0.25\text{\,}\mathrm{GPa}$ | glass-loading plate | 0.1
Poisson’s ratio |  | Coefficient of Restitution |
glass particles | 0.22 | glass-glass | 0.7
membrane particles | 0.25 | glass-membrane | 0.5
loading plate | 0.25 | glass-loading plate | 0.5
Density |  | Coefficient of rolling friction |
glass particles | $2500\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$ | glass-glass | 0.01
membrane particles | $1100\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$ | glass-membrane | 0.01
|  | glass-loading plate | 0.01

(a)

(b)

(c)

(d)

Figure 2: (a) Initial and (b) final configuration of the mono-disperse
specimen comprising of $5\text{\,}\mathrm{m}\mathrm{m}$ particles and (c)
Initial and (d) final configuration of the segregated specimen comprising of
$5\text{\,}\mathrm{m}\mathrm{m}$ in the bottom half and
$3\text{\,}\mathrm{m}\mathrm{m}$ particles in the top half.

###  3.2 Simulating the membrane

The membrane used in triaxial tests is simulated by a cylindrical layer of
bonded membrane particles around the glass particles. The radius of the
membrane particles is chosen as $0.5\text{\,}\mathrm{m}\mathrm{m}$, which is a
third of the smallest particle in the specimen. The properties of the membrane
particles can be seen in Table 1. The membrane particles are arranged in a
hexagonal close packing structure to minimise inter-particle gaps, as seen in
Figure 3.

Figure 3: Visualisation of the hexagonal close packing arrangement of membrane
particles in a cylinder that will blanket the glass beads.

Immediately after the generation of membrane particles, bonds are generated
between particles in contact according to a bond formation radius, which was
set at 150% the physical radius of the particle. This allows for the bonding
of the immediate neighbours amongst the particles.The bond normal stiffness is
set to $300\text{\,}\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ and a tangential
stiffness of $150\text{\,}\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$. These values
have been obtained by performing uni-axial tensile loading of a
$$8\text{\,}\mathrm{m}\mathrm{m}$\times$4\text{\,}\mathrm{m}\mathrm{m}$$ strip
of membrane. This is followed by replicating the tensile test in DEM with our
particle membrane and calibrating the parameters to match the macroscopic
stiffness of the strip. The comparison between the experiment and simulation
can be seen in Figure 4.

Figure 4: Comparison between experiment and simulation for uni-axial loading
of a $$8\text{\,}\mathrm{m}\mathrm{m}$\times$4\text{\,}\mathrm{m}\mathrm{m}$$
strip of membrane.

###  3.3 Simulating confinement

Confinement has been modelled by applying a radially varying central force on
the membrane particles. This implementation of confinement is only suitable
for axisymmetric modes of membrane deformation. Such axisymmetric deformation
is seen only when the top platen is non-rotatable; hence, only those cases
have been simulated. This confinement force has been derived, assuming each
membrane particle represents a small area of the membrane. These area elements
can be assumed to be hexagonal, but the shape is irrelevant to the confinement
force. The total confinement force ($F_{c}$) due to the desired confinement
pressure ($P_{c}$) is calculated (Equation 6), from which the initial force
($F^{\text{init}}_{p}$) on each particle is calculated (Equation 7). Assuming
the membrane expands axisymmetrically, the area represented by each element
and, subsequently, the force on each particle will increase linearly with an
increase in the radial position of the particle from the central axis. Also,
the reduction in the height of the membrane is incorporated by proportionally
reducing the height represented by each particle proportionally. This results
in a proportional reduction of area and, hence, confinement force. The
position-dependent force at a time t ($F_{\text{p}}^{\text{t}}$) on each
particle is therefore given by Equation 8 below.

$F_{\text{c}}=P_{\text{c}}\times 2\pi
r_{\text{s}}^{\text{init}}h_{\text{s}}^{\text{init}}$ (6)
$F^{\text{init}}_{\text{p}}=\frac{F_{\text{c}}}{N_{\text{p}}}$ (7)
$F_{\text{p}}^{\text{t}}=\frac{h_{\text{s}}^{\text{init}}-\Delta\text{h}_{\text{t}}}{h_{\text{s}}^{\text{init}}}\times\frac{r_{\text{t}}}{r_{\text{s}}^{\text{init}}}\times
F^{\text{init}}_{\text{p}}$ (8)

where $r_{\text{s}}^{\text{init}}$ and $r_{\text{t}}$ are the initial and
current radial distance of the membrane particle from the axis of the membrane
cylinder, respectively. $h_{\text{s}}^{\text{init}}$ and
$\Delta\text{h}_{\text{t}}$ are the initial height and the change in the
height of the specimen at time $t$, respectively. $N_{\text{p}}$ corresponds
to the total number of membrane particles.

It should be noted that sudden application of this force will cause
instabilities in the system. Hence, the force is gradually applied over a time
of $1\text{\,}\mathrm{s}$. The rationale behind this model of force
application is as follows. It can be assumed that a membrane particle
deforming from the initial value of $r_{\text{s}}^{\text{init}}$ becomes part
of a new imaginary membrane of radius equal to the final value of
$r_{\text{t}}$. This can be seen in Figure 5, as the membrane particle
highlighted in blue moves laterally, it is assumed to be part of an imaginary
membrane with a radius $r_{\text{t}}$. The surface area of the imaginary
membrane would be scaled by ${r_{\text{t}}}/{r_{\text{s}}^{\text{init}}}$ when
compared to the original membrane. The force due to the confinement pressure
would be directly proportional to the surface area; hence it is scaled by
${r_{\text{t}}}/{r_{\text{s}}^{\text{init}}}$ as seen in Equation 8.

(a)

(b)

Figure 5: (a)Initial and (b) Final configuration of the membrane particles
indicating the initial radius $r$ and final radius $R_{c}$ from the central
axis for one of the particles.

This approach to modelling confinement has certain shortcomings. Firstly, the
simulations are limited to axisymmetric cases, which means cases with rotating
loading platens cannot be simulated. Apart from that, the assumptions made
about the deformation method of the area element represented by the membrane
particles are approximate and assume ideal deformation cases. The following
subsection deals with validating the confinement force applied to establish
its accuracy.

###  3.4 Validating confinement

The confinement model presented above does not involve any feedback mechanism
or area calculation. This means that the confinement pressure is not
controlled during the simulation. To ensure the reliability of the confinement
model, a retrospective post-processing procedure to measure the confinement
pressure is proposed. The first step in computing the confinement pressure is
obtaining the required simulation data. We need the contact forces transmitted
between the membrane particles and the specimen particles and the positions of
points of contact to reconstruct the membrane as a mesh. The next step is to
arrange rectangles along the membrane. The position of the rectangles is
represented in polar coordinates. The $z$ and $\theta$ components of the
position of the rectangles are fixed to span the whole membrane. The elements
are angled such that their surface normals are always radially outward. To
determine the radial component of the position of a rectangle, the radial
positions of all the contact points whose $z$ and $\theta$ components fall
within the limits of the rectangle are aggregated. This ensures the array of
rectangles tile the entire membrane. Subsequently, the centres of each of the
rectangles are connected to form a mesh composed of quadrilateral elements.
For a specific element in this mesh, the four nodes’ radial and $z$ components
are used as limits to find the contact points that lie ‘within’ the element.
The contact force vectors associated with each contact point are added to find
the total force transmitted through the element. This cumulative force vector
and the area vector of the element are used to calculate the pressure
associated with the element.

An illustrative example of the methodology explained above in 2D can be seen
in Figure 6. Figure 6(a) shows the membrane particles and the force the
particle is applying on the specimen. Tiling the membrane with rectangles is
seen in Figure 6(b), except the rectangles are replaced by vertical lines. The
centres of the lines are joined to generate a mesh, as seen in Figure 6(c),
wherein the area vectors, along with the cumulative force vector associated
with each mesh element, are also known. This pressure is averaged over all the
elements of the mesh to obtain the average confinement pressure being applied.

(a)

(b)

(c)

Figure 6: (a) The membrane particles and contact force vectors between the
membrane particle and specimen particles along with the contact points (b)
Placing the rectangles (lines in 2D) along the membrane. (c) Creating the
final mesh from the centres of the rectangles (lines) shown in (b), along with
the corresponding normal vectors in black and contact force vectors in red.

The average pressure in the initial, middle and final configurations of the
simulations can be seen in Figure 7. Both the mono-disperse cases, along with
the segregated case comprising $3\text{\,}\mathrm{m}\mathrm{m}$ and
$5\text{\,}\mathrm{m}\mathrm{m}$ particles, are tested. Other segregated cases
($3\text{\,}\mathrm{m}\mathrm{m}$ with $4\text{\,}\mathrm{m}\mathrm{m}$ and
$3\text{\,}\mathrm{m}\mathrm{m}$ with $4.5\text{\,}\mathrm{m}\mathrm{m}$) are
not tested as any discrepancy arising in the intermediary particle sizes would
also be observed in the extreme particle sizes that have been tested. The
desired confinement pressure is $100\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}$.
It can be seen that the confinement pressure has a very small positive
gradient, but within the $5\text{\,}\mathrm{m}\mathrm{m}$ stroke applied,
which corresponds to a macroscopic strain of 0.0625, the difference is less
than 5%. This difference is sufficiently small enough to ensure the accuracy
of simulations.

Figure 7: Evolution of average confinement pressure for the mono-disperse
assemblies and the segregated assembly comprising the
$3\text{\,}\mathrm{m}\mathrm{m}$ and $5\text{\,}\mathrm{m}\mathrm{m}$
particles.

##  4 Results and discussions

###  4.1 Load response

We first look at the load response of the simulations and, wherever available,
the corresponding experiments. This mainly aims to validate our simulation
with the corresponding experiment. A good match between the simulation and
experiment means the contact parameters that we have chosen are closer to the
realistic values. Figure 8(a) shows the response for the
$3\text{\,}\mathrm{m}\mathrm{m}$ mono-disperse assembly and the
$5\text{\,}\mathrm{m}\mathrm{m}$ mono-disperse assembly. We see a good
agreement between the experiments and the simulation for both cases. We see
that the $3\text{\,}\mathrm{m}\mathrm{m}$ assembly exhibits a higher load-
bearing capacity than the $5\text{\,}\mathrm{m}\mathrm{m}$ assembly, seen both
in experiment and simulation. Figure 8(b) shows the comparison of the load
response obtained for the segregated assemblies. This includes the
experimental result for the segregated assembly comprising the
$3\text{\,}\mathrm{m}\mathrm{m}$ and $5\text{\,}\mathrm{m}\mathrm{m}$
particles, which again is in good agreement with the simulation. All the
segregated assemblies have similar load responses. When compared to the mono-
dispersed assemblies, it is seen that the segregated assemblies have a load
response closer to the $5\text{\,}\mathrm{m}\mathrm{m}$ particles.

(a)

(b)

Figure 8: (a)The deviatoric load response for the two mono-disperse assemblies
comprising experiment and simulation. (b)The deviatoric load response for the
three segregated assemblies, comprising experimental response of the assembly
comprising $3\text{\,}\mathrm{m}\mathrm{m}$ and
$5\text{\,}\mathrm{m}\mathrm{m}$ particles.

The main aim of this comparison is to validate the load response of the
simulations and experiments. Good agreement between the simulation and
experiment is a decent indicator towards accurate inter-particular contact
behaviour. Especially since it has already been established that the correct
confinement is being applied, this, in turn, translates to accurate rheology
during loading. It should also be noted that we have experimentally verified
cases with extreme particle sizes ($3\text{\,}\mathrm{m}\mathrm{m}$ and
$5\text{\,}\mathrm{m}\mathrm{m}$). If the discrepancy in applying the
confinement pressure is not present for extreme particle sizes in a segregated
assembly, such a discrepancy would not be present for all the intermediate
particle sizes.

###  4.2 Strain

The particle position data is used to calculate a ‘per-particle average’
strain data, and the methodology is described in Figure 9. First, Delaunay
triangulation is used to mesh the domain using particle centres as nodes. This
results in a mesh comprising tetrahedral elements. Some of the elements are
deleted if they fail based on a maximum aspect ratio criterion to avoid
erroneous strain calculations. The value of the maximum allowable aspect ratio
was set to 5. Subsequently, mesh deformation is tracked via nodal displacement
obtained from DEM data. From this, constant strain tetrahedron formulation is
used to calculate element-wise strain using the Green strain formulation,
which can be seen in Equation 9. Green strain is used because the strain is
calculated with respect to the initial/original configuration. Hence, the
small-strain approximation cannot be made. The per-particle average strain
tensor is obtained by averaging the components of the strain tensor, which
belong to all the elements for which the particle is a node. More details
about strain calculation used here can be seen in this work [19].

$E_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}+u_{k,i}u_{k,j}\right)$ (9)

Figure 9: A simplified workflow for calculating an arbitrary component of the
strain tensor for a four-particle system

For the purpose of visual representation, the strain data is interpolated over
a mesh constructed by Delaunay triangulation of the particle centres. This
triangulation is only for visual representation and is separate, albeit
similar to the mesh used for strain calculation. A $y$-normal plane passing
through the centre of the assembly is used to slice the interpolated data, and
the $zx$-component of the strain is shown on this slice.

Figure 10 shows the $zx$ component of Green strain tensor distribution for the
$3\text{\,}\mathrm{m}\mathrm{m}$ and $5\text{\,}\mathrm{m}\mathrm{m}$ mono-
disperse assemblies respectively. Stain localization is seen distributed
throughout the assembly. It is difficult to demarcate the shear bands due to
the small value of the applied stroke. Figure 11 shows the same strain
distribution seen above this time for the segregated assembly. Here, an
interesting pattern emerges. Diffused strain distribution is seen in the
region comprising $3\text{\,}\mathrm{m}\mathrm{m}$ particles. Comparatively,
the regions comprising the $4\text{\,}\mathrm{m}\mathrm{m}$,
$4.5\text{\,}\mathrm{m}\mathrm{m}$ and $5\text{\,}\mathrm{m}\mathrm{m}$
particles are almost devoid of this strain component. This phenomenon is
unexpected as both layers are bearing the same load, leading to the
expectation that the assembly with lower load-bearing capacity would undergo
more shear strain. By that train of thought, at least for the assembly
comprising the $3\text{\,}\mathrm{m}\mathrm{m}$ and
$5\text{\,}\mathrm{m}\mathrm{m}$, we would expect the
$5\text{\,}\mathrm{m}\mathrm{m}$ layer to undergo more shear strain. This is
because, in the mono-disperse case, it has been observed that the
$5\text{\,}\mathrm{m}\mathrm{m}$ particle assembly was shown to have a lower
load-bearing capacity. One more observation, which implies the preceding
observation, is that the interface between the two layers of particles seems
to act as a barrier for strain localization. While this effect has almost been
taken for granted till now, we don’t have an explanation for why this happens
instead of strain being distributed relatively uniformly throughout the
assembly across the interface.

(a) $3\text{\,}\mathrm{m}\mathrm{m}$ monodisperse assembly

(b) $5\text{\,}\mathrm{m}\mathrm{m}$ monodisperse assembly

Figure 10: For the mono-disperse assemblies $zx$ component of the Green strain
tensor along a plane passing through the centre of the specimen and
perpendicular to the $y$-axis.

(a) $4\text{\,}\mathrm{m}\mathrm{m}$ particles in bottom half and
$3\text{\,}\mathrm{m}\mathrm{m}$ particles in top half

(b) $4.5\text{\,}\mathrm{m}\mathrm{m}$ particles in bottom half and
$3\text{\,}\mathrm{m}\mathrm{m}$ particles in top half

(c) $5\text{\,}\mathrm{m}\mathrm{m}$ particles in bottom half and
$3\text{\,}\mathrm{m}\mathrm{m}$ particles in top half

Figure 11: For the segregated assemblies $zx$ component of the Green strain
tensor along a plane passing through the centre of the specimen and
perpendicular to the y-axis.

To quantitatively compare this phenomenon, a histogram comparison of the
distribution of $zx$ component of the Green strain tensor can be seen in
Figure 12. The left y-axis corresponds to the number of small
($3\text{\,}\mathrm{m}\mathrm{m}$) particles, and the right y-axis corresponds
to the large particles. Both the y-axes have been scaled proportional to the
total number of particles of the particular kind. This ensures normalization
with respect to the number of particles. Figure 12(a) compares the strain
amongst the mono-disperse assemblies. We see that between the strain values of
0 and 0.03, the $3\text{\,}\mathrm{m}\mathrm{m}$ particles are marginally more
in proportion, and the inverse is true post the strain value of 0.3. That
being said, the strain distribution is fairly similar in the two cases.
However, the same cannot be said for the segregated assemblies. Figures 12(b),
12(c) and 12(d) compares the strain distribution between the
$3\text{\,}\mathrm{m}\mathrm{m}$ and the $4\text{\,}\mathrm{m}\mathrm{m}$,
$4.5\text{\,}\mathrm{m}\mathrm{m}$ and $5\text{\,}\mathrm{m}\mathrm{m}$
particles in the segregated assembly respectively. In all these cases, a
higher number of larger particles undergo strain less than 0.02. Above that
value, the $3\text{\,}\mathrm{m}\mathrm{m}$ particles occupy a significantly
larger proportion. We reiterate that the comparison is not between the actual
number of particles but the proportion of particles. This comparison is
possible in the graphs because the y-axes have been scaled in proportion to
the total number of particles of the particular kind. This re-confirms our
previous observation about the $3\text{\,}\mathrm{m}\mathrm{m}$ particles
preferentially undergoing larger shear strain compared to the larger
particles. The average shear strain of the large and small particles for the
various assemblies can be seen in Table 2. The difference in the average shear
strain between the small and large particles is also shown. For the mono-sized
assemblies ($3\text{\,}\mathrm{m}\mathrm{m}$ and
$5\text{\,}\mathrm{m}\mathrm{m}$), the larger particles have a marginally
higher average shear strain compared to the smaller particles by an amount of
0.6. Surprisingly, for the segregated assembly comprising the
$3\text{\,}\mathrm{m}\mathrm{m}$ and $5\text{\,}\mathrm{m}\mathrm{m}$ also,
the average shear strain difference is 0.6, except the
$3\text{\,}\mathrm{m}\mathrm{m}$ particles have the higher shear strain.
Amongst the segregated assemblies, the average shear stain difference seems
inversely related to the size of the larger particles.

Some visual observations can be drawn from the experimental triaxial test of
the segregated assembly comprising $5\text{\,}\mathrm{m}\mathrm{m}$ and
$3\text{\,}\mathrm{m}\mathrm{m}$. As the specimen is loaded, we see bulging
only in the region containing the $3\text{\,}\mathrm{m}\mathrm{m}$ particles
as seen in Figure 13. It is seen that $5\text{\,}\mathrm{m}\mathrm{m}$
particles undergo almost no lateral deformation.

(a) Mono sized assemblies

(b) $3\text{\,}\mathrm{m}\mathrm{m}$/$4\text{\,}\mathrm{m}\mathrm{m}$
segregated assembly

(c) $3\text{\,}\mathrm{m}\mathrm{m}$/$4.5\text{\,}\mathrm{m}\mathrm{m}$
segregated assembly

(d) $3\text{\,}\mathrm{m}\mathrm{m}$/$5\text{\,}\mathrm{m}\mathrm{m}$
segregated assembly

Figure 12: Histograms comparing the distribution of the $zx$ component of Green strain tensor. Table 2: Comparison of the average value of the $zx$ component of the Green strain tensor ($\bar{E}_{zx}$) between the large and small particles in all the samples that have been simulated. For the mono-sized samples, the average value corresponds to the average in the respective sample. The errors shown are the standard error values calculated from the standard deviation of the strain distribution. Assembly type |  $\bar{E}_{zx}$ small |  $\bar{E}_{zx}$ large | Difference(small-large)
---|---|---|---
Mono-sized(5mm/3mm) | 0.026 $\pm$ 0.0005 | 0.032 $\pm$ 0.0015 | -0.006 $\pm$ 0.0016
Segregated(4mm/3mm) | 0.035 $\pm$ 0.0008 | 0.022 $\pm$ 0.0014 | 0.013 $\pm$ 0.0016
Segregated(4.5mm/3mm) | 0.035 $\pm$ 0.0008 | 0.026 $\pm$ 0.0015 | 0.009 $\pm$ 0.0017
Segregated(5mm/3mm) | 0.033 $\pm$ 0.0007 | 0.027 $\pm$ 0.0015 | 0.006 $\pm$ 0.0017

(a)

(b)

(c)

(d)

Figure 13: Snapshots of the experimental triaxial test of the segregated
assembly at (a) initial state (b) $5\text{\,}\mathrm{m}\mathrm{m}$ stroke (c)
$10\text{\,}\mathrm{m}\mathrm{m}$ stroke (d) $15\text{\,}\mathrm{m}\mathrm{m}$
stroke. It can be observed that the smaller particles, occupying the upper
half, bulge out, while the larger particles do not experience such bulging.

##  5 Conclusion

The effect of segregation on the strain localization in granular assemblies
subjected to triaxial loading has been studied. The discrete element method
has been employed to simulate the triaxial loading of granular assemblies.
Hexagonally arranged bonded particles were used to model the membrane
enveloping the sample. A position-dependent force on the membrane particles
was used to apply confinement to the assembly. A mesh reconstruction method
has been employed to verify if the correct confinement has been applied to the
assembly. Some of the simulations have been verified with experiments by
comparing the macroscopic load response, which was found to be in good
agreement. The average particle Green strain was calculated using Delaunay
triangulation, from which the $zx$ component of the strain was shown along the
y-normal plane passing through the centre of the assembly. Finally, the
following conclusions have been drawn.

  1. 1.

The position-dependent body-force-based confinement can be used to accurately
apply confinement, especially in the initial phase of loading.

  2. 2.

When a segregated sample is subject to triaxial loading, it is observed that
only the smaller particles undergo lateral deflection, which is translated
into a higher amount of shear strain in the region being occupied by the
smaller particles. This effect was also visually observed in the
experimentally tested segregated sample.

  3. 3.

The difference between the average shear strain experienced by the particles
in the segregated assemblies decreases as the size of the larger particles
increases, mainly due to the increase in shear strain in the larger particles.

## Conflict of interest statement

The authors declare that there is no conflict of interest.

## Acknowledgements

The authors are grateful to Altair Engineering India for providing funding,
technical support and the software required to run the simulations shown in
this article.

## References

  * [1] Peter A Cundall and Otto DL Strack.  A discrete numerical model for granular assemblies.  Géotechnique, 29(1):47–65, 1979.
  * [2] Heinrich Hertz.  Ueber die berührung fester elastischer körper.  Journal für die reine und angewandte Mathematik, 1882(92):156–171, 1882.
  * [3] R. D. Mindlin.  Compliance of elastic bodies in contact.  Journal of Applied Mechanics, 16(3):259–268, 1949.
  * [4] R. D. Mindlin and H. Deresiewicz.  Elastic spheres in contact under varying oblique forces.  Journal of Applied Mechanics, 20(3):327–344, 1953.
  * [5] Hide Sakaguchi, Eiji Ozaki, and Tohru Igarashi.  Plugging of the flow of granular materials during the discharge from a silo.  International Journal of Modern Physics B, 07(09n10):1949–1963, 1993.
  * [6] Yan Qin, Chun Liu, Xiaoyu Zhang, Xingang Wang, Bin Shi, Yue Wang, and Shang Deng.  A three-dimensional discrete element model of triaxial tests based on a new flexible membrane boundary.  Scientific Reports, 11(1):4753, 2021.
  * [7] S. Fazekas, J. Török, J. Kertész, and D. E. Wolf.  Morphologies of three-dimensional shear bands in granular media.  Physical Review E, 74(3):031303, 2006.
  * [8] J. de Bono, G. Mcdowell, and D. Wanatowski.  Discrete element modelling of a flexible membrane for triaxial testing of granular material at high pressures.  Géotechnique Letters, 2(4):199–203, 2012.
  * [9] Kai Wu, Weichen Sun, Songyu Liu, and Xiang Zhang.  Study of shear behavior of granular materials by 3D DEM simulation of the triaxial test in the membrane boundary condition.  Advanced Powder Technology, 32(4):1145–1156, 2021.
  * [10] Y.-H. Wang and S.-C. Leung.  A particulate-scale investigation of cemented sand behavior.  Canadian Geotechnical Journal, 45(1):29–44, 2008.
  * [11] Seung Jae Lee, Youssef M.A. Hashash, and Erfan G. Nezami.  Simulation of triaxial compression tests with polyhedral discrete elements.  Computers and Geotechnics, 43:92–100, 2012.
  * [12] Kai Wu, Sébastien Rémond, NorEdine Abriak, Patrick Pizette, Frédéric Becquart, and Songyu Liu.  Study of the shear behavior of binary granular materials by DEM simulations and experimental triaxial tests.  Advanced Powder Technology, 28(9):2198–2210, 2017.
  * [13] Kai Wu, Patrick Pizette, Frederic Becquart, Sébastien Rémond, NorEdine Abriak, Weiya Xu, and Songyu Liu.  Experimental and numerical study of cylindrical triaxial test on mono-sized glass beads under quasi-static loading condition.  Advanced Powder Technology, 28(1):155–166, 2017.
  * [14] Gang Ma, Wei Zhou, Xiao-Lin Chang, and Wei Yuan.  Combined FEM/DEM modeling of triaxial compression tests for rockfills with polyhedral particles.  International Journal of Geomechanics, 14(4):04014014, 2014.
  * [15] Pu Yang, Edward Kavazanjian, and Narayanan Neithalath.  Particle-scale mechanisms in undrained triaxial compression of biocemented sands: Insights from 3D DEM simulations with flexible boundary.  International Journal of Geomechanics, 19(4):04019009, 2019.
  * [16] Julio M Ottino and Devang V Khakhar.  Mixing and segregation of granular materials.  Annual Review of Fluid Mechanics, 32(1):55–91, 2000.
  * [17] Fang-wei Yu, Xiong-zhi Peng, and Li-jun Su.  Influence of grain segregation on the behavior of sand in triaxial tests.  Journal of Mountain Science, 18(10):2776–2790, 2021.
  * [18] Patricia Ann Thomas.  Discontinuous deformation analysis of particulate media.  University of California, Berkeley, 1997.
  * [19] C. O’Sullivan.  Particulate Discrete Element Modelling: A Geomechanics Perspective.  CRC Press, 2011.