Numerical Investigation of Granular Media

Introduction

In recent years, there has been a renewed interest among physicists to investigate the properties of granular matter. On the surface, this seems odd as granular media seem to have little to do with the more commonly thought of and seemingly higher prestige frontiers of physics, and clearly construction efforts within deserts, on beaches, or other sandy areas can benefit from clearer understanding of the grains' interactions. However, granular matter is ubiquitous; its study is important to industry, shipping and storing everything from grains of food to pharmaceutical powders can benifit from greater knowledge of the granular media. Also, particulate studies are helpful to understand large scale, geological events. In this study, I have implemented a simplistic simulation of granular particles and investigate their behavior over time.

Simulation

Simulation Basics

There are two popular approaches to numerically simulate granular particles. The first is to represent the particles as hard disks or spheres. For this project, however, I have followed the second approach, and have taken to particles to be soft disks interacting via a Lennard-Jones potential:


Equation 1: Lennard-Jones Potential

This potential results in the force between two particles being given as:


Equation 2: Lennard-Jones Force

For simplicity, the simulation is performed using ε as the unit of energy, σ as the unit of length, and a particle's mass as the unit of mass (at this time, I am only considering mono-disperse samples). The unit of time is that appropriate to these units of energy, length, and mass; I will use "t" to represent this temporal unit. The forthcoming third edition of An Introduction to Computer Simulation Methods [1] includes a Java simulation of particles in two dimensions interacting via a Lennard Jones potential; this program was used as a starting point.

However, the Lennard-Jones potential can result in attractive forces which is clearly non-physical for the average particulate interaction. So, by adjusting CutOff Factor / 2^(1/6) the user can create a cut-off in the potential, and thus the force. The default value of this input is 21/6σ making the force always repulsive and to smoothly go to zero at the separation resulting in the minimum of a standard Lennard-Jones potential. This gives the grains an effective radius of rcut-off = 21/6σ. The potential is slightly altered to keep the energy from having a discontinuity at rcut-off. The resulting altered potential and force are:


Equation 3: Lennard-Jones Potential altered to deal with a cut-off radius

Equation 4: The force of the altered Lennard-Jones Potential

Simply using this repulsive force between particles could lead to trouble though. If two particles came too close within the situation, or glanced off at the correct angle, the simulation could get a huge influx of kinetic energy through the interaction and quickly become meaningless. To alleviate this, the interaction is made inelastic. Whenever two particles are within rcut-off of each other, they experience the repulsion of a Lennard-Jones potential, but they also undergo a damping force given by:


Equation 5: Frictional force added to make collisions inelastic

Here, rij is the displacement vector pointing from the ith particle to the jth, and vij is the relative velocity between the two. The user-input Damping Factor allows γ to be adjusted to alter the strength of this viscous damping. Hirchfeld et al. used a slightly more realistic damping force which differed in the normal and transverse directions by more than the magnitudes of the x- and y-separations [2]. At this time, though, the simplistic damping suggested within the book seems to be adequate.

To reduce computation time, only those particles within a distance of 3rcut-off are considered when calculating the acceleration of each particle. A list is kept of all particles within this range. Originally these lists were updated whenever any particle accelerated rapidly (its kinetic energy becomes quadruple the average initial kinetic energy) or enough time has passed that an average particle could have traveled a distance given by (2.7 - rcut-off)σ/2 which is slightly more conservative than that suggested by the book. Now, the user can specify how often to calculate these lists to be on the safe side and ensure that nothing untoward is happening. Since the particles all start with a random kinetic energy (and their velocities are in a random direction), the simulation begins in a granular analogy to the gas phase allowing for the use of a skeleton of a molecular dynamics program for this project.

A problem I have been having with this simulation is that occasionally, even with a small time step and frequent computation of which particles are near each other, two particles occasionally come to close to one another and they fly off with very large kinetic energies. This results in the temperature spiking. At this point, I have put in three user inputs for different solutions to this problem and they are described below.

Inputs

The program has a total of sixteen different user-inputs. The first is Number of Particles which should be self explanatory, the default value of this value is 64. The next two inputs, Lx and Ly, give the horizontal and vertical sizes of the area the particles are contained to; both are set to 20 by default. The initial kinetic energy per particle gives the average kinetic energy of the particles at the simulation's beginning; it's default value is 10ε. The Desired KE per particle is a scaling factor used to try to combat granular collapse. If Give Random Kicks? is set to true, then at each time step, a random particle is selected and its kinetic energy is forced to take the value Desired KE per particle, the default of which is 5ε; this input should be adjustable during the simulation, but this feature does not seem to work on all computers.

The CutOff Factor / 2^(1/6) gives the radius of a grain in units of 21/6σ the default is 1 × 21/6σ. The Damping Factor gives the strength of the viscous damping between particles; the default value suggested by the book is 100. The step size of each simulation step is given by dt; by default it takes the value of 0.001t. Find Neighbor Time allows the user to specify how often the neighbor lists for each particle should be calculated; this input is in numbers of time steps. If the input is too large, the calculated value for this time (given by how far the particles could have traveled) may be used instead.

To simulate grains within a silo, Silo Set-Up? should be set to true. For a silo, the aperture at the base of the silo is given by Aperture Size. This input is in grain diameters and its default value is 5 grain diameters, that is 5 × 21/6σ. The hole is represented by a blue bar at the base of the display window.

When Have Close CutOff? is set to true, then when particles get within 0.15σ of each other on either axis, they will feel no force. This is clearly not the best way to try to deal with the problem of the kinetic energy blowing up. This is because according to [3] when granular collapse occurs, it first begins in nucleating clusters. This clumping occurs in reality due to the inelasticity of collisions between particles. In the simulation, the damping force, which makes the collisions inelastic, causes the clustering since particles do not always have enough energy to move away from each other after interacting. By shutting off the potential between particles if they get too close, this behavior may be suppressed; thus this input is set to false as a default. Artificial Temperature Reduction? allows the user to implement a bootstrap method of preventing kinetic energy explosions. If this is set to true, the default is false, then anytime a particle has a kinetic energy three orders of magnitude greater than the initial temperature, then that particle will artificially have its speed instantly reduced to that corresponding to the initial temperature. The temperature and counters of kinetic energy and potential energy are zeroed so this results in a discontinuity in the temperature versus time graph followed by a steep climb back to the system's temperature. Finally, if Slow Time Down? is set to true; the system should slow the time step to 1/100 of the input value if two particles get within 0.15σ of each other on either the horizontal or vertical axis. When they move away from each other, the time step should be restored to its original value. However, I am not sure if this input is working quite the way it should, though the system does in fact seem to behave when it is set to true. To not show every iteration of dt, steps per display can be increased above the default value of 1, but this seems to be causing trouble so it is not recommended.

Outputs

The program generates six different plots. Lennard-Jones Particles plots the position of each particle at each time step. Mean Temperature plots the average granular temperature of the system; this is discussed in more detail below. The X- and Y-Velocity Histograms tabulate the total number of times various velocity ranges have been populated. The bin size for each histogram is (2 * "initial kinetic energy per particle" / "Number of Particles"). The Speed Histogram maintains a record of the speeds of each particle and how often a given speed has been seen. Finally Number of Particles is only used for the silo situation. It plots the number of particles in the silo at any time. This can drop below 0 because I am using a feeding mechanism so the total number of particles in the display remains constant. The important feature of this graph is the slope.

At each time step, the program outputs the current temperature. When the program is stopped, it outputs the density, simply the number of particles divided by the area It also outputs the time, in terms of the number of steps taken, T and E the temperature (average kinetic energy) and total energy (sum of kinetic and potential energies) of the system when it was stopped.

More Detailed Theory

The particles' motion is calculated using the Verlet algorithm in the ODESolver of the Open Source Physics library. For non-silo situations, the simulation uses periodic boundary conditions. At each time step, each particle's x- and y- positions and velocities are updated.

Temperature

Based on the speed of the particles, the "granular temperature" is calculated. This is simply the average of the kinetic energies of all the particles at each time step. This definition of temperature does not have quite the same significance for granular particles, especially as I am modeling them, as it does for gases. Defining temperature as proportional to the average kinetic energy of gas particles "has sense only if on average the kinetic energy is the same for each molecule in the gas (or for each degree of freedom), i.e., the equipartition of energy is respected" [4]. This leads to problems for polydisperse samples as the mass of each particle is not the same, and they may have different shapes. It can also result in problems if energy can be dissipated as is the case here. This energy dissipation means that it may not be possible to arrive at relations between temperature and other properties, such as pressure and number of particles, for granular media as is common for the equations of state in thermodynamics. To analyze the kinetic energy of the granular media over time, three long runs were performed using the initial conditions. The simulation should be allowed to continue until the temperature fell to 10% of its initial value; but in practice this was never quite accomplished.

Granular Collapse

As energy is taken out of the system during collisions, eventually the particles can undergo granular collapse. That is, given enough time in the simulation, the particles can eventually come to rest (though I haven't actually seen all of the particles stop moving, the fact that the temperature can drop steadily throughout the time I have run the program shows that the system does slow down). To try to prevent this, random kicks can be given to the particles. When Give Random Kicks? is entered as true, at each time step, a particle will be chosen at random and its kinetic energy is set to Desired KE. The direction of the particle's new velocity is chosen at random. The velocities are assigned using the Box-Muller method. The Box-Muller method is explained in Section 11.5 of the book. The basic theory behind it involves setting a particle's speed equal to:

√(- 2 * KEdesired * ln(random))
Equation 6: Formula to find random kick to give to particles

This formula comes from finding the cumulative probability of the speed being between 0 and s and setting that probability equal to random. After this, s is solved for in terms of random and the above formula is arrived at through the fact that the when integrated over all values, the probability density of the speed must yield one. The direction of the speed is determined at random. This process should give a Gaussian distribution of speeds. After the system stabilizes, the three histograms can be examined to see if it actually generates a Gaussian:


Equation 7: Gaussian that should result from the random kicks

Silo Set-Up

When Silo Set-Up? is set to true, the simulation takes on a simplistic model of a silo. The periodic boundary conditions are dropped and instead, if a particle collides with a wall, it rebounds elastically. The only exception is on the base of the container. If a particle hits the line at y = 0 between (Lx - (Aperture Size)/2) and (Lx + (Aperture Size)/2), then a variable numberParticles is decremented, and the particle is allowed to move up to the top of the plot area as it would have with periodic boundaries in place. This allows the number of particles graphed to be constant for the entire simulation and thus allows it to continue indefinitely. At the same time, the counter numberParticles can be graphed to find how the flow rate out of the silo varies with the diameter of the hole at the base. According to [2] the relation should be given by:


Equation 8: Flow rate quoted by [2] through a hole in the base of a granular matter's container

In this, C is a material dependent constant of proportionality, ρ "is the grain packing density near the aperture," D is the aperture's size, d the diameter of a particle, and k a constant which must be determined based on the materials making the silo and particles. In this simulation ρ is not readily adjusted as the code is relatively slow and unwieldy as it is without trying to measure or alter local densities. During the simulation of a silo, since particles not on the base of the container are losing gravitational potential energy and gaining kinetic energy, the temperature quickly grows in magnitude and so that aspect of the simulation is not entirely meaningful for a silo.

Results

Long Term Temperature Measurement

Unfortunately, my code does not run very efficiently at this point. By my calculations it would take the simulation until approximately 68t for the temperature to drop from 10ε to 1ε. On my computer this translates to approximately 38 hours of simulation time... Up to this point, I haven't been able to leave my computer in my apartment running that long, and, well the chances of my account working in the lab for 38 minutes isn't too great. So I haven't been able to get one run to follow the system through to the point where it reaches 10% of the initial temperature using the conditions suggested in the book. I am trying to run the code for longer, and will update the page if I have any successes. Below, is a plot of Temperature versus Time of the longest runs I have been able to perform so far.


Figure 1: The longest runs accomplished at this time. The blue, purple, red, and green curves show the actual data; the black, thick curves are the average of these runs (the shorter is the average of all four runs over the length of the purple run; the longer is the average of the three longest curves over the entire length of the blue run). The thin, black curves are linear fits to the averages; their equations are shown in the chart area.

Unfortunately, I have not been able to find any reference to what functional behavior the temperature should follow at long time. [5] shows the behavior of the velocity of particles as a function of their location within a pile, but that is the closest I have found so far. However, it is clear that it should be a monotonic decrease, which is shown by the average of the runs I have performed so far. This average appears to be a linear fit so I used Excel to perform a linear regression and it appears to be a reasonably good fit.

Though the inability to completely investigate the temperature as a function of time has been discouraging, the program has been helpful in showing the "clumping" behavior associated with granular collapse:


Figure 2: Here the formation of a nucleating clump showing the beginnings of granular collapse can be seen in the top left portion of the window.

As the simulation is run, even those times in which particles go wild, a clump can generally be seen if one waits long enough. This is heartening as it implies that the program is accurately simulating at least one aspect of granular motion. Ideally, with an extremely long run, one could see virtually all of the particles collapse into one clump. This is seen to a small extent within the silo situation and this is discussed below.

Battle Against Granular Collapse

To keep the particles from collapsing upon one another, random kicks can be given to particles. The fact that these particles are selected randomly is the closest that this project comes to performing a Monte-Carlo simulation. Unsurprisingly, or at least unsurprisingly in hind-sight, the energy necessary to give each particle for the system to stabilize at 5ε must be greater than 5ε (since we're fighting against the particles slowing they have to be given a little extra something beyond their target energy). The energy given to a particle within a kick can be anything from 0ε up to the value entered as Desired KE per particleε, since the actual kick is found via a random number generator. With initial temperatures of 10ε, 5ε, or 1ε, giving kicks with amplitudes up to 5.05ε stabilized the system at 5.0ε these runs are shown below:


Figure 3: Here are the attempts to stabilize the system at a temperature of 5ε starting at temperatures of (from left to right): 10ε, 5ε, and 1ε.

After the initial fluctuations, the temperature stayed within the range of (5.00 ± 0.05)ε. The hardest situation to stabilize was the set-up starting at the target temperature. Initial fluctuations caused by the random starting conditions resulted in this set up frequently dropping down dramatically or ballooning in temperature at the beginning and it was difficult for the system to recover.

The x- and y-components of the particles velocities were examined for the case in which the initial temperature was 10ε. Since the kicks should result in a Gaussian distribution of velocities, they are plotted against the Gaussian given in Equation 7 above. These are shown here:


Figure 4: Graphs of p(vx) and p(vy), the probability distributions of vx and vy, versus Gaussians. The jagged curves are the computed probability densities, the red curves are Gaussians generated using Origin.

The probability distributions appear to be following Gaussians reasonably closely. Origin used the following equation for its fits, and the resulting values for p(vx) and p(vy) are tabulated below the equation:


Equation 8: The formula for the Gaussians generated using Origin 7.5

Parameters generated by Origin for Gaussian plots:
p(vx) p(vy)
p0 0.00052 ± 0.0011 0.00083 ± 0.00088
<v> -0.17345 ± 0.05794 0.03111 ± 0.05501
σ 4.24126 ± 0.16196 4.3086 ± 0.14609
A 0.30493 ± 0.1427 0.29911 ± 0.01205
R2 0.96817 0.97135

I had not expected the fits to be so close to Gaussians. Before running the simulation, I thought that since the particles lose kinetic energy whenever they interact that the velocity distributions would loose the Gaussian shape quickly due to the inelasticity.

Silo Emptying

The final situation to be investigated is the silo set-up. Though the grains do have effective radii of rcut-off, this only applies when they are interacting with one another. In terms of their location, they act as point masses. This means that for this set-up, k in Equation 8 should be zero since there is no friction modeled between the silo surface and the grains. In reality, of course, there is friction between grains and their containers (and this is why in hourglasses sand flows through the base at a constant rate, the walls support some of the weight of the column of sand). The spatial distribution particles after a long time for a typical silo run is:


Figure 5: Particles within the silo simulation with the hole diameter at 5d.

The particles clumping at the bottom tend to sit there without moving more than some vibration. This means that eventually the simulation could stop since all the particles could clump to either side of the aperture. This clumping also gives a hint as to why my simulation does not quite follow Equation 8; ρ, the particle density near the aperture is not constant. It starts around the same value as the bulk density in a random set-up (i.e. 0.16 particles / σ2 using the default initial conditions). A typical graph of the number of particles, this also one for a whole with diameter 5d, is shown below:


Figure 6: Number of particles as a function of time within the silo simulation.

Once data was collected for holes of diameter 10d, 10d, 7.5d, 6d, 5d, and 4d, they were plotted in Excel to find a regression for their behavior in time. The raw data was plotted as a function of time, and a regression fit was performed for each. The first derivative of each of these fits was then taken to get a rough expression for ; since the flow rate equation is not a function of time (as it shouldn't be for a steady-state flow which is what [2] describes), I found linear fits to the data; they weren't terrific fits, but were still relatively good. What I actually compared was the opposite of the first derivative, since as the program outputs the number of particles is negative. The derivatives versus particle diameter are shown below:


Figure 7: The flow rate as a function of hole diameter. The diamonds represent the information taken from the program generated data; the curve is simply a + D3/2 with a a constant picked to make the theoretical and computed derivatives have similar magnitudes.

The computed flow rate does increase with hole diameter, but since for large and small apertures, the computed flow rate is below a + D3/2, my simulation does not seem to be generating the correct flow rate. This problem is most likely caused by the fluctuating density, since I could not have a large number of particles without the simulation hanging, or the fact that the simulation does not include friction. Also, [2] points out that the equation only holds for holes with diameters greater than 5d and when the container's width is at least 2.5 times greater than the aperture size. The lower limit is not an issue here since the particles are essentially point masses when they consider the aperture. However, the upper bound could be another factor making my simulation less than completely accurate. To try to extend this portion of the project, if time allows I will try to streamline the simulation to allow for more particles or a larger container without destroying the simulation. Another extension which will be interesting if time allows is to alter the base of the silo so that it is not flat, but funnels particles into or away from the hole.

Conclusion

This project's goal was to begin work with numerical simulation of granular material. Though the simulation is not quite as accurate as I would like, as shown by the occasional explosion of kinetic energy and the container and sample size being too small to accurately recreate a silo, the program does reasonably portray the clumping behavior which proceeds granular collapse from a granular "gas." In the future the groundwork I started here may help point me in new directions for my experimental work.

References

[1] H. Gould, J. Tobochnik, and W. Christian. An Introduction to Computer Simulation Methods: Preliminary Third Edition (Addison-Wesley, 2006 – Forthcoming).
[2] D. Hirchfeld, Y. Radzner, and D.C. Rapaport, "Molecular dynamics studies of granular flow through an aperture," Phys. Rev. E 56, 4404-4415 (1997).
[3] http://www.physics.georgetown.edu/~granular/collapse.html
[4] I. Ippolito, C. Annic, J. Lemaître, L. Oger, and D. Bideau, "Granular temperature: Experimental analysis," Phys. Rev. E. 52, 2072-2075 (1995).
[5] Leonardo E. Silbert, Deniz Ertas, Gary S. Grest, Thomas C. Halsey, Dov Levine, and Steven J. Plimpton, "Granular flow down an inclined plane: Bagnold scaling and rheology," Phys Rev. E 64, 051302-1-14 (2001).

Code

Unfortunately, I have not been able to get an applet of my code working; however, if you would like to download the code for this project to peruse it, feel free to use the following links below. Also, the .class files are archived here to make it easier to access the program; just download the .class files and run them. The programs do not have package lines, and the .class files were made by compiling the code using stp.jar.
LJParticles Code
SiloParticles Code
LJParticlesApp Code
LJParticles .class file
SiloParticles .class file
LJParticlesApp .class file
Older versions of the code can be found at this website. They occasionally blow-up, but seem to behave relatively well, though they do not have as many features as those linked to above.

Updated December 9, 2005
Daniel Costantino