Abstract
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 timescale 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 stiffnessdependency.
Merging fluid and solid granular behavior
[0.4cm]
Dalila Vescovi and Stefan Luding
[0.4cm]
Multi Scale Mechanics (MSM), CTW, MESA+, University of Twente, PO Box 217, 7500 AE Enschede, the Netherlands.
Email: d.vescovi@utwente.nl; s.luding@utwente.nl
[8pt] ‡ Present address: Department of Civil and Environmental Engineering, Politecnico di Milano, 20133 Milano, Italy. Email: 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 micromechanical 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 solidlike 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 solidlike behavior and viceversa? 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 scaling^{1}). 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 particles^{6, 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 solidlike 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 powerlaw 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 numerically^{42, 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 timescales 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 protocoldependent, 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 fluiddynamics and statisticalmechanics community. In particular, complex fluids composed of a dispersion of one material in a continuous phase show a transition from mechanically solidlike to fluidlike 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 yieldstress 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 powerlaw relations of the shear strain, in the fluidlike and solidlike 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 fluidlike 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 rheology^{15} 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 emulsiontype 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
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 coarsegraining presented in this paper, were undertaken using the opensource code MercuryDPM (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 LeesEdwards^{29} 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 springdashpot 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 nonhomogeneous 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: ^{2}^{2}2Note 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’ effect^{8}) and/or the flow is not homogeneous (‘nonlocal’ effects^{26}), the standard laws fail because different mechanisms come into play.
In the case of nonhomogeneous 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 nonlocality (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 stored^{31} 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 nonlocal 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.
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 rigidlimit 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 rigidparticle rheology cannot be due to nonlocal effects, since all simulations presented are carried out in uniform conditions, i.e. nonhomogeneous 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 noncollapse of the stress ratio with the inertial number. They performed steady state flow simulations using soft (slightly frictional) particles in a splitbottom geometry, in which gradients in stresses arise due to gradients in both strain rate and pressure. As a consequence, both nonlocality 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.
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 point^{40, 41, 49} which corresponds to a coordination number equal to 6.
In nonsheared 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 protocoldependence. 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 placeholder for either , or ):
(1) 
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
(2) 
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 
Then, the resulting relations for the scaled pressure, shear stress and granular temperature become, respectively:
(3)  
(4)  
(5) 
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 
As a consequence, for small , we obtain:
(6)  
(7)  
(8) 
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 scaling^{1}), 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 ratedependency 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 powerlaw relations with the scaled shear rate; when is rate independent, turns out to be affected by , and thus by , in the jammed regime.
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:
(9)  
(10) 
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:
(11) 
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):
(12)  
(13)  
(14) 
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 :
(15) 
which, following the same approach as above, results in an implicit, merging function for :
(16) 
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 timescale, 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.
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:
(17)  
(18) 
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:
(19) 
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).
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:
(20) 
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).
As previously stated, the critical exponents strongly depend on the range of scaled shear rates. Otsuki and Hayakawa ^{43} performed homogeneous simple shear volumecontrolled 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 rateindependent 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):
(21) 
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.
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 emulsionlike systems. In the following, we present the four models and compare them with ours and the numerical data.

Chialvo et al. ^{8} proposed a new definition of the scaled pressure and the macroscopic friction:
(22) (23) where:
, and define the scaled pressure in the three regimes: fluid (named inertial in Ref. ^{8}), solid (quasistatic) 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
(24) 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 nonBrownian suspensions.
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:(25) 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.

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:
(26) (27) 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:
(28) with , , , and .
Table 3: List of functions in the constitutive relations of Berzi and Jenkins ^{2} In Fig. 12 we show the comparison of our model (solid lines) with that of Berzi and Jenkins (dashdotted 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.

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 :
(29) (30) with given by Eq. (24), and constant parameters , and and . Note that the data in Ref. ^{50} have been obtained using the ringshear geometry (inhomogeneous), by local coarsegraining, 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 nonlinear 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.
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 .

Finally, Paredes et al. ^{45} have presented a microscopic twostate theory for yieldstress fluids, elaborated in more detail by Dinkgreve et al. ^{12}, to describe the transition between jammed and unjammed states. Yieldstress 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 yieldstress 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 HerschelBulkley ^{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
(31) 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 materialdependent. 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 yieldstress 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 yieldstress 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 HerschelBulkley (branch at ) and dashdotted 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 HerschelBulkley equation and the Cross equation lead to the same expression at the jamming transition .