A coupled continuum/discrete model of granular flow
(For more information on this project, see the slides I presented at the Bay Area Scientific Computing Day.)
The simulation of flowing dense granular materials is a big challenge in physics and engineering. Despite some recent advances, there is still no complete continuum theory for accurately predicting macroscopic flow features, and for information at the level of a single particle (such as diffusion or packing structure) there are very few models available.
Perhaps the only reliable option is the Discrete-Element Method (DEM), which involves simulating individual particles according to Newton's Laws with a frictional contact model. One example of a DEM simulation is the parallel LAMMPS code developed at Sandia National Laboratories. The contact model consists of a normal elastic repulsion, a history-dependent frictional interaction, plus viscous damping. To model particles made of realistic materials (such a glass beads or sand) the elastic repulsion has a very high spring constant, leading to stiff equations which require a small timestep to simulate accurately. The large number of particles and the small timestep make DEM extremely computationally intensive. For example, a complete cycle of a pebble-bed reactor with 440,000 particles could take around a month of computation on ninety processors.
Although DEM provides an extremely accurate picture, the amount of computation clearly makes it impractical for applications in engineering. It would be good to develop a model which rapidly provides reasonable flows, and roughly captures microscopic statistics, without accurately modeling every particle. This could be used in industrial process control, where near real-time information may be required.
Results from PhD thesis
My PhD thesis had two main themes:
- In a dense granular flow, particles form an amorphous packing, and an individual particle is strongly geometrically constrained by its neighbors. In experiments and simulations, the velocity of particles are correlated with their neighbors on a local length scale, typically 3d to 5d. I developed a simulation of granular flow in a hopper, by breaking down the motion the motion into diffusing “spots” of free volume, which induce correlated motion in the underlying particles. The simulation was much faster than DEM and produced reasonably realistic flowing random packings, giving good predictions about diffusion and packing structure. However, the concept of a spot was hard to generalize to an arbitrary flow, and lacked a reasonable physical basis.
- A large challenge in creating a physical continuum model is that continuum quantities, such as stress and deformation rate cannot easily be interpreted at the level of a single particle since particle motion and forces are highly inhomogeneous. However, we were able to show that these quantities can be interpreted in an approximate sense at the scale of spot, allowing us to verify and disprove a number of continuum theories about granular materials. Based on these results, and some recent papers in the literature, Ken Kamrin was able to create a generally applicable elastoplastic continuum model. The model makes very good predictions about a number of flows, although it lacks any information at the microscopic level, and also appears to lack some stochastic effects.
Neither of these ideas represent a complete model, but perhaps these two ideas can be can combined. Can a continuum elastoplastic framework be applied down to the correlation length scale, and then a spot-like motion be used to make the individual particles flow?
A coupled continuum/discrete simulation
I have recently written a simulation attempting to combine the above two ideas. The simulation attempts to recreate the DEM simulations in the wide silo featured in the previous study, using a 150d × 8d × 70d packing which is periodic in the 8d direction. The simulation consists of a continuum solver, coupled to a model for discrete particle motion.
Part 1 – the continuum model:
The continuum model is carried out on a rectangular grid, with one grid point at the center of each computational element in the previous study (a scale of 2.5d × 2.5d). On this grid, a three-dimensional stress tensor and three-dimensional velocity are solved according to the elastoplastic model
Here, v is the velocity, sigma is the stress tensor, and a subscript 0 denotes the deviatoric part. Non-dimensionalized units are used which are similar to those in a related study. This model is not as sophisticated as the continuum model employed by Kamrin, but it is simpler to implement and captures the main features. The plastic yield is a function of the local packing fraction φ, and the ratio of shear stress to normal stress μ, and is set by
For this work, the Drucker–Prager definition of yield stress has been used, whereby
where the numerator is the magnitude of the deviatoric stress tensor, and p is the pressure. A number of different models for the dependence on the local packing fraction have been tested, with a typical one being
Initially, the velocity is set to zero, and the stress is set approximately to the values from DEM. The top surface of the packing is described using the level set method, and the level set is initialized by tracing the top of the particle packing in DEM.
The derivatives are computed using centered differences, except for the convective terms, which are computed using one-sided differences. At the boundaries, where these derivatives refer to points which are outside the simulation, interpolated ghost values are computed so that boundary conditions are upheld on the boundary. On the walls, the boundary stress is computed according to a local Coulomb yield criterion, and the normal velocity is set to zero. At the top surface, the normal components of stress are set to zero. The system is timestepped using a first-order Euler method, with the level set reinitialized at fixed intervals using a fast marching method.
Part 2 – the discrete simulation
The second half of the simulation consists of a discrete representation of the particles. The particle positions are initialized from the static packing created by DEM. At each step, the particles are moved according to a two-step process:
- Advect the particles according to the linearly interpolated continuum velocity.
- Apply a relaxation to enforce geometrical constraints: if particles are overlapped by an amount z, then they experience a normal displacement of an amount -az for some positive constant a.
Since this requires a consideration of all the particles, it is the computational bottleneck. However, it can be applied much less frequently than the continuum model. For the simulations below, the discrete update was applied once for every 250 continuum updates.
After the discrete step has been applied, the local packing fraction is recomputed, which is the fed back into the continuum model.
Preliminary results
 |
 |
| Discrete-Element Method | Continuum/discrete |
| Movie: QuickTime / AVI | Movie: QuickTime / AVI |
The above images show a comparison between DEM and continuum/discrete simulations for granular drainage. The DEM simulation is taken from the previous study. The continuum/discrete simulation was run using the method above, with an additional boundary condition that at the orifice, all the components of the stress must vanish. It is clear from the above images that there is a high degree of qualitative agreement between the two simulations. A careful calibration of the continuum/discrete model could potentially improve the match further.
 |
 |
| Discrete-Element Method | Continuum/discrete |
| Movie: QuickTime / AVI | Movie: QuickTime / AVI |
The images above show the mesoscopic stresses in granular drainage. The DEM stresses are computed from the individual particle forces within each mesoscopic element, and show smooth patterns. The continuum/discrete stresses are taken from the continuum stress field, and closely match the DEM results.
 |
 |
| Discrete-Element Method | Continuum/discrete |
| Movie: QuickTime / AVI | Movie: QuickTime / AVI |
The two images above show the packing fraction during a pushing simulation, where a layer of particles on the right half of the base is moved upwards with constant velocity. Although the DEM simulation shows a significantly wider shear band, the images are qualitatively similar. The shear dilation in the continuum/discrete simulation has not been explicitly put into the model, and is a prediction, due to the microscopic particle rearrangement.
The DEM simulation is run in parallel, and took approximately eight hours of computation time on twenty-four processors. The corresponding continuum/discrete simulation ran in 7 min 20 s on a single processor, representing a dramatic speedup, and offering the possibility of near real-time computation, when approximate results are needed.
Current problems, and future directions
While the results are promising, there are many problems that need to be addressed, and potential directions for improvement:
- An improved continuum integration – The current continuum model, featuring one grid point per element, was chosen for its conceptual simplicity. However, it is difficult with this approach to implement boundary conditions and ensure stability. Switching to a second order scheme, perhaps making use of a staggered grid arrangement, would improve the speed and reliability of the simulation.
- A better continuum model – The elastoplastic system considered above provides reasonable results, but some aspects of it are unrealistic. In particular, it can support tension, which is unphysical in a granular material. This may not be a serious issue in the current simulations, as the material always appears to be in compression, but it needs to be fixed for a general system. Using sections of the Kamrin elastoplastic model may help.
- A rigorous foundation – The current simulation has been developed from previous studies and physical principles. However, at this stage, it could greatly help to put it on a more rigorous mathematical foundation, to help improve it and understand its accuracy. Due to the complex nature of a granular system it is perhaps too much to hope for many rigorous mathematical results, but it may be possible to establish statistical bounds, when combined with some computational studies. There appear to be some parallels with other techniques. For example, for elliptic PDE problems the multigrid method creates a dramatic speedup by solving a coarse problem and then refining. For the current problem, perhaps the key issue is that coarse-grained stresses in a granular material are roughly smooth. DEM spends a large amount of time by calculating these at the particle level, which seems a waste, when their local averages are behaving quite smoothly.
- Alternative methods of coupling – The current model features only a very weak coupling between plasticity and packing fraction. However, the continuum/discrete framework offers a lot of scope for more advanced coupling effects. It is particularly promising because the elastoplastic framework currently cannot explain all features of granular flow: since flow only occurs when a yield criterion is exceeded, the model predicts plastic regions of flow, and elastic regions of no flow. This distinction is hard to make in real granular materials; for example, shear bands typically have exponentially decaying velocity profiles, with no well-defined elastic/plastic boundary, suggestive of stochastic effects. This could potentially come from a discrete coupling effect. One possibility is to enforce that the continuum velocity is perfectly tied to total motion the particles. For example, if particles are advected by v and then the relaxation causes an additional displacement w, then the continuum field should be corrected to v+w.
- A better discrete step – The particle relaxation could be viewed as a very coarse DEM step, but it is perhaps better to think of it as the enforcement of a geometrical constraint. While the current normal repulsion suffices to create valid random packings, it is unclear how closely the particle motion matches DEM results. In particular, studies of shearing suggest that the normal displacement alone will create some particle diffusion, but not enough. Perhaps accounting for friction is needed. The shearing study also shows that gravity is significant at a microscopic level, since particles can fall into gaps when they open up; this may also need to be accounted for in the particle relaxation.