1. Introduction
Dune migration, landslides, avalanches, and silo instability are a few examples of systems where granular materials play an important role. Furthermore, handling and transport of these materials are central to many industries such as pharmaceutical, agricultural, mining, and construction and pose many open questions to the researchers. In spite of their ubiquity, understanding and predicting the flow behavior of granular materials is still a major challenge for science and industry. Even in a seemingly simple system such as dry sand, the presence of large numbers of internal degrees of freedom leads to highly nonlinear effects making it difficult to relate the microscopic grain‐level properties to the macroscopic bulk behavior.
Granular systems can show properties commonly associated with either solid or liquid. They can behave like a fluid, that is, yielding under an applied shear stress. On the other hand, they can also behave like solids, being able to resist applied stresses without deforming, showing also interesting anisotropic structure (contact‐and force‐networks) [1, 2]. Lucretius (ca. 98–55 B.C.) was among the first ones to recognize this interesting behavior of soil‐like materials, when he wrote “One can scoop up poppy seeds with a ladle as easily as if they were water and, when dipping the ladle, the seeds flow in a continuous stream” [3].Granular materials exhibit solid‐like behavior if the particles are packed densely enough and a network of persistent contacts develops within the medium, resulting in a mechanically stable jammed structure of the particles. On the other hand, when the grains are widely spaced and free to move in any direction, interacting only through collisions, the medium is unjammed and behaves like a fluid [4].
Due to their microscopic, discrete nature and their interesting macroscopic, bulk behavior response, granular materials are studied using both discrete and continuum mechanics frameworks. In the realm of the discrete approach, several numerical techniques that are able to reproduce the single particle motions with the given micromechanical properties of the grains have been developed. In such an approach, the dynamic behavior is studied by integrating the Newton’s equations of motion for each grain using micromechanical properties and specific interaction law. Following the pioneer work by Goldhirsch [5, 6], several numerical techniques have been developed to obtain continuum fields from discrete particle data.
Using these numerical methods, one can study the flow behavior of the idealized grains, characterized by some specific micromechanical properties, which might not exist in the nature, but is helpful in understanding the underlying physics of their global behavior. In spite of their versatile applicability and benefits, these numerical methods have limitations such as excessive computational requirements, round off or truncation errors, and an intrinsic dynamic that is sometimes not reflecting the experimental reality. On the other hand, continuum models give a macroscopic view to investigate granular material behaviors. Continuum mechanics theories solve the conservation equations for the whole medium, that is, the balance of mass, momentum, and when necessary, energy. Although the balance laws are easily deducible, defining the constitutive relations poses the bigger challenge. The latter relate stresses and strains taking into account the physics of the grain‐grain interaction.
The goal of the present book chapter is to study the constitutive behavior of granular systems using particle, numerical simulations, and micro‐macro transition. In particular, we focus on the different mechanical responses of a granular material in dense and dilute conditions, corresponding to the fluid and solid behaviors, respectively. In order to systematically analyze the influence of some crucial material parameters, which affect the flow behavior, we focus on an idealized material composed of frictionless, spherical particles, in the absence of any interstitial fluids. Moreover, in order to concentrate on the rheology of particulate systems, disregarding boundary effects, we have considered two system setups which allow simulating steady and homogeneous flows.
This chapter is organized as follows. Section 2 introduces the general rheological framework to describe the flow behavior of granular materials. In the same section, we also briefly review some existing granular rheological models. The particle simulations along with micro‐macro transition are introduced in Section 3, where different system setups that are used to study the steady and homogeneous granular flows are shown. Finally, in Section 4, we present a comprehensive comparison of the existing simulation data with frictionless particles in dilute and dense regimes. In the same section, we highlight the effect of various micromechanical properties (coefficient of restitution, polydispersity, and particle stiffness) on the macroscopic fields (stresses and volume fraction). We present a comparison of these results with the theoretical models in two regimes: the kinetic theory in the dilute regime, and a recently proposed generalized rheological model in the dense regime.
2. Granular rheology
2.1. A micromechanical based continuum approach
Despite the fact that granular materials are discontinuous media, their behavior is commonly described by a continuum approach. Continuum mechanics theories solve the conservation equations of the whole medium, that is, the balance of mass, momentum, and when necessary, energy. Although the balance laws are easily deducible, the big challenge is the definition of the constitutive relations, that is, the rheology. The latter captures the macroscopic behavior of the system, incorporating the microscale grain‐grain interaction dynamics.
A granular flow can undergo different behaviors depending on both properties at the particle level and the macroscopic characteristic of the flow (i.e., velocity and concentration). At the microscopic level, each particle is characterized by its shape, dimension, material, and contact properties. For the sake of simplicity, in this chapter an assembly of identical spheres, of diameter d, density ρ_{p}, and equivalent linear contact stiffness k_{n} is considered. The density of the continuum medium can be computed as the product of the particle density and the volume fraction, ν, defined as the fractional, local volume occupied by the spheres: ρ = ρ_{p} ν. Given that each grain i moves with velocity v_{i}, the macroscopic velocity of N‐particles flow in a volume V can be defined as the average
In the framework of continuum mechanics, dimensionless numbers are often introduced in order to describe the material behavior. These dimensionless numbers are defined as the ratio of different time scales or forces, thus signifying the relative dominance of one phenomenon over another.
In the case of granular flows, the macroscopic time scale associated with the shear rate parallel to the flow plays an important role. Then, it is convenient to scale all the quantities using the particle diameter, particle density, and shear rate
2.2. Continuum models
In the early modeling attempts, granular flow is envisaged as existing in either dense solid‐like or loose gas‐like regimes. Early works using shear cell experiments observed these regimes by varying the shear rate and allowing the bed to dilate or compact. Granular materials exhibit solid‐like behavior if the particles are packed densely enough and a network of persistent contacts develops within the medium, resulting in a jammed mechanically stable structure of the particles. On the other hand, when the grains are widely spaced and free to move in any direction, interacting only through collisions, the medium is unjammed and behaves like a fluid [7].
In the fluid‐like limit, the system is very dilute and the grains interact mainly through binary, instantaneous, uncorrelated collisions. One of the first rheological models for granular flows in this regime was proposed in 1954 by Bagnold [8]. This empirical model, derived from experiments in two‐dimensional plane shear flows, basically states that the stresses are proportional to the square of the strain rate. This simple law, now known as “Bagnold scaling,” has been the first to understand the physics of granular dynamics at large deformations and has been verified for dry grains in a number of experimental and numerical studies [9–12]. In the fluid‐like regime, the generalization of kinetic theory of granular gases provides a meaningful hydrodynamic description.
On the other hand, when the system is very dense, its response is governed by the enduring contacts among grains, which are involved in force chains; the deformations are extremely slow because the entire network of contacts has to be continuously rearranged (jammed structure). In these conditions, the granular material behaves like a solid, showing an elastic response in which stresses are rate independent. The corresponding flow regime is usually referred to as quasi‐static. Slowly deforming quasi‐static dense granular material has been mainly investigated in the framework of geo‐mechanics. There, the majority of the constitutive models are based on the theories of elasto‐plasticity and visco‐plasticity [13–16], and many of them have been conceived by starting from the well‐known critical state theory [17, 18].
In the transition phase, the grains interact via both force chains and collisions. None of the models cited above is able to deal with this phase‐transition of granular materials from a solid‐like to a fluid‐like state and vice‐versa. Intensive studies of the granular rheology at the phase transition have been conducted in the last decades, for example, by Campbell [19], Ji and Shen [20, 21], and Chialvo et al. [22] using 3D simulations of soft frictional spheres at imposed volume fractions. In these works, the authors derived a flow‐map of the various flow regimes and analyzed the transition areas. In particular, they found that, for a collection of particles, the solid‐fluid transition occurs in the limit of zero confining pressure at the critical volume fraction ν_{c}. Then the solid‐like regime, in which stresses are independent of shear rate, occurs for volume fractions ν > ν_{c}, whereas, at volume fractions ν < ν_{c} the system shows a fluid‐like behavior with stresses scaling with the square of the shear rate. In the proximity of the critical volume fraction, a continuous transition between the two extreme regimes takes place, for which the rheological behavior is still not fully understood.
More recently, new theories have been developed to model the phase transition. The French research group GDR‐MiDi [23] has suggested that dense granular materials obey a local, phenomenological rheology, known as μ(I)‐rheology, that can be expressed in terms of relations between three nondimensional quantities: volume fraction, shear to normal stress ratio, usually called μ, and inertial parameter I. The latter is defined as the ratio of the time scales associated with the motion perpendicular and parallel to the flow:
Below we present a summary of the two continuum theories that well describe the flow behavior in the limits and their extension to the intermediate regime. Kinetic theory in its standard form (SKT) provides a meaningful hydrodynamic description for frictionless particles in the very dilute regime, while μ(I)‐rheology holds for both frictionless and frictional particles for dense flows. It is important to mention that both theories work only for ideal systems, made of rigid, perfectly elastic, monodisperse particles. Finally, the extension of μ(I)‐rheology to deal with soft and deformable particles is also introduced.
2.2.1. Standard kinetic theory (SKT)
This section is largely based on the notable works of Brilliantov et al. [33], Garzo et al. [34, 35], Goldhirsch [6, 36], and Pöschel et al. [37].
The term “granular gas” is used in analogy with a (classical) molecular gas, where the molecules are widely separated and are free to move in all directions, interacting only through instantaneous, uncorrelated collisions. The main differences between molecular and granular gases are that in the latter case part of the energy is irreversibly lost whenever particles interact and the absence of strong scale separation. These facts have numerous consequences on the rheology of granular gases, one of which being the sizeable normal stress differences [38].
Analogous to the molecular gases (or liquids), the macroscopic fields velocity and mass density are defined for granular systems [6]. An additional variable of the system, the granular temperature, T, is introduced as the mean square of the velocity fluctuations of the grains, in analogy with molecular gases, quantitatively describing the degree of agitation of the system.
Following the statistical mechanics approach, the kinetic theory of granular gases rigorously derives the set of partial differential equations given by the conservation laws of mass, momentum, and energy (the latter describing the time development of the granular temperature) for the dilute gas of inelastically colliding particles.
In this section, we summarize the standard kinetic theory (SKT) for the case of steady and homogeneous flows for a collection of ideal particles, that is, they are rigid, monodisperse, frictionless with diameter, d, and density, ρ_{p}. In this case, the mass balance is automatically satisfied, the momentum balance trivially asserts that the pressure, p, and the shear stress, τ, are homogeneous and the flow is totally governed by the balance of energy, which reduces to
where Γ is the rate of energy dissipation due to collisions and γ is the shear rate. The constitutive relations for p, τ, and Γ are given as [39]
where, f_{1}, f_{2}, and f_{3}, are explicit functions of the volume fraction ν and the coefficient of restitution, e_{n}, (ratio of precollisional to postcollisional relative velocity between colliding particles in the normal impact direction), and are listed in Table 1.

Further, by substituting the constitutive relations for τ and Γ into the energy balance, the granular temperature drops out, so that the pressure becomes proportional to the square of the strain rate (Bagnold scaling [8])
SKT was rigorously derived under very restrictive assumptions. In particular, the granular system is assumed to be monodisperse and composed of spherical, frictionless, and rigid particles, interacting only through binary and uncorrelated collisions [7, 40, 41]. Several modifications to the SKT have been introduced in the literature accounting for different effects: interparticle friction [4, 7, 42–44], nonsphericity [45], or polydispersity [46]. As one example, Jenkins [47, 48] extended the kinetic theory to account for the existence of correlated motion among particles at high concentration.
2.2.2. Traditional µ(I) rheology
A convincing, yet simple phenomenological model that predicts the flow behavior in moderate‐to‐dense regime is the µ(I) rheology. Once again, this rheological law is based on the assumption of homogeneous flow of idealized rigid, monodisperse particles, though the extra constraint of frictionless particles can be dropped. According to this empirical model, only three dimensionless variables are relevant for steady shear flows of granular materials: the volume fraction ν, the shear stress to normal stress ratio µ = τ/p, and the inertial number I [9, 23, 28]. The collaborative study GDR‐Midi showed the data collapse for various shear geometries such as inclined plane, rotating drum, and annular shear when analyzed in terms of the inertial number. µ(I) rheology in the standard form is given by
with µ_{0}, µ_{∞}, and I_{0} being dimensionless, material parameters which are affected by the micromechanical properties of the grains [49].
To account for the polydispersity of particles, the generalized inertial number taking into account the average diameters of the particles was introduced by [50]. Traditional µ(I) rheology had been successful in describing the flow behavior of homogeneous flows (both dense and fast). But it has failed to capture the slow and nonhomogeneous flow, where a shear rate gradient is present. Researchers have made significant efforts into developing nonlocal models for granular flows [51].
2.2.3. Soft µ(I) rheology
When particles are not perfectly rigid, instead they have a finite stiffness (or softness), the binary collision time is nonzero and hence presents an additional timescale, which is ignored in the standard inertial number phenomenology. A dimensionless number signifying the finite softness of the particles is the dimensionless pressure
with the dimensionless pressure p^{*} being the characteristic pressure at which this correction becomes considerable.
The other dimensionless number needed for the full flow characterization is the volume fraction ν. In case of rigid particles under shear, the packing will dilate and hence ν depends only on the inertial number I. On the other hand, a packing made up of soft particles will dilate due to shear, at the same time pressure will lead the compression of the particles. Hence ν depends on both I and p^{*} as
where I_{c} and p_{c}^{*} are material dependent dimensionless quantities [49, 52] and ν_{c} is the critical volume fraction, governing the fluid‐solid transitions introduced in the previous section. Its dependence on the polydispersity of the system will be discussed in Section 4.
3. Numerical simulations
Since a few decades, dynamic particle simulations have been a strong tool to tackle many challenging issues related to understanding the flow behavior of particulate systems.
The molecular dynamics or discrete element methods (DEM) is the term given to the numerical procedure, which is used to simulate assemblies of discrete particles. Molecular dynamics (MD) was originally introduced to simulate the motion of molecules [53–55]. It is essentially the simultaneous numerical solution of Newton’s equation for the motion of individual particles, for which the position, velocity, and acceleration are computed at each time step. Through averaging of positions, velocities, and forces of the particles, the macroscopic fields of the whole system, such as the density, mean velocity, and stresses can be obtained in terms of the micromechanical properties. This helps in revealing insights of the behavior of granular materials, which cannot be captured by experiments. In particular, with MD methods, the role of micromechanical properties of the grains on the macroscopic collective behavior of the system can be analyzed.
Particle simulation methods include three different techniques: The discrete element method (DEM), the event‐driven (ED), and the contact dynamics method (CD). All these methods simulate the inelastic and frictional nature of the contacts among grains through microscopic coefficients (i.e., the coefficients of restitutions and the interparticle friction coefficient). In DEM, deformations of particles during contacts are modeled allowing a finite overlap between grains, whereas in the other two methods, the particles are assumed to be infinitely rigid. Since the results presented in this chapter are obtained by using DEM simulations, below we briefly present an overview of DEM. Readers interested in the latter two methods are referred to Refs. [56–58].
3.1. Discrete element method (DEM)
The discrete element method (DEM) is a family of numerical methods for simulating the motion of large numbers of particles. In DEM, the material is modeled as consisting of finite number of discrete particles, with given micromechanical properties. The interactions between particles are treated as dynamic processes with states of equilibrium developing when the internal forces balance. As previously stated, the granular material is considered as a collection of discrete particles interacting through contact forces. Since the realistic modeling of the deformations of the particles is extremely complicated, the grains are assumed to be nondeformable spheres which are allowed to overlap [58]. The general DEM approach involves three stages: (i) detecting the contacts between elements; (ii) calculating the interaction forces among grains; and (iii) computing the acceleration of each particle by numerical integrating the Newton’s equations of motion while combining all interaction forces. This three‐stage process is repeated until the entire simulation is complete. Based on the fundamental simulation flow, a large variety of modified codes exist and often differ only in terms of the contact model and some techniques used in the interaction force calculations or the contact detection.
In this chapter, we focus on the standard linear spring‐dashpot (LSD) model. Considering two particles, i and j, of diameter d and density ρ_{p} (i.e., mass m = ρ_{p}πd^{3}/6), their contact leads to the normal (in the direction connecting the centers of the two particles in contact) and tangential components of forces as
where
Collisions may be described using the coefficients of normal and tangential restitution, e_{n} and e_{t}, respectively, relating the pre‐collisional and post‐collisional relative velocities. For the spring‐dashpot model, the following relations between the coefficients of restitution, the spring constants and the damping coefficients hold [59]
3.2. Micro‐macro transition
A research goal in the granular community is to derive macroscopic continuum models based on relevant micromechanical properties. This means to bridge the gap between the microscopic properties and the macroscopic mechanical behavior. The methods and tools for this so‐called micro‐macro transition are often applied to small so‐called representative volume elements (RVEs), where all particles can be assumed to behave similarly. Note that both time‐ and space‐averaging are required to obtain reasonable statistics, the latter being appropriate in the case of steady states.
As previously introduced in Section 2.1, the average velocity of N particles in the RVE V gives the macroscopic velocity u, while the strain‐rate tensor involves the velocity gradient of the particles
being v_{i} the velocity of particle i. For the particular case of granular systems with mean flow in the x‐direction only and subjected to shear in the y‐direction, the shear rate is introduced as
The stress tensor is of particular interest for the description of any continuum medium. In the case of granular assemblies, previous studies have proposed stress‐force relationships for idealized granular systems that relate average stress in the assembly to fundamental parameters that are explicitly related to statistical averages of inter‐particle load transmission and geometrical arrangement. When referring to a homogeneous volume element V, the macroscopic stress tensor σ can be calculated as
where F_{ij} is the contact force and l_{ij} the branch vector in between connecting the centers of particles i and j, and V_{i} = v_{i} ‐ u is the velocity fluctuation of particle i. The first and second terms in the previous equation represent the dynamic and static contributions, respectively [5, 60]. The pressure and shear stress are finally defined as
where
3.3. Simulation setups
There are two popular ways to extract continuum quantities relevant for flow description such as stress, density, and shear rate from the discrete particle data. The traditional one is ensemble averaging of “microscopic” simulations of homogeneous small samples, a set of independent RVEs. A recently developed alternative is to simulate a nonhomogeneous geometry where dynamic, flowing zones and static, high‐density zones coexist. By using adequate local averaging over equivalent volume (inside which all particles can be assumed to behave similarly), continuum descriptions in a certain parameter range can be obtained from a single simulation.
In Section 4 we will combine results from (a) simple shear RVE and (b) split‐bottom shear cell. The setups are briefly introduced and shown in Figure 1 (see Refs. [30, 49] for more details) and relevant numerical parameters are reported in Table 1. When dimensionless quantities (see Section 2.1) are matched and averaging zones are properly selected, the behaviors from different setups are comparable and a wide flow range can be explored.
3.3.1. Simple shear RVE
The collection of spheres of mean diameter d and density ρ_{p}, sheared under steady conditions is considered. Here and in the following, x and y are taken to be the flow and the shearing directions, respectively, and variations along the transversal direction z are ignored. We also introduce the polydispersity w as the ratio between the maximum and the minimum particle diameter. In this simple configuration, the flow is assumed to be one‐dimensional such that the horizontal velocity u_{x} is the only nonzero component, and the stress tensor reduces to two scalars; the pressure p and the shear stress τ. In the steady state, the mass balance equation is automatically satisfied and the divergence of the velocity is zero. The momentum balance equation, in absence of external forces, indicates that both pressure p and shear stress τ are constant. Simple shear flows are homogeneous if the horizontal velocity of the medium varies linearly along the gradient direction and the dominant kinematic variable is its first spatial derivative, the shear rate,
Variables governing the problem are the volume fraction ν (also known as density/concentration defined as the fraction of volume occupied by the spheres), the pressure p, and the shear stress τ. Using DEM simulations, we have performed simulations by using two types of simple shear experiments, that is, (i) constant pressure (here refers to normal stress) or (ii) constant volume boundary conditions. In the former (Figure 1b), pressure and strain rate are held constant, hence density and shear stress are outputs and the system is free to dilate/compact based on the initial volume fraction of the packing. In case of constant volume (Figure 1a), volume fraction and shear strain rate are held constant, so that pressure and shear stress are the outputs. Constant pressure is one of the traditional methods used in the soil mechanics to estimate the shear strength of the material, while constant volume method is used often to understand the flow behavior close to the jamming transition. Shearing under constant‐volume conditions is difficult to realize experimentally due to the fundamental characteristic of the behavior of granular materials, however, a pertinent experiment would be the undrained shear test on water‐saturated sand where the volume of the whole specimen can be kept constant within the range of experimental error [18]. On the other hand, dense granular flows under constant stress are present under experimental or natural conditions, for example, sand or/and powders sheared in different shear cells [61] or in an avalanche [62].
Constant‐volume steady simple shear samples are placed in a cuboid box (Figure 1a). The height of the computational domain as H = 20d, with d particle diameter, is fixed before we compute the x‐ and z‐size L according to the chosen, fixed, volume fraction ν. Simulations have been performed using a monodisperse system (w = 1) by systematically changing both the volume fraction ν, ranging from dilute to dense regime and the particle stiffness k_{n} such that the dimensionless shear rate γ(ρ_{p}d^{3}/k_{n})^{1/2} ranges from 3 × 10^{−2} to 3 × 10^{−4}.
In the case of RVE under constant normal stress condition (Figure 1b), granular systems with polydispersity w = 2 and w = 3 are considered. The initial length of side is set to L, along with the center point in x‐y‐plane (marked as O), where one always has zero mean field shear velocity during the whole simulation. The normal stress σ_{yy} is kept constant along y‐direction. In this way, the sample is free to dilate/compact along y‐direction and smoothly reaches its steady state. In order to investigate the sheared granular flow behavior with different inertia and particle stiffness, we systematically vary both the confined normal stress σ_{n} and shear strain‐rate γ such that the dimensionless stress/softness σ_{yy}(d/k_{n}) ranges between 10^{−3} and 10^{−1} and the dimensionless shear strain‐rate
3.3.2. Split‐bottom ring shear cell
A common feature of natural slow granular flow is the localization of strain in shear bands, which are typically of few particle diameters width. A specialized geometry proposed recently which allows one to impose an external deformation at constant rate is so‐called split‐bottom geometry (Figure 1c). In this geometry, stable shear bands of arbitrary width can be achieved allowing for a detailed study of microstructure associated with the flow of granular materials in the steady state. Unlike the previous setups, in the split‐bottom geometry, the granular material is not sheared directly from the walls, but from the bottom. The bottom of the setup that supports the weight of material above it is split in two parts, the two parts move relative to each other and creates a wide shear band away from sidewalls. The resulting shear band is robust, as its location exhibits simple and mostly grain independent properties.
In this geometry, due to inhomogeneous flow, granular packings with contrasting properties and behavior coexist, that is, high‐density static to quasi‐static areas and dilated dynamic flowing zones are found in the same system. A superimposed grid meshes the granular bed and averaging is performed within each grid volume. Inside a grid volume all particles are assumed to behave similarly and information for a wide parameter range can be obtained using a single numerical experiment, for example, at increasing pressure levels along the depth of the cell. In the following sections, when presenting data from split‐bottom cell simulations, only grid‐points in the center of the shear band will be considered, where the shear rate γ is higher than a given threshold (see Refs. [3, 30–32, 63] for details on the data processing). Data in center of the shear band are not affected by boundary effects, so that flow gradients can be neglected and the system can be considered as locally homogeneous. In the split‐bottom geometry, the shear rate γ is computed as a function of the relative angular velocity Ω between inner and outer cylinders. Details on the geometry setup and numerical parameters adopted for the simulations described in the following section are reported in Table 2.
4. Rheological flow behavior
In this section, we compare the results from various flow setups discussed above for low‐to‐high volume fractions. We vary various particle and contact properties to understand how the particle micromechanical properties influence the macroscopic flow behavior. We have compared different datasets from different setups and/or authors, and numbered as follows: [A] Peyneau et al. [64]; [B] Chialvo and Sundaresan [65]; [C] Shi et al. (unpublished); [D] Singh et al. [30, 63], and [E] Vescovi and Luding [49]. Unless specified, we will only use the data labels in the following discussion for the sake of brevity.
4.1. Influence of coefficient of restitution
Figure 2 presents a data collection from two different setups and plots the dimensionless pressure against volume fraction. It shows data with constant pressure simulations from data [A] together with the constant volume simulation results of data [B], for frictionless monodisperse rigid particles. As expected, the data from the two setups are in good agreement. We observe that the restitution coefficient e_{n} affects the dimensionless pressure strongly for volume fractions ν < 0.6, which increases with increase in e_{n}. However, in the high volume fraction limit, the data for different e_{n} collapse on the limit curve diverging at ν_{c}, that is, ν ranging between 0.6 and the critical volume fraction ν_{c}.
For the dilute case, a granular gas with high restitution coefficient, for example, e_{n} = 0.99 will behave nearly like an ideal gas, that is, almost no energy loss during each particle‐particle collision. Hence, the system will reach equilibrium with higher fluctuation velocity (proportional to the dimensionless pressure) for each particle. In the other extreme, for a restitution coefficient equal to 0, the particles lose all their energy at one collision. Such strong dissipation leads to a rather small pressure in the system. As ν approaches the critical volume fraction, for rigid spheres, the mean free path available for particles decreases making it more difficult to move the particles by imposing shear. The frequency of the collisions and thus the pressure both increase since the free path decreases, diverging in the limit case. Once one reaches the critical volume fraction limit, the system is jammed, hence shear movement of particles without further deformation is not possible. The increase of the pressure for decreasing volume fraction (below 0.1), as the probability of collisions is reduced in the dilute case, is due to the collisional energy loss with a higher steady state pressure. As for the standard kinetic theory prediction, it captures the behavior below volume fractions 0.5 well, but fails for higher volume fractions. This is expected because the standard kinetic theory (SKT) does not take the critical volume fraction into account and thus leads to an underestimation of the pressure for high volume fractions.
4.2. Influence of polydispersity
Figure 3 shows the variation of the nondimensional pressure with volume fraction for different polydispersity for constant pressure (data [A] and [C]), constant volume (data [B]) homogeneous shear flow simulations, together with the local shear band data from nonhomogeneous shear flows (data [D]). We observe that for low‐to‐moderate volume fractions, pressure is weakly increasing with volume fraction. The data from different shear setups and different polydispersity collapse and agree with the predictions of SKT. However, for higher volume fractions (ν > 0.55), pressure increases when approaching ν_{c}. However, different polydispersity yields different ν_{c} [66], so that the pressure decreases with increase in polydispersity, due to the increase in free space available for particle movement for higher polydispersity (in the cases studied here). In some cases, the small particles can move into the gaps between larger particles and form rattlers (rattlers do not contribute to the pressure as for mechanically stable contacts). Therefore, the critical volume fraction ν_{c} increases with increase in polydispersity as shown by the vertical dashed lines, consistent with previous studies [66–68]. Note that the shear band data from nonhomogeneous split‐bottom setup (data [D]) has more scattered than the others, due to the fluctuations of the local averaging over small volumes. But most of the data still follow exactly the same trend as the homogeneous shear data for same polydispersity. We also note that some data points, for example, for polydispersity w = 3, go beyond the critical volume fraction due to the fact that DEM particles are not infinitely rigid (they have large but finite stiffness). This softness (and hence possibility of deformation) leads to flow above ν_{c} and will be elaborated next.
4.3. Effect of particle stiffness
In Figure 4, we show the dimensionless pressure as a function of volume fraction for various values of dimensionless particle stiffness, ranging from 10^{3} to 10^{7}. The vertical dashed line shows the monodispersed critical volume fraction as in Figure 3. For the sake of comparison, rigid cases (data [A] and [B]) are also plotted. As expected, for the rigid case, pressure diverges close to the critical volume fraction. For soft particles, the deviation from the rigid case is a function of particle stiffness and depending on the system volume fractions (even for the softest particles the deviation from the rigid limit is small for volume fractions smaller than 0.55). When decreasing the volume fraction below 0.5, all different stiffness data tend to collapse. The solid line is the same standard kinetic theory as in Figure 3 where the assumption of rigid particle breaks down for volume fractions ν > 0.5. And the horizontal dashed line is the prediction from extended rheological model in Eq. (6) using the fitting parameters taken from Ref. [49] for the data with dimensionless particle stiffness 10^{5}. Our new extended dense rheological model smoothly captures the soft particles behavior even beyond the critical volume fraction and works perfectly between volume fraction 0.3 and 0.7.
4.4. Combining both particle stiffness and polydispersity in the dense regime
Figure 5 displays dimensionless pressure plotted against volume fraction for both constant volume (data [E]) and normal stress (data [C]) setups with three polydispersities and dimensionless contact stiffnesses, in the moderate to dense volume fraction regime. Diamonds represent constant volume simulation for monodisperse particles while stars and triangles refer to the constant pressure simulation data for polydispersity 2 and 3, respectively, and different color represent different particle stiffness. For ν < 0.55, the data points from the two setups collapse and following the same trend as for the rigid case (Figure 3, data [A]). Interestingly, for the data above the critical volume fraction ν_{c}, the pressure data for different polydispersity are found to collapse with a given dimensionless stiffness (both for 10^{5} and 10^{7}). This indicates that once the system is jammed, the particle stiffness (deformation) determines the pressure without much effect of the polydispersity of particles. The solid and dashed lines are the same lines as in Figure 4, but is given there as a guidance to the eye representing a reference to the connections. We observe the SKT solid line is not predicting the behavior at all while the extended dense rheology dashed line is qualitatively capturing the behavior even for volume fractions ν > 0.7, but with considerable deviations. Note that there are small differences between the data from two setups and it is due to the small differences in the particles stiffness, and this will be elaborated in the next section.
4.5. From dilute to dense, from “liquid” to “solid,” universal scaling
Figure 6 shows the pressure nondimensionalized in two possible ways (a) using shear rate and (b) using particle stiffness (as introduce in Section 2.1) plotted against the distance from the critical volume fraction for the data from different simulations using frictionless particles. Figure 6a shows a good data collapse for the volume fractions below the critical volume fraction (unjammed regime), or the so‐called fluid regime. In the special case of nearly rigid particles or small confining stress, the scaled pressure diverges at the critical volume fraction, which indicates that the granular fluid composed of rigid particles under shear cannot reach a denser shear jammed state. For the data with softer particles, flow is possible even above the critical volume fraction. For low to moderate volume fractions, the agreement of our data with the rigid case is excellent, while for high volume fractions (especially close to the critical volume fraction) deviations are considerable. The data collapse in the low volume fraction regime shows that the Bagnold scaling relationship between pressure and volume fraction is not strongly affected by particle stiffness, polydispersity, and shear setups, but was influenced by the restitution coefficient (see Figure 2). The “fluid” experiences the energy loss more prominent due to collisions.
For larger volume fractions, the scaling does not collapse the data. Note the deviation between constant volume (data [E]) and constant pressure (data [C]) due to the small difference in the dimensionless stiffness as shown in the legend.
Figure 6b shows the same data but only the soft particle simulations ([C] and [E]) with pressure nondimensionalized by the particle stiffness. In this way, we observe a data collapse for high volume fractions, ν > ν_{c}, in agreement with the rate independent behavior as observed in other studies. This collapse of data for ν > ν_{c} indicates that above the critical volume fraction the steady state rheological behavior of soft granular media under shear is dominated mostly by particle stiffness, while the influences of polydispersity and restitution coefficient (e_{n} = 0.8 in data [C] and e_{n} = 0.7 in data [E]) are of minor importance. In this regime, the higher the volume fraction the more solid like the behavior, and hence the less influences come from other microparameters than stiffness. It is also important to mention that even though we presented the analysis for pressure only, the shear stress shows a similar quantitative behavior [49].
4.6. So much for the granular rheology
While up to now, the focus was on understanding the relation between pressure and volume fraction, a granular rheology also must consider the shear stress.
Figure 7 shows the steady state shear stress ratio, μ = τ/p (scaled by pressure, mostly referred as macroscopic friction), against inertial number for all the data discussed from Figure 6a (with different polydispersity, restitution coefficient, particle stiffness, as simulated in diverse numerical setups). It is important to realize that though both shear stress and pressure diverge close to the critical volume fraction point, their ratio does not. We observe the traditional μ(I)‐rheology as a basic trend. For low inertial number, μ is almost independent of, I, and increases with increasing, I, for intermediate to large, I. Interestingly, although the qualitative trend of all the data is predicted by the traditional rheology, we still observe the deviations from the prediction in Figure 7. There are still many unveiled folders in the granular rheology like nonlocal behaviors, small shear rates diffusion, particle softness influence, etc., not to mention the complexity of including the frictional and cohesive granular media or/and with liquid bridges and suspensions. And also, the missing link between the dilute and dense granular rheological models is still a great challenge in the future.
5. Conclusion
This chapter gives an overview of recent progress in understanding and theoretically describing the collective mechanical behavior of dissipative, deformable particles in different states, both fluid‐like and solid‐like. Particulate systems and granular matter display collisional, dilute and solid, mechanically stable states, either switching forth and back, or both at the same time. In which state the system resides depends not only on material properties like, for example, their discrete nature (elastic stiffness), the dissipation (restitution coefficient) or the size distribution (polydispersity) of the particles, but also on the density of the system and balance between the energy input by (shear) stress or strain‐rate and the energy dissipation by collisions or plastic deformations. Realistic material properties like friction and cohesion as well as nonsphericals particles go beyond the scope of this chapter.
One extreme case of low and moderate density collisional flows (for weak to moderate dissipation and arbitrary polydispersity) is well described by standard kinetic theory (SKT) up to system volume fractions about 0.5, beyond which the elastic behavior of longer‐lasting contacts becomes dominant. Open challenges involve very soft particles for which basic theoretical assumptions of kinetic theory fail, for example, due to multiple contacting particles.
The other extreme case of quasi‐static flow of elastic, mechanically stable solid‐like structures are approximately described by the classical μ(I)‐rheology in the limit of rigid particles, but require a softness correction for comparatively large confining stresses. Remarkably, dissipation, as quantified by the coefficient of restitution, dominates the collisional flows in the dilute regime, while the particle stiffness, the polydispersity, and the friction (data not shown here) are the controlling microparameters for denser quasi‐static and jammed flows.
The mystery of bridging the gap between the collisional, dilute, and the denser quasi‐static, elastic solid‐like regimes is not completely solved yet. The particulate, microscopic states are well understood by particle simulations that via so‐called micro‐macro transition can guide the development of macroscopic, continuum constitutive relations that allow to predict the state and characteristics where a granular system resides in. A unified description that ranges from dilute to dense, from rapid to slow, from soft to rigid, etc., is still one of the great challenges of today’s research.
This chapter provided a few methods and some phenomenology, as well as an overview of recent literature in this field, with theories that can describe the extremes. Various recent works attempted to combine those limit‐cases and provide first combined, generalized theories that go beyond the classical states. However, due to dissipation, friction, cohesion, and nonsphericity of realistic materials, this poses still plenty of challenges for today’s research. Our own ongoing research focuses on providing simple unified/generalized theories, also for systems with attractive forces and with anisotropic microstructures, which were not addressed in this chapter.