Transcript Granular Materials: A window to studying the Transition

Managing A Computer
Simulation of Gravity-Driven
Granular Flow
Scientific Computing Seminar, May 29, 2007
The University of Western Ontario
Department of Applied Mathematics
John Drozd and Dr. Colin Denniston
Event-driven Simulation
Collision Time Algorithm
• Tree data structure for collision times
– Ball with smallest collision time value kept at
bottom left
• Sectoring
– Time complexity O (log n) vs O (n)
– Buffer zones
Domain Sectoring
Collision Times: O(log n) vs O(n)
Racking balls (use MPI or OpenMP )
Staggered
in front
and back
as well
We only update interacting balls locally in adjacent sectors,
and we only do periodic global updates for output.
Ball-Ball Collisions
• Treat as smooth disk collisions
• Calculate Newtonian trajectories
• Calculate contact times
 
bij  rij  vij
rij t  tij   rij  vijtij    vij2 tij2  2 bijtij  rij2   2  0
• Adjust velocities




| |
2 
 vi   v j   bij /  rij  vij
pi  pi  C rij , pj  p j  C rij
pi2  p 2j  pi2  pj2 , pi  p j  pi  pj
 pi  p j  pi  pj
pi  p j  C rij  pi   C rij  p j 
pi  p j  C 2 2  C rij  p j  C rij  pi  pi  p j
0  C   C rij   p j  pi   C 
2
2
pi  pi  C rij  m vi  m vi 
 vi   bij /  rij
2
mj
mi  m j
 m rij  vij
 m bij

1   
2

2

 m bij
2
rij   vi  vi  vi  
bij

2
rij
A Smooth Ball Collision
Numerical Tricks
Avoid catastrophic cancellation, by rationalizing
the numerator in solving the quadratic formula for:
a t  b tij  c  0  v t  2 bijtij  r    0
2
ij
2 2
ij ij
2
ij
2
 b  b 2  4ac   b  b 2  4ac 
 2c
for b  0 : tij 

  b  b 2  4ac  b  b 2  4ac
2a


 b  b 2  4ac   b  b 2  4ac 
 2c
for b  0 : tij 

  b  b 2  4ac  b  b 2  4ac
2a


When comparing floating point numbers,
take their difference and compare to float epsilon:
if x1  x2  if x1  x2  float
Coefficient of Restitution
• µ is calculated as a velocity-dependent restitution
coefficient to reduce inelastic collapse and overlap
occurrences as justified by experiments and defined below

1  B vn , vn  v0
 vn   
 , vn  v0

• Here vn is the component of relative velocity along the line
joining the disk centers, B = (1)v0,  = 0.7, v0=g
and  varying between 0 and 1 is a tunable parameter for
the simulation.
1
a
0.98
0.96
0.94
0.92
0.9
6
b
5
vn
4
3
2
1
0
50
100
150
200
250
Savage and Jeffrey
J. Fluid Mech. 130, 187, 1983.
Collision rules for dry granular media as
modelled by inelastic hard spheres
Velocity
 r1'   r1  1   v n    m 2
   

 r '   r  m  m   m
1
1
2 
 2  2
r1
q
r2
m 2   r1  q 
 
 q
 m1   r2  q 
T he velocitydependentcoefficient of restit ution  v n 
determinestheenergy loss :
0.7

v 
1  1   0  n  , v n  v 0
 v n   
 v0 

0 ,
vn  v0


v2  v1   qˆ
0 
v 2  v1   qˆ
As collisions become weaker
(relative velocity vn small),
they become more elastic.
 1 , v 0  ga , a  radiu sof particle
Energyis dissipated.
C. Bizon et. al., PRL 80, 57, 1997.
300 (free fall region)

Donev et al PRL96
"Do Binary Hard Disks
Exhibit an Ideal Glass
Transition?"
250 (fluid region)
monodisperse
polydisperse
vy
200 (glass region)
dvy/dt
150
y
P
P  D  k BT1   / c 
1
y
x
z
Polydispersity means
Normal distribution
of particle radii
The density in the glassy region is a constant.
In the interface between the fluid and the glass does the
density approach the glass density exponentially?
0 = 0.9
0 = 0.95
Interface width seems to increase as 0  1
0 = 0.99
0.7
a
0.6
0.5

0.4
0.3
0.2
0.1
20.
10.
fluid
glass
vy
vy
free fall
b
15.
5.
20.
dvy dt
15.
10.
5.
0.
2.5
10
1
y
0.1
0.01
0
c
  c  A exp y
100 200 300
How does  depend on (1  0) ?
Density vs Height in Fluid-Glass
Transition
  c  A exp y
0
1
0
c
0
Log10
2
0
0
3
0
0
0
4
0
0
160
180
200
220
y
240
260
280
300
0.8
0.9
0.95
0.96
0.97
0.98
0.99
0.995
0.999
Length Scale in Transition
  c  A exp y
  1  0 
0.44
Slope = 0.42 poly
Slope = 0.46 mono
"interface width diverges"
300 (free fall region)
Y Velocity
Distribution
vy
3
2
Poiseuille flow
1
a y 250
250 (fluid region)
0 5 10 15 20 25 30
x
2
vy
1.5
1
Plug flow
snapshot
200 (glass region)
0.5 b y 200
0 5 1015202530
x
2
1.5
vy
Mono-disperse
(crystallized)
only
4
1
Mono
kink
fracture
150
y
0.5 c y 150
0 5 1015202530
x
x
z
300 (free fall region)
y 250
35
vx2
vy
20
15
vy
25
2
30
10
fluid
5
250 (fluid region)
0
5
10 15 20 25 30
x
y 235
235 (At Equilibrium
Temperature)
vx2
1.2
1
vy
0.6
vy
2
0.8
0.4
200 (glass region)
equilibrium
0.2
0
5
vx2
0.2
vy
0.3
2
0.4
vy
Granular
Temperature
10 15 20 25 30
x
y 200
glass
x
0.1
0
5
10 15 20 25 30
x
150
y
z
Fluctuating and Flow Velocity
Experiment by
N. Menon and
D. J. Durian,
Science, 275, 1997.
5
 v v
2
16 x 16
vy
1
v
2/3
flow
0.5
32 x 32
0.2
In Glassy Region !
Simulation results
0.1
0.5
1
5
10
J.J. Drozd and
C. Denniston
vf
Europhysics Letters, 76 (3), 360, 2006
"questionable" averaging over nonuniform regions gives 2/3
vy
 v  v  Tg  v y  v c 
   1 in fluid glass transition
2
x
For 0 = 0.9,0.95,0.96,0.97,0.98,0.99
Subtracting of Tg and vc and not averaging over regions of different vx2
1.5
Down centre
1.25
1
Slope  = 1.0
log10
vx2
0
0
1
0
0
0
0
2
0
0.8
0.9
0.95
0.96
0.97
0.98
0.99
0.995
Tg
Tg
0
0.75
vx2
0
0.5
log10
1
0.25
0
0.25
3
50
100
150
200
y
250
300
350
400
1.5
1
log10 vy
0.5
vc
0
0.5
"Particle Dynamics in
Sheared Granular Matter"
Physical Review Letters 85,
Number 7, p. 1428
(2000)
Experiment:
(W. Losert, L. Bocquet,
T.C. Lubensky and
J.P. Gollub)
vy
h
Velocity Fluctuations
vs. Shear Rate
240
1.7
1.6
1.5
1.4
1.3
1.2
   x v y
U
5
10
h
15 20 25
x
240 , y g 185 ,
30
0
 x y   
0.9
v2x
Tg 1 2
1.6
1.8
2
2.2
2.4
3.5
3
log 10
2.5
xv y
U
2
log 10
1.5
1
U
Slope = 0.406  0.018
From simulation
Must Subtract Tg !
Experiment Slope = 0.4
Physical Review Letters 85, Number 7, (2000)
Shear Stress
 xy
1
1
ˆ q
ˆ  xˆ q
ˆ  yˆ 

 1    r1  r2   q

t collisions 2
1
ˆ q
ˆ  xˆ q
ˆ  yˆ 
  f c 1   r1  r2   q
2
q
y
x
y 100
y 100
2
0.02
1
0.01
q.x q.y
xy
3
0
1
0
0.01
2
0.02
3
5
10 15 20 25 30
x
slope 0.15554674  x xy
1
ratio of slopes   fc 1    (r1  r2 )  q
2
5 10 15 20 25 30
x
slope  0.00165986
Viscosity vs Temperature
…can do slightly better…
   x v y
 x y   
x 16
0.5
xvy
1
log10
log10
yx
1.5
Slope ~ 2  1
1.92  0.084
   xy / x vy
2
2.5
3
0
0.9
0
0.95
0
0.96
0
0.97
0
0.98
0
0.99
3.5
1.5
1
0.5
0
2
log 10 v x
0.5
1
1.5
Tg
…"anomalous" viscosity.
Is a fluid with "infinite" viscosity
a useful description of the interior phase?
Transformation from a
liquid to a glass takes
place in a continuous
manner. Relaxation
times of a liquid and its
Shear Viscosity increase
very rapidly as
Temperature is
lowered.
Experimental data from the book:
“Sands, Powders, and Grains: An Introduction to the Physics of Granular Materials”
By Jacques Duran.
Related to Forces:
Impulse Distribution
Simulation 
Impulse defined:
Magnitude of
momentum
after collision
minus momentum before
collision.
Experiment 
(Longhi, Easwar)
Quasi-1d Theory 
(Coppersmith, et al)
Most frequent collisions
contributing to smallest
impulses
Power Laws for Collision Times
Collision time
= time between
collisions
10
P
0.1
0.001
0.00001
1.
10
7
1) spheres in 2d
2) 2d disks
3) 3d spheres
0.01 0.02
P   

0.05
0.1
0.2
0.5
15% polydispese
r grains
  2.81 0.06 2.75to 2.87in glassy region
Similar power laws for 2d and 3d simulations!
Comparison With Experiment
: experiment 1.5 vs. simulation 2.8
Discrepancy as a result of
Experimental response time and
sensitivity of detector.
 Experiment
Pressure Transducer
“Spheres in 2d”:
3d Simulation with
front and back
reflecting walls
separated one
diameter apart 

Figure from experimental paper:
“Large Force Fluctuations in a Flowing Granular Medium”
Phys. Rev. Lett. 89, 045501 (2002)
E. Longhi, N. Easwar, N. Menon
Probability Distribution for
Impulses vs. Collision Times (log scale)
P    
2
0
ln I
2
4
6
8
10
15
12.5
10
7.5
ln
5
2.5
0
 = 2.75
 = 1.50
random packing
at early stage
 = 2.75
Is there any difference
between this glass and
a crystal?
Answer: Look at
Monodisperse grains
crystallization 
at later stage
 = 4.3
Disorder has a universal effect on
Collision Time power law.
Summary of Power Laws
Radius
Polydispersity
0%
(monodisperse)
2d disks
Spheres in 2d
4
4.3
15 %
(polydisperse)
2.75
2.85
3d spheres
4
2.87
Conclusions
• A gravity-driven hard sphere simulation was used to study
the glass transition from a granular hard sphere fluid to a
jammed glass.
• We get the same 2/3 power law for velocity fluctuations
vs. flow velocity as found in experiment, when each data
point is averaged over a nonuniform region.
• When we look at data points averaged from a uniform
region we find a power law of 1 as expected.
• We found a diverging length scale at this jamming (glass)
to unjamming (granular fluid) transition.
• We compared our simulation to experiment on the
connection between local velocity fluctuations and shear
rate and found quantitative agreement.
• We resolved a discrepancy with experiment on the
collision time power law which we found depends on the
level of disorder (glass) or order (crystal).
Momentum Conservation
kik+gi = 0
Normal Stresses Along Height
 x yx   y yy   g
0
2
yy
4
6
8
10
12
0
Weight not supported
by a pressure gradient.
100
 y yy
200
y
x
300
0
400
Momentum Conservation
3
y y y,
2
1
0
1
x yx
xy
g
y 100
 x yx   y yy   g
2
3
5
 x xy
10 15 20 25 30
x
y
y 100
0.25
0.2
0.15
0.1
0.05
5
10 15 20 25
x
 0 Weightsupportedby shear stress