Ljubljana 2010

Simulation of Granular
Materials and Coupling with
Solids and Fluids

Peter Wriggers, Christian Wellmann and Bircan Avci

 Simulation of Granular Material 2
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Motivation

• complex material behavior

• phenomenological approach:
– assume continuum
– elasto-plastic constitutive equation
• drawbacks:
– numerous parameters
– stability, accuracy and
technical issues for
large deformations
• multi-scale model:
particle model for critical regions,
elsewhere continuum
• homogenization to check particle
model and adapt continuum
model

pictures from [White & Bolton,

 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

motivation: DEM

• continuum approach issues

– appropriate constitutive equation not
obvious
– treatment of discontinuities complex
– lack of scale separation in hot-spot
regions might invalidate continuum
approach
• discrete approach
– modeling of the microstructure
– microscopic constitutive equation:
inter-particle contact law
– arbitrary deformations
– computationally expensive
– correct bc‘s need DEM-FEM coupling
 Simulation of Granular Material 4
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Outlook

• discrete element model (DEM)

– geometry, contact

• cubical samples for homogenization

– periodic boundaries, load application

• reference material: Leighton Buzzard Sand

– adapt DEM parameters

• generation of random samples

• triaxial test series and results

 Simulation of Granular Material 5
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

DEM, Geometry

• superellipsoid particle geometry [Barr, 1981]

0.5 1 1.5

0.5

1
• smooth & convex for 0 < εi < 2
• arbitrary aspect ratios & roundness / angularity

1.5
 Simulation of Granular Material 6
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Contact Formulation

• admit small overlap

• Hertz contact theory:
– consider overlap as elastic
deformation at contact points
– requires principal curvatures

• split contact force

 Simulation of Granular Material 7
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Tangential Part

• integration over time

• elastic compliances [Mindlin,1949]

• compute trial force & check Coulomb crit.

 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

contact detection

• global search via Verlet-list [Verlet,1967]

• local search for contact points as 2D
unconstrained optimization problem
• parameterize contact direction

• determine points where surface normal is (anti-)parallel to

• contact points correspond to

• combined Newton and

Levenberg-Marquardt
method

C. Wellmann, C. Lillie and P. Wriggers (2008)

 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

spheres superellipsoids
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

parallelization with OpenMP

• resultant forces & torques are accumulated in a

loop over all neighbor-lists
• standard parallelization leads to race condition
• sort particles according to spatial position

CPU 1 CPU 2

chunk chunk
1 2

chunk chunk 2
1
Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

performance of parallel method

spheres

superellipsoids
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

performance of parallel method

spheres

superellipsoids
 Simulation of Granular Material 13
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Cuboidal Unit Cell

• periodic boundaries: no spurious effects

from particle-wall contacts
• average stress from boundary contact forces

• arbitrary triaxial loading by controlling Li’s:

– strain-controlled:

– stress-controlled by considering boundaries

as walls of thickness t
 Simulation of Granular Material 14
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Leighton Buzzard Sand Fraction B

• size: 0.6 – 1.18 mm

• silica sand
• high resistance against crushing
• widely used in research
• DEM material parameters:
– E=50 GPa, ν=0.2
[Kingston et al., 2008]
– friction [Rowe, 1962], [Ishibashi, 1994]: μ=0.24

reference density, g/ Фmin Фmax d10, d50, d60, U=d60/ shape

cm^3 mm mm mm d10
Rowe, 1962 0.55 0.67
Stroud, 1971 0.56 0.67 0.65 0.85
Palmeira, 1987 2.66 0.56 0.67 0.8 1.3 angular
Schnaidt, 1990 2.65 0.56 0.68 0.8 1.3 subangular-
subrounded
Lings & Dietz, 2.65 0.56 0.66 0.64 0.78 0.81 1.27 rounded
2004
White & Bolton, 2.65 0.56 0.66 0.84 rounded
2004
Kingston,... 2008 2.65 0.57 0.66 rounded-subrounded
 Simulation of Granular Material 15
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Particle Geometry
“sphericity”
• particle dimensions
[Clayton et al., 2009]:

• superellipsoid radii
with normal distribution
[Clayton et al.,
2009]

rounded subangular very angular

“roundness”
• no agreed shape
measure, classification by
comparison charts:
rounded-subrounded
• superellipsoid exponents

well-rounded subrounded angular

 Simulation of Granular Material 16
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Generation of Random Samples

• random addition of particles to unit cell, Ф ~ 25%

• algorithm to tighten package:

while Ф < Ф*
loop over particles in random order
determine distances to nearest neighbors
choose by random either
translate & scale particle
or
rotate particle
end
end
end
• compress package by application
of low pressure (1 kPa)
Simulation of Granular Material 17
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Leighton
Buzzard:
Фmin = 0.55
Фmax = 0.68
Simulation of Granular Material 17
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Leighton
Buzzard:
Фmin = 0.55
Фmax = 0.68
Simulation of Granular Material 17
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Leighton
Buzzard:
Фmin = 0.55
Фmax = 0.68
 Simulation of Granular Material 18
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Triaxial Test Series

• triaxial tests for 4 confining stresses

• 4 cell sizes

L 10 15 20 30

ns 10 8 5 4

np 1,7 k 6,0 k 14,0 k 45,0 k

• initial solid fraction Ф ~ 0.67
• corresponds to
relative density of 37%
=> loose to medium state
 Simulation of Granular Material 19
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Loading Velocity

• strain controlled compression

• quasi-static deformation modeled via

contact and global damping
• how slow for quasi-static?
• [Cruz et al., 2005]
• quasi-static: I < 10-4 – 10-5
[Agnolin & Roux, 2007]
• tests for I from 10-5 - 10-3
• choose I = 10-4
 Simulation of Granular Material 20
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Results for L = 10

• stress ratio & volumetric strain

independent of confining stress
• high initial modulus
• constant stress ratio at ~ 10%
axial strain
• initial contraction followed by
constant dilation though low
initial relative density of 37%
• standard deviation of the order
of 5% of mean
Simulation of Granular Material 21
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010
 Simulation of Granular Material 22
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Simulation vs. Experiment

• experiments from [Schnaid, 1990]:

– triaxial tests on loose, medium
and dense Leighton Buzzard
14/25 sand
– strain results:
– good agreement
• raise of inter-particle friction
coefficient yields no significant
increase of bulk shear strength
• variation in shear strength might
be due to non-convexity of real
grains
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

coupling with FEM: surface discretization

• coupling via contact between particles

and FE-surface mesh
• original FE-surface hard to handle
• discretize FE-surface by triangular mesh
• contact formulation analogous to
inter-particle contact
• static equivalent force at fe nodes
• not applicable if
– edge intersects particle
– node lies inside particle
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

integration of coupled system

 Simulation of Granular Material 25
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Numerical Model

• glass beads modeled by DEM with Hertzian contact forces

• sample generation
1. random packing with initial solid fraction of 0.5 and uniform size
distribution
2. shrink rigid cylinder (without friction and gravity) until p=1 kPa
3. replace rigid cylinder by rubber membrane modeled by FE
and rigid top and bottom plates
4. apply lateral confining pressure
5. move top and bottom plates at constant velocity
• rubber membrane discretized by 16 x 12 4-node
membrane elements [Gruttmann & Taylor, 1992]
– large elastic deformations
– incompressibility
– using three-term Ogden material
 Simulation of Granular Material 26
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Recorded Data

• resultant force on top and bottom plate

=> stress in axial direction
• volume inside membrane =>
volumetric strain
• cylinder height => axial strain
• interpretation in terms of
– principal stress ratio (axial / lateral)
– volumetric strain
– solid fraction
 Simulation of Granular Material 26
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Recorded Data

• resultant force on top and bottom plate

=> stress in axial direction
• volume inside membrane =>
volumetric strain
• cylinder height => axial strain
• interpretation in terms of
– principal stress ratio (axial / lateral)
– volumetric strain
– solid fraction
 Simulation of Granular Material 27
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Numerical Model for Sperical Glass Beads

• glass beads modeled by DEM with Hertzian contact forces

• sample generation
1. random packing with initial solid fraction of 0.5 and uniform size
distribution
2. shrink rigid cylinder (without friction and gravity) until p=1 kPa
3. replace rigid cylinder by rubber membrane modeled by FE
and rigid top and bottom plates
4. apply lateral confining pressure
5. move top and bottom plates at constant velocity
• rubber membrane discretized by 16 x 12 4-node
membrane elements [Gruttmann & Taylor, 1992]
– large elastic deformations
– incompressibility
– using three-term Ogden material
 Simulation of Granular Material 28
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Solid Fraction M, p=25

kPa
• local solid fraction from spherical RVEs
• RVE centers on regular grid
• RVE size: approx. 100 grains per RVE
• exact analytical calculation of solid
fraction for each RVE
• densification at loading plates
 Simulation of Granular Material 28
Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Solid Fraction M, p=25

kPa
• local solid fraction from spherical RVEs
• RVE centers on regular grid
• RVE size: approx. 100 grains per RVE
• exact analytical calculation of solid
fraction for each RVE
• densification at loading plates
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

outlook

• concurrent multiscale model

of granular materials
FEM
• use overlapping domain
like in Arlequin method
[Dhia, 1997] for smooth
transition
• plenty of methods proposed
for MD – FEM coupling:
– quasicontinuum method
[Tadmor et al., 1996]
– bridging scale decomposition DEM
[Wagner & Liu, 2003]
– bridging domain method
[Xiao & Belytschko, 2004]
– … FEM

overlapping domain
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

[10-6 N]
preliminary: cone example

• reduction of fluctuations of
cone force
• no more “out of bounds” forces
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

analogous for superellipsoids

• rigid cone, non-flat cap
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Coupling of Discrete and finite Elements

 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Motivation
• Erosion and sedimentation processes in river
engineering
 Multiscale methods for sediment transport
 Homogenization strategies

• Scour genesis in front of bridge piers and

bank constructions
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Model of the Fluid – Solid Particle System

• N solid particles particles (p = 1,...,N) in a fluid (density ρ ,
viscosity μ)
• domain occupied by the fluid at time t:
• domain occupied by the particle p at time t:

Fluid phase:
• nonstationary incompressible Navier-Stokes
• Finite Element Method
Particle phase:
• Newton-Euler equations for rigid particles
• Discrete Element Method
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Fluid phase
Navier-Stokes equations in

with the stress tensor

Particle phase

Newton-Euler equations in

Translation:

Rotation:
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Coupling between fluid and particle

• Implicit coupling („Distributed Lagrange Multiplier / Fictitious
Domains“) [Glowinski et. al. (1999); Patankar et. al. (2000)]
- combined weak coupling
- mutual forces cancel

• Explicit coupling („Fictitious Boundary Method“)‫‏‬
[Duchanoy & Jongen (2003); Turek, Wan & Rivkind (2003)]

compute flow phase

compute forces on particle

compute particle phase next time step

correct flow phase

 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Forces and torque acting on the particles

Hydrodynamic
• Force:

• Torque:
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Forces and torque acting on the particles

Hydrodynamic
• Force:

• Torque:

Particle – particle contact

• Contact of spheres:
• Normal forces for contact of viscoelastic particles:

[Hertz (1881), Brilliantov et. al. (1996)]

 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Tangential forces for contact of the particles

- total relative velocity at the contact surface:

- tangential velocity:

- static / dynamic friction: ,

- projection of the tangential spring into the

actual tangential plane:

- tangential test force:

[Wriggers (1987), Luding (2007)]

 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Tangential forces for contact of the particles

- total relative velocity at the contact surface:

- tangential velocity:

- static / dynamic friction: ,

- projection of the tangential spring into the

actual tangential plane:

- tangential test force:

[Wriggers (1987), Luding (2007)]

 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Interaction between particle and fluid

• No-slip boundary condition at the interface fluid-particle

- additional constraint to the Navier-Stokes equations

 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Efficient computation of particle-fluid flow

• Instationary flow solver FEATFLOW:
- multigrid techniques [Turek (1998)]
• Computation of the particle-particle interactions
- efficient force computation, e.g. Verlet lists [Verlet
(1967)]
- efficient integration scheme, e.g. Gear's algorithmus
[Gear (1971); Allen & Tildesley (1987)]
• Computation of the hydrodynamic forces
- fast assignment of the elements of the flow field to the
particles using the different multigrid level meshes

Level 3

Level 2

Level 1
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Numerical Examples
Sedimentation of one particle in a channel
• system dimension: 2 / 2 / 6
• particle diameter: d = 0.25
• fluid density: ρf = 1.0
• particle density: ρp = 1.5
• kinematic viscosity: ν = 0.01
• gravity acceleration: g = 980

Computation
• multigrid: level 1: 7.496 elements ( 181.764 dof‘s)
level 2: 139.968 elements ( 1.426.896 dof‘s)
level 3: 1.119.744 elements (11.306.304 dof‘s)
• timestep: flow Δtf = 1.0E-4, particle Δtp = 1.0E-6
• computation time: ≈ 10h
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Sedimentation of one particle in a channel

Flow velocity Pressure field
Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010
Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Comparison with other results

1
3D-Fall Velocity
Wriggers & Avci
-1
2D-Fall Velocity Wan
-3 & Turek [2007]

-5

Velocity
-7

-9

-11

-13

-15
0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 0,45 0,5
Time
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Sedimentation of 64 particles in a channel

• system dimension: 2 / 2 / 6
• particle diameters: d = 0.25
• fluid density: ρf = 1.0
• particle denisty: ρp = 1.5
• kinematic viscosity: ν = 0.01
• gravity acceleration: g = 980
Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Driven Cavity Flow with 144 particles

• system dimension: 2 / 2 / 6
• particle diameters: d = 0.25
• fluid density: ρf = 1.0
• particle density: ρp = 1.01
• kinematic viscosity: ν = 0.1
• gravity acceleration: g = 980
• velocity at the top: u = 10.0

u = 10.0
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

Driven Cavity Flow with 144 particles

Flow field Pressure field

Side view
Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010
 Simulation of Granular Material Ljubljana
Peter Wriggers, Christian Wellmann and Bircan Avci 2010

The Team

C. Wellmann, C. Lillie and P. Wriggers,

Homogenization of granular material modeled by a three-
dimensional discrete element method, Computers and
Geotechnics, 35, 395-405, 2008

C. Wellmann, C. Lillie and P. Wriggers,

Comparison of the macroscopic behavior of granular materials
modeled by different constitutive equations on the microscale,
Finite Elements in Analysis & Design, 44, 259-271, 2008.

C. Wellmann, C. Lillie and P. Wriggers,

A contact detection algorithm for superellipsoids based on the
common-normal concept, Engineering Computations, 25,
432-442, 2008.

• B. Avci
• C. Wellmann I. Temizer & P. Wriggers:
A Multiscale Contact Homogenization Technique for the Modeling of
Third Bodies in the Contact Interface, CMAME, 198, 377-396, 2008