1 Introduction

Simple homogeneous shear flows of frictionless, deformable particles are studied by particle simulations at large shear rates and for differently soft, deformable particles. The particle stiffness sets a time-scale that can be used to scale the physical quantities; thus the dimensionless shear rate, i.e. the inertial number (inversely proportional to pressure), can alternatively be expressed as inversely proportional to the square root of the particle stiffness. Asymptotic scaling relations for the field variables pressure, shear stress and granular temperature are inferred from simulations in both fluid and solid regimes, corresponding to unjammed and jammed conditions. Then the limit cases are merged to unique constitutive relations that cover also the transition zone in proximity of jamming. By exploiting the diverging behavior of the scaling laws at the jamming density, we arrive at continuous and differentiable phenomenological constitutive relations for the stresses and the granular temperature as functions of the volume fraction, shear rate, particle stiffness and distance from jamming. In contrast to steady shear flows of hard particles the (shear) stress ratio does not collapse as a function of the inertial number, indicating the need for an additional control parameter. In the range of particle stiffnesses investigated, in the solid regime, only the pressure is rate independent, whereas the shear stress exhibits a slight shear rate- and stiffness-dependency.

Merging fluid and solid granular behavior


Dalila Vescovi and Stefan Luding


Multi Scale Mechanics (MSM), CTW, MESA+, University of Twente, PO Box 217, 7500 AE Enschede, the Netherlands.

E-mail: d.vescovi@utwente.nl; s.luding@utwente.nl

[8pt] ‡ Present address: Department of Civil and Environmental Engineering, Politecnico di Milano, 20133 Milano, Italy. E-mail: dalila.vescovi@polimi.it

1 Introduction

In the recent past, the flow of granular materials has been the subject of many scientific works; this is due to the large number of natural phenomena (i.e., landslides and debris flows) and industrial processes involving solid particles flowing. Granular materials can behave very differently depending on both the micro-mechanical properties of the particles and the macroscopic characteristic of the system (i.e., geometry, gradients of velocity, density). In particular, granular systems, as well as other soft matter such as suspensions of soft particles or polymers, foams, or emulsions, can behave like fluids, meaning that they yield under shear stress, or like solids able to resist applied stresses without deforming.48 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.
As a consequence, a key question concerns how to model the transition from fluid- to solid-like behavior and vice-versa? This transition is reflected in the relation between stresses and deformation rates, that is the rheology. Jammed structures can show a rate independent behavior, but in unjammed granular systems the stresses are proportional to the square of the strain rate (Bagnold scaling1). In the proximity of jamming, a continuous transition between the two extreme regimes takes place, for which the correct rheology is still not fully understood.32
Numerical simulations performed by using perfectly rigid spheres (i.e., infinite stiffness)37, 7 have shown that pressure, shear stress and granular temperature (defined as the mean squared velocity fluctuations, representing a measure of the degree of agitation of the system) diverge when the density approaches jamming, as predicted by the kinetic theory of granular gases.14, 54, 3 As a consequence, jammed flows of rigid particles are not possible at constant volume. In contrast, numerical results with soft particles6, 23, 8 show that steady, constant volume flows are possible also at densities above jamming. Their softness allows the particles to deform and to attain steady shear at very dense, jammed configurations. Experiments performed on concentrated suspension 52 of hard and rigid particles have revealed different physical processes in the straining and the flow behavior of hard and soft spheres suspension, in particular in the unsteady regime, due to the development of permanent contacts of the particles in the case of deformable spheres at large volume fractions.
Although several constitutive models have been proposed in the literature to account for the softness of the particles,8, 4, 31, 38, 2 most of them have some important limitations. In particular, most of them “match” the limits of fluid- and solid-like behavior and many do not provide continuous and differentiable equations for all the variables of the system (i.e., pressure, shear stress and granular temperature). Many models are not able to quantitatively predict steady shear flows at all densities.
For a collection of particles, the jamming transition occurs in the limit of zero confining pressure at the critical volume fraction (where, for a granular material composed of identical particles, the local volume fraction is defined as the ratio of the local material density to the particle density, and the subscript J refers to the jamming transition). Recently, several authors have shown that, under homogeneous steady shearing, stresses and granular temperature can be expressed as asymptotic power-law relations of the shear strain, if scaled by powers of , the distance to jamming.42, 20, 43, 44 However, different authors report different critical exponents in the solid (jammed) and fluid (unjammed) regimes, both numerically42, 20, 8 and theoretically.43 The critical exponents, as well as the jamming volume fraction, are affected by several parameters: the polydispersity of the system, 39, 28 the friction of the particles, the force contact model adopted in the simulations (linear spring or Hertzian) and the ratios of relevant time-scales like the inverse shear rate or square root of the particle stiffness per mass. Even though the jamming density is known to be pressure-, and material- as well as protocol-dependent, see Ref. 32 and references therein, for the purpose of scaling we have to work with a constant . Jamming transition and scaling of quantities with respect to the distance from jamming are also very popular subjects in the fluid-dynamics and statistical-mechanics community. In particular, complex fluids composed of a dispersion of one material in a continuous phase show a transition from mechanically solid-like to fluid-like states, similarly to granular systems, when the shear stress is increased above some critical value, the yield stress.45 Emulsions, colloids, foams, gels and suspensions of (soft) particles or polymers belong to this category. By performing experiments on different kinds of yield-stress fluids, Paredes et al. 45 and Dinkgreve et al. 12 have found that, by appropriate scaling with the distance to jamming, such kind of complex fluid systems allowed a data collapse onto universal power-law relations of the shear strain, in the fluid-like and solid-like regimes. The presence of the “continuum phase”, i.e. the liquid in which droplets or particles are dispersed, leads to critical exponents for the case of complex fluids which differ from those of dry collections of soft particles, especially in the fluid-like state, but the intrinsic rheology is actually surprisingly analogous.
This work focuses on the simple shear flow of an ideal granular material, composed of identical, frictionless spheres, under steady conditions, at constant volume. A series of Discrete Element Method simulations are performed in order to investigate the role of particle stiffness in a wide range of volume fractions, both well below and above jamming. The goal is to propose phenomenological constitutive relations for granular materials merging dilute and dense flow conditions. First, we analyze our numerical results to derive appropriate scaling laws in both unjammed and jammed states, far from the jamming transition (with special functional forms that are needed further on). Then, we use the scaling laws that collapse the data in either solid or fluid regimes to postulate a phenomenological model for the stresses and the granular temperature. In particular, the two regimes are merged in a unique function which is (i) continuous and differentiable at any point and (ii) able to predict the behavior even close to the jamming transition.
This paper is organized as follows: in Section 2 we describe the simulation method and the flow configuration; in Section 3 the difference between rigid and soft particles granular systems are discussed in terms of the rheology15 and the jamming density is evaluated based on the behavior of the coordination number as function of volume fraction and particles stiffness. Section 4 is devoted to the scaling relations obtained from the numerical data; in Section 5 the merged constitutive model is presented and compared with numerical results; in Section 6 we discuss the comparison of the proposed model with the rheologies proposed by Chialvo et al. 8, Berzi and Jenkins 2, Singh et al. 50 and Paredes et al. 45; in particular, the model in Ref. 45 was proposed for emulsion-type systems, and is here revised to adapt to the case of dry assemblies of soft spheres, according to our numerical data. Finally, the results are summarized and concluding remarks are given in Section 7.

2 DEM numerical simulations

Figure 1: Sketch of the constant-volume simple shear flow configuration. A granular material composed of frictionless, soft particles is homogeneously sheared at constant shear rate , with flow in horizontal direction (walls move in opposite horizontal directions with velocity ), by using Lees Edwards boundary conditions. Colors indicate speed, from blue (zero velocity of the particles in the core of the domain) to red (maximum velocity of the particles at the boundaries).

Simple shear flows are homogeneous if the horizontal velocity of the medium is linearly changing along a line perpendicular to the shearing direction and the kinematic variable which affects the problem is its first spatial derivative, the shear rate, , that is then constant along the flow depth. The other variables governing the problem are the volume fraction (or density, or concentration, defined as the fraction of volume occupied by the spheres), the pressure , the shear stress and the granular temperature . The latter is defined as one third of the mean squared particle velocity fluctuations and represents a measure of the degree of agitation of the system, as introduced in the framework of kinetic theories of granular gases. 22, 36, 5, 14, 16
We have performed DEM numerical simulations of steady simple shear flow of frictionless spheres to investigate the role of the particle stiffness , especially at large volume fraction. The simulations and the coarse-graining presented in this paper, were undertaken using the open-source code Mercury-DPM (www.mercurydpm.org).51, 55 The simulations are done under constant volume with a uniform shear in a rectangular box of dimensions , respectively in , and directions (Fig. 1). The shear is applied using Lees-Edwards29 periodic boundary conditions in the direction and periodic boundary conditions are employed in the and directions. In all simulations, we use 2000 particles of diameter , density , i.e., mass , and fix the height of the computational domain as , before we compute the and size accordingly to the chosen volume fraction . The system is sheared applying a constant velocity gradient (shear rate) with mean flow in the direction only. We use a linear spring-dashpot model, then the normal force between particles at contact is computed as , with overlap .34, 30 is the damping coefficient and is related to the normal coefficient of restitution: . In the simulations, we use normal coefficient of restitution , tangential coefficient of restitution , interparticle friction coefficient and normal spring stiffness . The collision time can be analytically obtained as . Simulations have been performed by systematically changing both the volume fraction , ranging from 0.2 to 0.68, and the particle stiffness , such that the dimensionless quantity ranges from to .
At large volume fraction (i.e., ) and large particle stiffness (i.e., ), we observed crystallization of the monodisperse systems, which results in non-homogeneous flows such as shear bands. To avoid crystalline structures, in such situations we use slightly polydisperse systems, with mean particle diameter and uniform polydispersity , where is the ratio of the maximum to the minimum particle diameter. Such a small value of the polydispersity does not affect significantly the measured quantities, and, in particular, does not vary so much the critical volume fraction.28
Here and in the following, we use standard notation to refer to dimensional variables, whereas the star is adopted to denote quantities scaled using the particle diameter , density and stiffness . Then, the scaled granular temperature, pressure, shear stress and shear rate are given, respectively, as: , , and . We also introduce three time scales: the microscopic relaxation time scale associated with the transversal motion of a particle submitted to a pressure : ; the macroscopic time scale associated with the shear rate parallel to the flow: ; and the time scale associated with the particles deformability (contact duration): .

3 Influence of the particle stiffness

A convincing yet simple phenomenological model, which has been used frequently in the literature in the last years, is the rheology.15, 10, 25 According to this model, only three dimensionless variables are relevant for steady shear flows of granular materials: the volume fraction , the stress ratio and the inertial number . The inertial number represents the ratio between the microscopic relaxation and the macroscopic shear time scales: 222Note that some authors define the inertial number using the bulk density instead of the particle density: .. On the planes and , data obtained on different flow configurations collapse in the limit of rigid particles.13 Two main assumptions at the basis of the rheology are: (i) perfectly rigid particles and (ii) homogeneous flow (local rheology). Under these two assumption, and are the only time scales influencing the problem. When the particles are soft (‘softness’ effect8) and/or the flow is not homogeneous (‘non-local’ effects26), the standard laws fail because different mechanisms come into play.
In the case of non-homogeneous flows, boundaries affect the system and flow gradients cannot be neglected anymore. In such situations, Kamrin and Koval 26 have proposed a new “diffusive” state parameter, the granular fluidity, with a diffusive evolution equation, to account for some non-locality (see also Ref. 19, 11).
On the other hand, when the particles are not perfectly rigid, contacts are not instantaneous but take a finite time during which a part of the energy, the elastic potential energy fraction, is stored31 due to the persistent deformations of the particles.
When a steady flow is unjammed and sheared very slowly, rigid and deformable particles exhibit very similar behavior, controlled only by and , because is much smaller than either; when the transition to the solid, jammed phase takes place (and this happens only to soft spheres), the time scale associate to the particles deformability (contact duration) starts to influence the system, and softness effects appear. In this work, we investigate homogeneous steady shear flow of soft particles, but still can disregard non-local effects.
In Fig. 2 we plot the stress ratio (a) and the volume fraction (b) versus the inertial number, measured from our simulations, for different values of the dimensionless particle stiffness , that corresponds to (or, equivalently, ). In particular, we use ranging from to (circles). Also plotted in Fig. 2 are the numerical results of Otsuki and Hayakawa 43, who adopted more rigid particles having from to (squares, here we consider their data obtained using 20000 monodisperse spheres), and those of Mitarai and Nakanishi 37 and Peyneau and Roux 46, which were performed using rigid particles (triangles). Note that in our simulations and in Ref. 37 the coefficient of normal restitution is 0.7, whereas in Ref. 43 and Ref. 46 it is very close to zero. This, however, does not affect the main features of the plots on the and planes, as discussed next.

Figure 2: Stress ratio (a) and volume fraction (b) versus the inertial number, for different values of the dimensionless particle stiffness. Circles, squares and triangles represent, respectively, our simulations, the numerical measurements of Otsuki and Hayakawa 43 and the data for rigid particles obtained by Mitarai and Nakanishi 37 and Peyneau and Roux 46.

As predicted by the rheology, in the case of infinitely hard particles, for vanishing , the stress ratio tends to a constant, asymptotic value, often interpreted as the yield stress ratio. At the same time, the volume fraction saturates to the rigid-limit jamming value (here as detailed in the following). In the case of deformable particles, the data deviate from the rigid limit curves and behave very differently depending on the softness. Data with large stiffness agree with the rigid limit data in a large range of inertial numbers and deviate for very slow shear conditions (i.e., very low due to either small and/or large ). For decreasing and constant stiffness, on the plane, Fig. 2(a), curves do not saturate but continuously decrease; the smaller the stiffness, the larger the when the softness effect arises. The transition to solid overcompressed jammed states ( in Fig. 2b) is connected to this decrease in the stress ratio occurring for soft spheres. Note that the deviation from the standard rigid-particle rheology cannot be due to non-local effects, since all simulations presented are carried out in uniform conditions, i.e. non-homogeneous flow data are disregarded here. Steady flows of very soft particles are possible only at large inertial numbers: for example, for , never reaches .
In Ref. 50 and Ref. 47, the authors observed a similar non-collapse of the stress ratio with the inertial number. They performed steady state flow simulations using soft (slightly frictional) particles in a split-bottom geometry, in which gradients in stresses arise due to gradients in both strain rate and pressure. As a consequence, both non-locality and softness can affect the system and distinguishing between the two effects is not possible. In Ref. 8 homogeneous simple shear flows of deformable, frictional particles have been carried out, and the deviation of the data from the standard rheology was reported only at large inertial numbers, . Conversely, for frictionless spheres, we have observed deviations in the whole range of inertial numbers investigated.

Figure 3: Coordination number versus volume fraction, for different values of the dimensionless particle stiffness.

Fig. 3 shows the mean coordination number (average number of contacts between all particles) in our simulations plotted versus the volume fraction. As observed for frictional spheres,23 at low volume fraction, the coordination number increases with decreasing particle stiffness. Decreasing stiffness is equivalent to an increasing scaled shear rate . This dependency inverts at higher volume fraction. Intuitively, when the system is unjammed the softness makes the contact duration longer resulting in a larger , at a fixed density. On the other hand, when the system is jammed very deformable particles (small stiffness) allow to reach denser configurations than more rigid particles for a fixed value of the average number of contacts.
All the curves intersect each other at , where . In theory, the jamming transition occurs at the isostatic point40, 41, 49 which corresponds to a coordination number equal to 6. In non-sheared isostatic packings, particles that do not belong to the force network (i.e. with exactly zero contacts) do not contribute to the coordination number. In shearing conditions, there may be particles having a finite number of contacts for some short time, which do not contribute to the mechanical stability of the packing. Frictionless particles with less than 4 contacts are considered “rattlers”, since they cannot be mechanically stable and hence do not contribute to the contact network.18, 28 The “corrected” coordination number has been introduced as the ratio between the total number of contacts of the particles with at least 4 contacts, and the rattler fraction is defined as , where is the total number of particles in the system. In Ref. 17, the authors have shown that the corrected and the standard coordination number are related as . In particular, at jamming the corrected coordination number is 6, as in isostatic conditions, and , then the standard coordination number results . This value of is consistent with that at which all the curves intersect and change convexity in Fig. 3, denoting the transition between solid, jammed structures and fluid, unjammed systems. As a consequence, we assume the corresponding value of to be the jamming volume fraction: .
In general, the jamming volume fraction is not a constant (even if pressure or density are held constant),27, 32 but shows a history-, material-, and protocol-dependence. We cannot exclude a priori that also particle stiffness and shear rate affect the jamming density. In our simulations, variability of is small relative to the full range of data investigated, but, in general, it could make a difference, especially for soft particles. Nevertheless, in this work we assume, for simplicity, the jamming volume fraction to be constant and equal to 0.634.

4 Scaling plots

As previously stated, two extreme flow behaviors, commonly interpreted as flow regimes, are associated with different power laws scaling with : stresses scale as in the fluid, unjammed regime, that is below the critical volume fraction (), but show almost no rate dependence in the solid, jammed regime (). At , a continuous transition between the two extreme regimes takes place (jamming transition).
In order to find out the scaling functions for the stresses and the temperature, we introduce the following scaling relation for the generic scaled quantity (with being a place-holder for either , or ):


where , and are constant coefficients, whereas the exponent depends on the flow regime . In this framework, we do not consider the scaling relation to describe the jamming transition. Assuming in the solid and in the fluid regime; far from the jamming volume fraction, Eq. (1) can be expressed as


where , and and are dimensionless constitutive parameters, such that the scaling is encompassed by and the density scaling by .

Next, we determine the critical exponents and the constant coefficients for the granular temperature and the stresses by using the results obtained from our numerical simulations with frictionless spheres. The constant coefficients are very sensitive to both the critical volume fraction and the range of scaled shear rate analyzed.43, 8 In our simulations, the scaled shear rate ranges from to () and, for all the quantities, we find good collapses using , the jamming volume fraction obtained in the previous Section. The extrapolated values of the coefficients for the stresses and the granular temperature resulting from our DEM results are reported in Tab. 4. The values in Tab. 4 are estimated by fitting our numerical data, and are quite similar to those obtained in Ref. 20 for bidisperse mixture of frictionless particles.

6/5 9/5 1 2 0 12/5 6/5
6/5 8/5 1/2 2 1/6 2 14/15
2 3/2 2 2 1 1 1/2
Table 1:  Scaling exponents appearing in Eq. (1) for pressure, shear stress and granular temperature, as inferred from the numerical simulations

Then, the resulting relations for the scaled pressure, shear stress and granular temperature become, respectively:

Figure 4: Collapse of scaled (a) pressure, (b) shear stress and (c) granular temperature plotted against scaled shear rate, for different values of the (dimensionless) particle stiffness as given in the inset. Dash-dotted lines represent the scaling laws in the unjammed regime whereas dashed lines those in the jammed regime. The right most data are those for and are not supposed to collapse on the limit cases.

Figs. 4(a)-(c) show, for the three quantities, the collapse of numerical data obtained using Eqs. (3)-(5), respectively, using . From Figs. 4(a)-(c), we obtain the values of the critical exponents , and in both the fluid (unjammed) and the solid (jammed) regime, by interpolating the collapsing curves at small enough . In particular,

  • when (fluid regime)

  • when (solid regime)

0.634 0.0075 0.60 0.0105 0.12 0.0090 0.05 1.4 0.2
Table 2:  Constant coefficients appearing in Eqs. (6)-(8) as inferred from the numerical simulations

As a consequence, for small , we obtain:


The dimensionless constitutive parameters appearing in Eqs. (6)-(8) can be inferred fitting the power laws in Figs. 4(a)-(c) and are reported in Tab. 2.
We point out that:

  • in both the extreme regimes, , and depend not only on and , but also on some power of the volume fraction (, and , respectively).

  • In the fluid regime, all the three quantities scale quadratically with the scaled shear rate (Bagnold scaling1), but with different (negative) powers of .

  • Differently from what was shown in other works, in the jammed regime only the pressure is rate independent, whereas the shear stress appears to scale with but only in a very limited range. The granular temperature scales linearly with , as already suggested in Ref. 20, 43.

The slight rate-dependency of the shear stress in the solid regime becomes evident when plotting the stress ratio , also known as macroscopic friction, as a function of the volume fraction. As shown in Fig. 5, at small volume fraction, is not much affected by the particle stiffness (i.e. by ): the numerical measurements of obtained using different particle stiffnesses collapse together with the rigid particle data obtained by Mitarai and Nakanishi 37 and Peyneau and Roux 46, deviating only for the most soft particles. This means that pressure and shear stress scale with the same power () of in the unjammed regime. Therefore, the stress ratio in unjammed configurations of homogeneous steadily sheared systems depends only on the volume fraction. For increasing volume fractions, the curves deviate more and more from the rigid particle limit, already below the jamming volume fraction, and remain well separated at . Furthermore, at larger volume fractions, the separation factor is almost the same between all the curves, and equal to , suggesting that in the solid regime. As a consequence, pressure and shear stress must have different asymptotic power-law relations with the scaled shear rate; when is rate independent, turns out to be affected by , and thus by , in the jammed regime.

Figure 5: Stress ratio as a function of the volume fraction for different values of the particle stiffness. Also plotted (black triangles) are the rigid particle data obtained by Mitarai and Nakanishi 37 and Peyneau and Roux 46.

5 Model

The goal of this Section is to propose phenomenological constitutive relations for , and based on the asymptotical analysis discussed in the previous Section. The idea is to merge the equations describing the two scaling regimes with a unique function which must be (i) continuous and differentiable at any point and (ii) able to predict the behaviour even at and around the jamming transition.
In the fluid and solid regimes, pressure, shear stress and granular temperature are given by Eqs. (6)-(8); these equations show the same form, as summarized in Eq. (2). Both branches of Eq. (2) can be solved explicitly for the distance to jamming:


By defining the distance to jamming as the sum of the two right hand side contributions of Eqs. (9) and (10), we obtain a unique continuous and differentiable equation, valid for any value of the volume fraction:


Eq. (11) does not include singularities, but cannot, in general, be solved explicitly for . We can derive implicit functions which relate pressure, shear stress and temperature to the volume fraction, through the scaled shear rate, in the same form as Eq. (11):


It is important to notice that the form of the merging function Eq. (11) is general and can be used for any quantity for which the asymptotical behavior in the fluid and solid regime is known and which is given in the form of Eq. (2). In particular, different values of the critical exponents of can be adopted for , and , as well as different power of the volume fraction.
Moreover, combining the asymptotical equations for and , Eqs. (6)-(7), we obtain the limit equations for the stress ratio :


which, following the same approach as above, results in an implicit, merging function for :


with parameters and .

In the following Eqs. (12)-(14) and (16) are solved and compared with the numerical results of our simple shear simulations.
Fig. 6 (a) and (b) depict the predicted scaled pressure and shear stress, compared with the results of our numerical simulations of simple shearing for a range of scaled shear rates from to . Also shown are the simple shear results of Chialvo and Sundaresan 7 (asterisks) for frictionless spheres having and . The agreement is good for both variables in the whole range of volume fractions investigated and, in particular, in the transitional regime, i.e. around the jamming volume fraction . Only at very small volume fractions () the shear stress is underestimated by the model. In the solid, jammed regime (), the measured pressure data tend to collapse, in agreement with the rate independent behavior observed in other works,24, 20, 43, 8 as predicted by Eq. (6), except for the softest particles (). Conversely, when the contact duration gets into the order of magnitude of the shear time-scale, a slight systematic dependence of the shear stress on the shear rate is observed even at large volume fractions (Fig. 6b), which is well captured by the proposed model.

Figure 6: Scaled (a) pressure and (b) shear stress as functions of the volume fraction for different values of the scaled shear rate (or, equivalently, stiffness: ). Symbols represent the numerical data obtained from numerical simulations and solid lines the proposed model Eqs. (12)-(13). The black asterisks are data from Ref. 7.
Figure 7: Dimensionless (a) pressure and (b) shear stress as functions of the volume fraction for different dimensionless particle stiffness . Symbols represent data from the simulations and solid lines the proposed model Eqs. (12)-(13). The black asterisks are data from Ref. 7; the black triangles are rigid data from Ref. 37, 46 and the dotted lines are Eqs. (17) and (18).

In Fig. 7 (a) and (b), the (dimensional) pressure and the shear rate are made dimensionless using the particle diameter, density and the (dimensional) shear rate . Different colors, as in Fig. 6, indicate different values of the dimensionless particle stiffness , ranging from to . The results of the numerical simulations with rigid particles Ref. 37, 46 are added to the figures, as well as the limit curves for rigid particles (, i.e. , black dotted lines); these equations are identical to those of the unjammed regime in Eqs. (6)-(7) and diverge at the jamming volume fraction:


The predicted curves of and for soft particles collapse on the rigid ones in the fluid, unjammed regime, and start to deviate at different values of the volume fraction depending on the particle stiffness (even for the softest particles the deviation from the rigid limit is smaller than 15% for , for both pressure and shear stress). The softer the particles, the lower the transition density. As observed in Fig. 6(b), the shear stress is strongly underpredicted at very small densities, , where the factor in Eq. (18) strongly affects the expression of , and is too small, not catching the functional behavior in the regime where kinetic theory holds.53 However, the model allows for very good predictions for all the other values of volume fractions.

The model predictions for the granular temperature are compared with the numerical results in Fig. 8 in terms of (a) and (b). Like for the stresses, data with different stiffness collapse on the rigid curve in the fluid regime, at , if scaled with the particle diameter and the shear rate, see Fig. 8(b). According to Eq. (8), the expression of in the rigid case is:


On the other hand, Fig. 8(a) highlights the linear rate dependence (inversely proportional to ) of in the solid regime, where the data do not collapse on a single curve (for large volume fractions). In Fig. 8(b), for the softest particles the deviations from the rigid limit are smaller than for , whereas they are larger than for . Hence, for the granular temperature, conspicuous softness effects arise already at volume fractions smaller than for pressure and shear stress, and are not properly reproduced by Eq. (14).

Figure 8: Scaled and dimensionless granular temperature (a) and (b) as a function of the volume fraction for different values of the dimensionless particle stiffness or, equivalently, scaled shear rate [See legends in Fig. 6 for (a) and in Fig. 7 for (b)]. Symbols represent the data obtained from simulations and solid lines are the proposed model Eq. (14). The black asterisks and triangles are data from Ref. 7 and Ref. 37, respectively; the dotted line is Eq. (19).

Finally, in Fig. 9, we compare the theoretical expression for the stress ratio Eq. (16) with the numerical measurements. The macroscopic friction varies over a narrow range so we can use a linear scale instead of a logaritmic one, differently from what was done for the other quantities, so that the quality of the predictions can be appreciated more accurately. From Eq. (16), the stress ratio in the rigid limit reads:


In Fig. 9(a), the proposed model captures well the stress ratio, especially in the transition and solid regime. Some disagreements still remain in the fluid regime at volume fractions between 0.2 and 0.6, where the model does not capture well the slight -dependence of . Nevertheless, the disagreements are smaller than 10% in the range of volume fractions investigated, see the quality factor in Fig. 9(b) where is the theoretical prediction of given by Eq. (16). At , the model is unable to reproduce the increase of for decreasing , and the stress ratio nullifies when , Eq. (20), in contrast with the data. This discrepancy is due to the bad predictions in the stress ratio at very small volume fractions, and, in particular, due to the inappropriate exponent of the volume fraction in Eq. (18).

Figure 9: (a) Numerical (symbols) and theoretical (lines, Eq. (16)) stress ratio plotted versus volume fraction, for different values of the dimensionless particle stiffness or, equivalently, scaled shear rate. (b) Numerical stress ratio scaled by the proposed prediction , Eq. (16). The disagreements are smaller than 10% in the fluid range of volume fraction investigated, and are considerably better for the transitional and solid regime.

As previously stated, the critical exponents strongly depend on the range of scaled shear rates. Otsuki and Hayakawa 43 performed homogeneous simple shear volume-controlled simulations of monodisperse, frictionless particles having stronger dissipation and using a range of smaller from to . They found exponents different from ours for the quantities, as well as a slightly different value of the jamming volume fraction (0.639 in contrast to our 0.634). Our numerical data thus do not scale with their theoretically predicted exponents, in neither fluid nor solid states. This discrepancy can be due to the different values of the coefficient of restitution adopted, which highlights the sensitivity of scaling results and exponents, in particular, to the magnitude of dissipation. Furthermore, the shear stress measured by Otsuki and Hayakawa exhibits a strictly rate-independent behavior in the solid state, in disagreement with our dependence on , which, however, encompasses both rate and softness effects. We have previously discussed that a rate independent at implies that the stress ratio is only a function of the volume fraction, independent on , in contrast with what is shown in Fig. 5. On the other hand, if we try to predict their numerical results with our model, we do not succeed at very small (or, equivalently, very high stiffness). The choice to relate , and through a power law in the solid, jammed state, while working for our intermediate regime of stiffnesses and shear rates, is probably not the right approach to reproduce flows in a wider range of shear rates, including their very slow and very stiff cases. Nevertheless, the form of the merging function proposed is suitable with any function of the shear rate, provided that it is proportional to some power of the distance to jamming. For example, if in the solid regime the scaled shear stress can be expressed as , with an arbitrary function of and , then the constitutive equation for simply reduces to . The lack of numerical data at larger stiffnesses (and small shear rates) in jammed conditions, due to very long computational time required, prevents us from finding a more appropriate relation for the shear rate in the solid regime, tending to a rate independent behavior in the limit of very large stiffness.

In Fig. 10 we depict, for completeness, the comparisons between the model and the data on the planes (a) and (b), considering curves at constant (that is constant ). The volume fraction can be easily expressed as a function of and by substituting into Eq. (12):


The relation between , and can be obtained combining Eqs. (21) and (16), and involves also the volume fraction. Fig. 10 confirms the good agreement of the model with the data. The discrepancies in the stress ratio at large inertial number, , correspond to those at in Fig. 9(a), while the relation between and density is captured well throughout.

Figure 10: Stress ratio (a) and volume fraction (b) versus the inertial number, for different values of the dimensionless particle stiffness. Symbols represent data from the simulations and solid lines the proposed model.

6 Comparison with other models

In this Section, we compare our model with other theories developed for deformable particles, from the literature. Here we do neither consider rheologies which apply only to rigid particles, such as the model or similar,15, 10, 25, 46 nor standard kinetic theory based models.14, 53, 54 In particular, we focus on the constitutive models proposed by Chialvo et al. 8, Berzi and Jenkins 2 and Singh et al. 50 Moreover, we consider also the model proposed by Paredes et al. 45 derived for emulsion-like systems. In the following, we present the four models and compare them with ours and the numerical data.

  1. Chialvo et al. 8 proposed a new definition of the scaled pressure and the macroscopic friction:



    , and define the scaled pressure in the three regimes: fluid (named inertial in Ref. 8), solid (quasi-static) and intermediate, respectively; is the correction to the standard rheology to account for softness effects; , , , and are dimensionless model parameters. Furthermore, for systems of frictionless, monodispersed particles, the authors in Ref. 8 estimated the jamming volume fraction to be 0.636. The standard rheology is here denoted as and given by


    with and , , for frictionless particles. The same expression for has been adopted by Ness and Sun in Ref. 38, where the authors also add a stress contribution for the intermediate viscous fluid in order to extend the rheology to non-Brownian suspensions.

    Figure 11: Scaled pressure (a) and stress ratio (b) plotted against the volume fraction, for the proposed constitutive model (solid lines, Eqs. (12) and (16)), and the constitutive model of Chialvo et al. 8 (dashed lines, Eqs. (22) and (23)). Symbols represent our numerical data.

    Fig. 11 depicts, for three cases , and , the comparison of the present model (solid lines) with the theory of Chialvo et al. 8 (dashed lines), which underpredicts the pressure and overpredicts the stress ratio in both fluid and solid regimes. Besides the shift, the trend of the stress ratio at small volume fractions, , is qualitatively better captured by the model of Chialvo et al., which predicts the slight rate dependency of even in the fluid regime (Fig. 11b).
    Concerning , the bad prediction of Eq. (22) is mainly due to the powers of which govern the fluid and solid response, which are and , respectively. If we consider only the fluid and solid regimes in the theory of Chiavo et al., and , then we can merge the two equations using our approach:


    The curves obtained using Eq. (25) completely overlap with those of Chialvo et al. in Fig. 11(a) (note that the comparison is not shown in the figure for the sake of clarity) meaning that the two approaches result in very similiar predictions for . The form of Eq. (25) is simpler than Eq. (22), because it (i) does not need the definition of the intermediate regime pressure and (ii) does not require an “if” condition to distinguish below and above the jamming volume fraction, given that any divergence is avoided. Anyway, the critical exponents adopted in our model (-12/5 for the fluid and 6/5 for the solid regime), together with the extra dependence of on and the choice , allow a better quantitative agreement with most of our numerical data. Note however, that we did not fit the parameters of Chialvo et al. to our data, so that the slight shift must be disregarded, while the qualitative behavior can be appreciated.

  2. Berzi and Jenkins 2 have extended the standard kinetic theory to account for the deformability of the particles. Their model can be summarized as follows:


    Similarly to Eq. (22), scaled pressure and shear stress are given by using three contributions: rigid (subscript rig), deformable (subscript def) and elastic (subscript el):

    Here, , , , for the case of frictionless particles having , whereas and are functions of the volume fraction and the coefficient of restitution, derived from kinetic theory, and summarized in Tab. 5. For frictionless particles, the authors have used . The scaled granular temperature is computed as solution to the balance of fluctuation energy and, in simple shear conditions, results in:


    with , , , and .

    Table 3: List of functions in the constitutive relations of Berzi and Jenkins 2
    Figure 12: Scaled pressure (a) and stress ratio (b) plotted against the volume fraction, for the proposed constitutive model (solid lines, Eqs. (12) and (16)), and the constitutive model of Berzi and Jenkins 2 (dash-dotted lines, Eq. (26) and ratio of Eqs. (27) to (26)). Symbols represent our numerical data.

    In Fig. 12 we show the comparison of our model (solid lines) with that of Berzi and Jenkins (dash-dotted lines). For the latter, the stress ratio is computed just dividing Eq. (27) by Eq. (26). The major limitation of the theory in Ref. 2 is that the stress ratio is not affected by the particle stiffness (Fig. 12b). This is due to the assumption that the scaled granular temperature scales with even at , Eq. (28); this assumption, besides being in disagreement with what is shown by the numerical data (Fig. 8a), makes the stress ratio constant in the solid regime:

    given that the model parameters, in simple shearing, are such that . On the other hand, it must be noticed that the theory of Berzi and Jenkins is generally developed for any kind of flow configurations, being described by conservation laws of mass, momentum and fluctuation energy and a full, more complete set of constitutive relations, including also the coefficient of restitution. Moreover, at very low densities, , the model in Ref. 2 corresponds to the standard kinetic theory in Ref. 14 which has been proven to quantitatively predict all the variables, and, in particular the shear stress, differently from the model proposed here.

  3. Singh et al. 50, as rephrased in Ref. 35, 33, have generalized the standard relation accounting for the influence of the particle stiffness through the scaled pressure :


    with given by Eq. (24), and constant parameters , and and . Note that the data in Ref. 50 have been obtained using the ring-shear geometry (inhomogeneous), by local coarse-graining, and slightly frictional and polydisperse particles, which implies a small shift in and some other parameters. In particular, in order to compare our results with the model of Singh et al., here we adopt and suitable for monodisperse systems. Furthermore, the additive corrections in Ref. 50 have been rewritten (identical in first order approximation) as multiplicative correction terms in Ref. 33 that imply an higher order non-linear correction neglected in Ref. 50. Using Eq. (30), can be computed implicitly as function of and (or, equivalently, ). The present theory is compared with the model of Singh et al. in Fig. 13.

    Figure 13: Scaled pressure (a) and stress ratio (b) plotted against the volume fraction, for the proposed constitutive model (solid lines, Eqs. (12) and (16)), and the constitutive model of Singh et al. 50 (dotted lines, Eqs. (30) and (29)). Symbols represent our numerical data.

    As it has been already stated, we have calibrated the parameters in the model of Singh et al. by using our numerical results, as a consequence the model results in a good quantitative agreement with the data, except for very soft particles (). As for the model of Chialvo et al. 8, the stress ratio predicted by the model of Singh et al. is able to qualitatively capture the rate dependency of at .

  4. Finally, Paredes et al. 45 have presented a microscopic two-state theory for yield-stress fluids, elaborated in more detail by Dinkgreve et al. 12, to describe the transition between jammed and unjammed states. Yield-stress fluids are complex fluids composed of dispersions of one material (particles, drops or bubbles) in a liquid (or continuous phase), whose mechanical behavior is characterized by the emergence of a yield stress for volume fractions higher than some critical value (the jamming volume fraction) and Newtonian flow for lower volume fractions, with shear thinning in either case for high shear rates. The measurements of the shear stress obtained in the experiments of Paredes et al. 45 and Dinkgreve et al. 12 with different kinds of yield-stress fluids, collapse when rescaling all data as: - , with and scaling parameters, similar to our Fig. 4(b). According to this collapse, the scaled shear stress has been found to obey: (i) the Herschel-Bulkley 21 equation when , with yield stress expressed as a power law in the distance to jamming, and (ii) the Cross equation 9 when , where the Newtonian viscosity satisfies a power law in . This result can be summarized as


    where , and are adjustable dimensionless parameters and the coefficients , and are given in Fig. 14. The yield stress and the Newtonian viscosity are given by and , respectively.
    As previously stated, the critical exponents, as well as the jamming volume fraction, are material-dependent. The original model of Paredes et al. has been derived for soft matter systems, and although it is compatible in spirit with the granular rheology, is completely different from our granular fluid in all quantitative numbers and exponents, especially below jamming. As a consequence, it requires careful calibration of the parameters, differently from the models previously analyzed, specifically derived for granular materials. In particular, the exponent , appearing in the unjammed phase, is measured equal to 1 for any kind of yield-stress fluid considered. implies that, in the fluid regime, far from jamming, the shear stress varies linearly with the shear rate. This linear dependence, typical of viscous liquids, does not apply to dry granular systems, where (Bagnold scaling), implying . In order to adapt the model of Paredes et al. to dry granular systems, we use in Eq. (31) and calibrate the critical exponents and the dimensionless parameters of the model by collapsing our numerical data on the plane - , Fig. 14. We obtain , and (whereas Dinkgreve et al. 12 have estimated, for different kinds of yield-stress fluids, , and ). These values coincide with those reported in Tab. 4 for , in fact Fig. 14 corresponds to Fig. 4(b) unless for (actually this last provides a better collapse of the data in the unjammed regime). In Fig. 14, Eq. (31) is also plotted: the dashed line represents the Herschel-Bulkley (branch at ) and dash-dotted line represents the Cross equation () with . By fitting the data, we infer the dimensionless parameters appearing in Eq. (31): , and . Note that both branches of Eq. (31) are defined in proximity of jamming and fit well the data in this region. Moreover, both the Herschel-Bulkley equation and the Cross equation lead to the same expression at the jamming transition .

    Figure 14: (a) Collapse of scaled shear stress plotted against scaled shear rate as suggested in Ref. 45, for different values of the (dimensionless) particle stiffness. Dash-dotted line () and dashed line () represent the two branches of Eq. (31) where , , , , and . (b) Scaled shear stress plotted against the volume fraction, for the proposed constitutive model (solid lines, Eq. (13)), and the model of Paredes et al. 45 (dash-dash-dotted lines, Eq. (31)). Symbols represent our numerical data. The inset shows the scaled shear stress predicted by Eq. (31) for the case