Extended kinetic theory applied to pressure-controlled shear flows of frictionless spheres between rigid, bumpy planes

Abstract

We numerically investigate, through discrete element simulations, the steady flow of identical, frictionless spheres sheared between two parallel, bumpy planes in the absence of gravity and under a fixed normal load. We measure the spatial distributions of solid volume fraction, mean velocity, intensity of agitation and stresses, and confirm previous results on the validity of the equation of state and the viscosity predicted by the kinetic theory of inelastic granular gases. We also directly measure the spatial distributions of the diffusivity and the rate of collisional dissipation of the fluctuation kinetic energy, and successfully test the associated constitutive relations of the extended kinetic theory, i.e., a kinetic theory which includes the role of velocity correlations. We then phrase and numerically integrate a system of differential equations governing the flow, with suitably modified boundary conditions. We show a remarkable qualitative and quantitative agreement with the results of the discrete simulations. In particular, we study the effect of (i) the coefficient of collisional restitution, (ii) the imposed load and (iii) the bumpiness of the planes on the profiles of the hydrodynamic fields, the ratio of shear stress-to-pressure and the gap between the bumpy planes. Finally, we predict the critical value of the imposed load above which crystallization occurs, based on the value of the solid volume fraction near the boundaries obtained from the numerical solution of the kinetic theory. This notably reproduces what we observe in the discrete simulations.

---

1 Introduction

In most industrial and geophysical applications involving granular fluids, e.g., hopper discharge, transport on conveyor belts, landslides, the presence of solid boundaries induces strong spatial heterogeneities and controls the flow.

The discrete nature of granular media suggests the adoption of particle-based numerical methods (discrete element methods, DEM).¹ However, the high computational costs for solving realistic boundary-valued problems lead to consider continuum approaches as the best suited for dealing with large granular assemblies. Nevertheless, particle simulations may provide physical insights into both the microscopic and the macroscopic dynamics of the granular system, and numerical data can be used to test the continuum models. The latter require appropriate boundary conditions to solve boundary-valued problems and provide the spatial distribution of the associated hydrodynamic fields.

Among the various continuum models proposed in the literature for rigid particles, the most widely adopted is the μ − I rheology,^2,3 which relates phenomenologically and algebraically the stress ratio (μ) – the particle shear stress made dimensionless with the particle pressure – to the inertial number (I) – the dimensionless particle shear rate. Extensions to address nonlocal effects arising from interactions with boundaries have been proposed.⁴ For soft particles, continuum models have also been introduced to account for their finite stiffness.^5–9 It has been recently shown¹⁰ that μ − I models can be considered as special cases of the kinetic theory of granular gases,^11–14 extended to account for correlations in the velocity fluctuations.¹⁵ The main consequence of the development of velocity correlations is that the rate of collisional dissipation of the fluctuation kinetic energy is over-predicted by classical kinetic theories. This is because the single particle velocity distribution begins to differ from the distribution of the relative velocity between two neighbours,¹⁵ and the measure of the single particle agitation is not anymore indicative of the frequency of collisions. The introduction of a correlation length^16–18 in the expression of the rate of dissipation permits the reproduction of the dependence of the granular temperature (one-third of the mean square of the particle velocity fluctuations) on the solid volume fraction in homogeneous shearing flows of rigid spheres.¹⁹ This extended kinetic theory can successfully reproduce measurements of solid volume fraction, mean velocity, granular temperature and stresses in heterogeneous granular systems, such as gravity driven flows,^20,21 flows through vertical channels²² and Couette flows,²³ once appropriate boundary conditions are adopted.

In the case of flows over rigid, bumpy planes, boundary conditions for the rate of supply of momentum and fluctuation kinetic energy to nearly elastic, frictionless spheres have been derived by Richman.²⁴ Empirical modifications for the rate of supply of momentum in the case of inelastic, frictionless spheres have been proposed,^23,25 by using the results of DEM simulations of volume-imposed Couette granular flows.

In this work, we first perform DEM simulations of steady, pressure-imposed, heterogeneous flows of identical, frictionless, inelastic spheres sheared between rigid, parallel bumpy planes, in the absence of gravity. We employ a value of the particle stiffness large enough to not affect the results of our simulations. We investigate the roles of the particle inelasticity, the imposed pressure and the bumpiness of the planes, where the first is a material property, while the latter two are associated with boundary conditions. We then use the measurements of stresses, solid volume fraction, granular temperature and velocity gradient to assess the validity of the constitutive relations for the particle pressure and viscosity of the kinetic theory, as previously done in different configurations.^20,23 For the first time in 3D systems, to our knowledge, we also test directly the constitutive relations for the flux of fluctuation kinetic energy and for the rate of collisional dissipation. Measurements of the same quantities were conducted in 2D molecular dynamics simulations of inelastic disks.^26–29

Then, we check the validity of the boundary condition for the rate of supply of momentum proposed by Berzi and Vescovi²⁵ and suggest a simple modification to the rate of supply of fluctuation kinetic energy derived by Richman²⁴ to account for the increase in the energy dissipation for inelastic particles.

We phrase the set of differential equations governing the steady, pressure-imposed shearing flow, using the balance of fluctuation kinetic energy and the constitutive relations of extended kinetic theory, that we numerically solve using the aforementioned boundary conditions, given the total amount of particles in the system. The quantitative agreement between the measurements in the DEM simulations and the predictions of extended kinetic theory are remarkable, in terms of profiles of solid volume fraction, granular temperature, mean velocity and fluctuation energy flux. The theory is also able to notably predict the distance between the moving boundaries and the macroscopic friction (ratio of the particle shear stress to the particle pressure), which represents a measure of the shear resistance, and is a key quantity to test the capability of the granular material to act (or not) as a lubricant between the moving boundaries.

Finally, we observe in the DEM simulations the transition from a random, fluid state to an ordered, crystalline state when the imposed pressure exceeds a critical value. The extended kinetic theory does not apply in those situations. Nevertheless, by assuming that crystallization occurs when the solid volume fraction at the bumpy boundaries reaches the freezing point, we can fairly predict the dependence of this phase-transition on the particle inelasticity. Although not tested against DEM simulations, we also discuss the dependence of the phase-transition on the bumpiness of the boundaries and the amount of particles in the system.

2 DEM simulations

We study the motion of a collection of identical, frictionless spheres, of diameter d and mass density ρ_p, sheared between two parallel bumpy planes, in the absence of gravity (Fig. 1(a)). The two bumpy boundaries move at constant relative velocity V in the direction tangent to the planes themselves. We take x and z to be the flow and shearing directions, respectively, and neglect variations in the transversal direction y. The planes are made bumpy by gluing spheres of diameter d_w in a regular hexagonal fashion, at close contact and aligned with the direction of the flow, as shown in Fig. 1(b). We take z = 0 to be at a distance of d/2 from the edges of the bottom boundary particles. The top of the flowing particles is at z = H, at a distance of d/2 from the edges of the top boundary particles. This choice is motivated by previous studies which looked at boundary conditions for granular flows over rigid bumpy bases.²⁴

Fig. 1 (a) Snapshot of a typical pressure-imposed simulation once the steady state is attained. Also shown is the frame of reference. (b) Top view of the bumpy plane composed of regularly arranged spheres. (c) Temporal evolution of the instantaneous (light grey line) and the smoothed (black line) stress ratio for a typical simulation in steady state (p = 2.9 × 10⁻⁵, d_w = 1 and e_n = 0.5).

We use the open-source code LAMMPS† (ref. 30) to perform discrete element simulations on a collection of N = 3150 spheres interacting via linear spring-dashpots. Periodic boundary conditions are employed in the x- and y-directions. The size of the simulated domain is L_x = 20 and particle diameters along x and y, respectively. The mass of the flowing particles per unit basal area of the rigid boundary is the mass hold-up, h = ρ_p(π/6)d³N/(L_xL_y). The boundary particles may have a different diameter from the moving particles. The bumpiness of the boundaries can be measured through the penetration angle ψ, defined as sin ψ = (d_w + l)/(d_w + d),²⁴ where l is the distance between the edges of two adjacent boundary spheres. In this work, as mentioned, we set l = 0 and vary d_w in order to investigate the effect of the bumpiness.

In the simulations, all the variables are made dimensionless using particle mass density and diameter and the relative velocity between the boundaries: then e.g., lengths, velocities and stresses are given in units of d, V and ρ_pV², respectively. In the case of frictionless spheres, the collision does not alter the relative velocity in the plane tangent to the point of contact. The particles are characterized by the stiffness k of the linear spring, which is set equal to 2 × 10⁵ for both moving and boundary particles. The collision between two spheres is completely characterized by the coefficient of normal restitution, e_n, that is the negative ratio of the relative velocities of the two colliding particles after and before the collision along the direction normal to the point of contact. In each simulation, we set the same value of e_n for both moving and boundary particles.

The integration time step is set equal to t_c/50, where t_c = min(t^ij_c) and t^ij_c is the duration of the collision between particle i and particle j which for the linear spring-dashpot contact model is

where m_ij = πd_i³d_j³/[6(d_i³ + d_j³)], being d_i and d_j the diameter of particle i and j, respectively, equal either to d for moving particles or d_w for boundary particles. We save data every 5000 simulation time steps, that is the saving time step of the simulations is set equal to 100t_c.

For a given total number of flowing particles in the system, the simulations can be carried out by either (i) imposing the relative distance between the moving planes and measuring the normal force acting on them (volume-imposed simulations) or (ii) imposing the normal force per unit area – the pressure p, if we neglect differences in the normal stresses – acting on the planes and measuring the distance between them (pressure-imposed simulations).

Unlike previous works which consider volume-imposed shear flows,^{23,25,31–33} here we focus on pressure-imposed conditions, which seem more realistic in view of practical applications of granular materials as lubricants. However, we have performed a large number of discrete simulations and checked that, in the steady state, volume-imposed and pressure-imposed simulations would give exactly the same results, if the gap between the planes and the normal force on them are matched. We point out that, in the steady state, no dilation occurs even in pressure-imposed situations. In order to systematically analyse the role of the coefficient of restitution, the imposed pressure and the bumpiness, we have carried out 25 simulations, with e_n = {0.2, 0.4, 0.5, 0.7, 0.8}, p = {2.9 × 10⁻⁵, 2.9 × 10⁻⁴, 2.9 × 10⁻³} and d_w = {1/4, 1/2, 1}.

In general, DEM simulations provide the complete microscopic description of the system at each time step, i.e., position and velocity of each particle, as well as inter-particle forces at contact. Macroscopic, hydrodynamic fields can be determined by appropriate averaging,³⁴ here consisting of both spatial and temporal averaging.

The continuum fields, as described in the next section, are: the x-component of the mean velocity, u; the solid volume fraction, ν; the pressure, p; the shear stress, s; the granular temperature, T (one third of the mean square of the particle velocity fluctuations); the rate of collisional dissipation of the fluctuation kinetic energy (the kinetic energy associated with the fluctuations around the particle mean velocity), Γ; and the z-component of the flux of fluctuation kinetic energy, Q.

To obtain the local, continuum fields from discrete measurements, we adopt the averaging method described in ref. 35 and summarized in Appendix A, in which the coarse-graining (weighting) function^36,37 is given by the Dirac delta function. We have checked that using a Gaussian rather than a Dirac delta weighting function does not affect the coarse-graining.

In the absence of gravity, the shear stress and the pressure are spatially uniform. The periodic boundary conditions along x and y ensure that, in the steady state, the remaining variables of the problem can only change along the z-direction. We uniformly divide the domain in slices, parallel to the x–y plane, of thickness equal to one diameter in the z-direction, and we coarse-grain within each slice to obtain local profiles of the continuum variables.

As observed in previous works,^38,39 the measured inter-particle contact forces show large temporal fluctuations, even in the steady state, as contacts are strongly intermittent and the system is still relatively small. This is reflected in corresponding large fluctuations of stresses, rate of collisional dissipation and fluctuation energy flux. As an example, Fig. 1(c) shows the time evolution of the instantaneous stress ratio, s/p, in a simulation with p = 2.9 × 10⁻⁵, d_w = 1 and e_n = 0.5, in steady state.

The instantaneous continuum fields are smoothed using a moving-average filter. We set the size of the moving-average window equal to 5000 saving time steps, that is the minimum above which the average is independent of the size of the window itself. The smoothed data collected in the steady state, over a total of 20 000 saving time steps are then employed to obtain time-averages and standard deviations. This permits to add error bars to the measurement points in the following plots.

3 Continuum model

As mentioned, in the absence of gravity, the momentum balances along x and z imply that, in the steady state, the pressure and the shear stress are uniform in the domain. The balance of fluctuation kinetic energy reduces to a balance between diffusion (due to collisions), production (through the work of the shear stress) and dissipation (due to the inelasticity of contacts)²³

Q′ = su′ − Γ. (1)

Here and in what follows, a prime indicates the derivative with respect to the z-direction. Eqn (1) is the only balance equation governing the problem, and must be supplied with the constitutive relations for p, s, Γ and Q, and appropriate boundary conditions.

3.1 Constitutive relations

We adopt the constitutive relations of the kinetic theory of inelastic granular gases,¹³ modified to account for the development of correlations in the velocity fluctuations at solid volume fractions larger than 0.49^15,17,19 (extended kinetic theory):

p = [1 + 2(1 + e_n)νg₀]νT; (2)

Here, the coefficients J and M are explicit functions of the coefficient of normal restitution, the solid volume fraction and the radial distribution function at contact, g₀, and are reported in Table 1. Such expressions were derived by Garzó and Dufty¹³ by generalizing to inelastic particles the standard kinetic theory of nearly elastic granular gases. The latter would underestimate s and Q by 30% and 90%, respectively, in the case of e_n = 0.5. We employ the expression of the function g₀, which accounts for excluded volume and hindering effects and depends on the volume fraction, proposed in ref. 23 and reported in Table 1. This expression combines the well known formulae of Carnahan and Starling⁴⁰ and Torquato:⁴¹ like the latter, the radial distribution function of Table 1 diverges when ν approaches the random close packing, ν_rcp = 0.636. Unlike the Torquato's expression, however, its derivative with respect to ν is continuous over the entire range of admissible values of solid volume fraction, thus facilitating the numerical integration of the differential equations. Other expressions for the radial distribution function at contact in both 2D⁴² and 3D⁴³ can be found in the literature.

Table 1 List of auxiliary functions in the constitutive relations of kinetic theory, in which Θ(·) is the Heaviside function

---

In eqn (4), L is the correlation length of extended kinetic theory,^16,17 a phenomenological quantity which has been introduced to reduce the collisional rate of dissipation whenever correlations in velocity fluctuations are present. Its expression, obtained from discrete numerical simulations of homogeneous shearing, reads:¹⁹

We point out that the only difference between standard kinetic theory (i.e., kinetic theory based on the Enskog equation and generalized to inelastic collisions, as summarized in ref. 13) and the extended kinetic theory here adopted (originally developed in ref. 16 and 17), lies in the definition of L into the constitutive relation of the rate of collisional dissipation of the fluctuation kinetic energy, eqn (4). In particular, L = 1 in standard kinetic theory, whereas it is given by eqn (6) in its extended form.

The expression for the flux of fluctuation kinetic energy, eqn (5), should also contain a term proportional to the gradient of the solid volume fraction.¹³ However, this term is negligible in dense flows.⁴⁴ We will show later that eqn (5) permits the satisfactory reproduction of the measurements in our DEM simulations even in dilute conditions.

We now test the constitutive relations of eqn (2)–(5) against the DEM simulations described in the previous section. In Fig. 2 and 3, we show the local measurements performed in five DEM simulations, taken at different z-location in the domain, with e_n = 0.5 and different values of the imposed pressure and the bumpiness. Similar agreement, not shown here for brevity, is obtained for all the simulations.

Fig. 2 Measurements in DEM simulations (symbols) against predictions of kinetic theory (eqn (2) and (3), lines) of (a) scaled pressure and (b) scaled shear viscosity as functions of the solid volume fraction, for e_n = 0.5.

Fig. 3 Measurements in DEM simulations (symbols) against predictions of kinetic theory (eqn (4) and (5), lines) of (a) the scaled rate of collisional dissipation and (b) the scaled diffusivity as functions of the solid volume fraction, for e_n = 0.5. The solid and dashed lines in (a) are eqn (4) when the correlation length L is given by eqn (7) and when L = 1, respectively.

Fig. 2(a) and (b) confirm that the kinetic theory notably reproduces the behaviour of the measured scaled pressure, p/T, and shear stress (shear viscosity), s/(T^1/2u′), over the entire range of solid volume fraction less than the random close packing, even in pressure-imposed heterogeneous shearing flows. The constitutive relations for the pressure and the shear stress have been already tested in previous works on homogeneous shearing flows,⁴⁵ volume-imposed heterogeneous shearing flows²³ and inclined flows over rigid, bumpy beds.²⁰

Having directly measured in the DEM simulations the rate of collisional dissipation and the flux of fluctuation energy, we are able to test, for the first time to our knowledge, their corresponding constitutive relations. Fig. 3(a) depicts the scaled rate of collisional dissipation, Γ/T^3/2, as a function of the volume fraction. Once again, the measurements are well reproduced by the expression of kinetic theory, eqn (4). At volume fractions larger than 0.49, with L = 1, the standard kinetic theory would overestimate the measured rate of collisional dissipation (dashed line in Fig. 3(a)). If L of eqn (6) is accounted for, the predictions of the kinetic theory are well within the error bars of the measurements (solid line Fig. 3(a)). For simplicity, in Fig. 3(a), we employ the expression of the correlation length valid in the case of steady, homogeneous shearing flows,⁴⁵ that is when Q′ is neglected in eqn (1):

Finally, Fig. 3(b) shows the dependence of the negative of the scaled flux of fluctuation energy (the scaled pseudo-thermal diffusivity), −Q/(T^1/2T′), on the solid volume fraction, as measured in the DEM simulations and predicted by eqn (5). The agreement is satisfactory, and suggests that the flux of fluctuation energy mainly depends on the temperature gradient even at solid volume fractions as small as 0.05.

3.2 System of differential equations and boundary conditions

Using eqn (1)–(5), and introducing the partial hold-up, (the particle mass per unit area of the boundary contained below z), we phrase the following set of first-order differential equations:

m′ = ν. (11)

In eqn (9), ∂g₀/∂ν is the partial derivative of the radial distribution function at contact reported in Table 1 with respect to the solid volume fraction. The correlation length in eqn (10) is that of eqn (6), with u′ determined from eqn (8). The granular temperature is determined algebraically by inverting eqn (2), as T = p/[ν + 2(1 + e_n)ν²g₀].

Boundary conditions for collisional flows over rigid, bumpy planes can be expressed in terms of the slip velocity, u_w, and the flux of fluctuation energy, Q_w, injected into the flow.²⁴ Then, the set of differential equations, eqn (8)–(11) is subjected to the following six boundary conditions:

m(z = 0) = 0; (12)

Q(z = 0) = Q_w; (14)

m(z = H) = h; (15)

Q(z = H) = −Q_w. (17)

The two additional boundary conditions permit the determination of the shear stress s and the height H as parts of the solution.

Using statistical mechanics arguments, Richman²⁴ derived the flux of fluctuation energy from a rigid, bumpy base in the case of nearly elastic, frictionless spheres as

where the first term on the right hand side is the energy produced by the work of the shear stress through the slip velocity, and the second term is the energy dissipated in collisions. There, T_w is the value of the granular temperature at the boundary. Richman²⁴ obtained B = 1 for nearly elastic spheres, but we have found that B = 3 permits to better reproduce the results of the DEM simulations, probably because we employ highly dissipative particles. By inverting eqn (18), and using the measurements of s, u_w, Q_w and T_w from the present DEM simulations permit to plot the coefficient B as a function of the coefficient of restitution, e_n (Fig. 4(a)). Unfortunately, there is no clear trend, likely because of the difficulties in measuring with sufficient accuracy the energy flux near the boundaries. However, B = 3 approximates the average of the data inferred from the simulations.

Fig. 4 Coefficients (a) B as a function of the coefficient of restitution and (b) C as a function of the scaled slip velocity, measured in DEM simulations of pressure-imposed (open squares) and volume-imposed (filled circles, after ref. 25) shearing flows between rigid, bumpy planes. Also shown are the expressions B = 3 and eqn (20) (solid lines).

For the slip velocity, a modified version of the original expression of Richman²⁴ was proposed by Jenkins et al.,⁴⁶

where C is a function of the scaled slip velocity, ψu_w/T_w^1/2. Fitting against DEM simulations of volume-imposed, heterogeneous shearing flows²⁵ suggested a power-law dependence, with an exponent equal to 3/2. This scaling is confirmed in Fig. 4(b), in which we plot the coefficient C, obtained by inverting eqn (19) and using the measurements of s, u_w and T_w in the present DEM simulations, as a function of the measured scaled slip velocity. The data can be fitted by the following expression:We point out that eqn (20) is different from the expression reported in ref. 25, because the scaled slip velocity was there inadvertently multiplied by the distance H. With this error corrected, the measurements in volume-imposed and pressure-imposed shearing flows between rigid, bumpy planes are in agreement (Fig. 4(b)).

We solve numerically the system of differential eqn (8)–(11) with the boundary conditions (12)–(17), and eqn (18)–(20), using the ‘bvp4c’ routine in Matlab. The inputs are the imposed values of p, d_w (or alternatively the bumpiness ψ), e_n and the mass hold-up h. The outputs are the profiles of u, ν, Q, and m, the shear stress s and the distance H. From the pressure and the profile of solid volume fraction, the profile of the granular temperature can also be determined from eqn (2).

4 Results and comparisons

In this section, the predictions provided by the theoretical model are tested against the results of the present DEM simulations. We consider separately the effect of (i) the coefficient of restitution, (ii) the imposed pressure and (iii) the bumpiness, on the profiles of mean velocity, solid volume fraction, granular temperature and fluctuation energy flux.

In Fig. 5, we compare the results of the simulations and the model predictions when the restitution coefficient changes from 0.5 to 0.8, with p = 2.9 × 10⁻⁴ and d_w = 1. The predictions (lines) are in excellent qualitative and quantitative agreement with the simulations (symbols), for all values of e_n. The presence of the bumpy planes induces a strong heterogeneity, with significant spatial gradients in both solid volume fraction and granular temperature. For all values of the coefficient of restitution, the particles form a dense core region, surrounded by two more dilute layers, -shear bands- (Fig. 5(b)). This is due to the presence of the bumpy planes, where we observe large slip velocity (Fig. 5(a)) and granular temperature (Fig. 5(b)). Then, the bumpy planes act as local sources of fluctuation energy which originate the two “hot” (high T) and low density zones, where most of the shear rate is localized (Fig. 5(a)), sandwiching the “cold” and dense quiescent core, in which the energy flux vanishes (Fig. 5(d)). This is very much reminiscent of the Leidenfrost effect in molecular fluids.⁴⁷ We emphasize that shear bands have been previously observed in 2D discrete simulations of shearing flows even in the absence of physical boundaries.³³

Fig. 5 Measured (symbols) and predicted (lines) profiles of (a) mean velocity, (b) solid volume fraction, (c) granular temperature and (d) normalized fluctuation energy flux, Q/p, for p = 2.9 × 10⁻⁴, d_w = 1, h = 7.94, and different values of the coefficient of restitution e_n (see the legend in (a)).

As already observed in volume-imposed, heterogeneous shearing flows,⁴⁵ the size of the dense core and the temperature-difference with the shear zones increase as the coefficient of restitution decreases. At the smallest restitution coefficient, e_n = 0.5, the dense core occupies half of the domain and it is three orders of magnitude colder than the boundaries.

The capability of the kinetic theory to reproduce the presence of the shear bands, at least in case of volume-imposed bounded shear flows, was already shown in other works by means of linear stability, bifurcation and nonlinear analysis.^48–53 Nevertheless, in this work we want to probe the ability of the theory to provide predictions which are quantitatively reliable. In this perspective, the extension of the kinetic theory to account for the correlation length plays an important role, making the predictions much more accurate. As an example, Fig. 6 shows that, if the velocity correlation is not taken into account (L = 1), the solid volume fraction in the dense core region is overestimated by as much as 7%, whereas the granular temperature is underestimated by one order of magnitude.

Fig. 6 Profiles of (a) solid volume fraction and (b) granular temperature predicted by the theory when the correlation length is given by eqn (6) (solid lines) and L = 1 (dashed lines), for the case p = 2.9 × 10⁻⁴, d_w = 1 and e_n = 0.7. Also plotted are the numerical measurements obtained in the simulation (symbols).

In Fig. 7, we show the effect of varying the imposed pressure over two orders of magnitude when e_n = 0.5 and d_w = 1. The pressure has little influence on the profiles of mean velocity (Fig. 7(a)) and normalized flux of fluctuation energy (Fig. 7(d)). Perhaps intuitively, the size of the dense core increases with the pressure, while both the solid volume fraction in the dense core (Fig. 7(b)) and the granular temperature at the boundaries (Fig. 7(c)) are independent of p. Conversely, given the relation between ν and T through the equation of state (eqn (2)), the granular temperature in the dense core and the solid volume fraction at the boundaries strongly increase when the pressure increases. As the pressure increases, the distinction between the dense core and the shear zones fades, and the granular system is much less heterogeneous. We anticipate that the increase of the solid volume fraction at the boundaries with the imposed pressure can explain the transition from a random to a crystalline state that we observe in the DEM simulations for sufficiently high values of p. Once again, the continuum model can remarkably reproduce the results of the DEM simulations.

Fig. 7 Measured (symbols) and predicted (lines) profiles of (a) mean velocity, (b) solid volume fraction, (c) granular temperature and (d) normalized fluctuation energy flux, Q/p, for e_n = 0.5, d_w = 1, h = 7.94, and different values of the imposed pressure p (see the legend in (a)).

Finally, in Fig. 8 we investigate the role of the bumpiness by changing the diameter of the particles which constitute the rigid boundaries, d_w, from 1 to 1/4, that is, ψ varies between 0.52 and 0.20, with e_n = 0.5 and p = 2.9 × 10⁻⁵. The reduction of the bumpiness yields a larger slip velocity at the boundaries (Fig. 8(a)), and a wider (Fig. 8(b)) and colder (Fig. 8(c)) dense core region. This is due to the fact that boundaries that are less bumpy are less effective in transferring particle momentum from the flow, x-, to the gradient, z-direction, with subsequent reduction in the flux of fluctuation energy from the boundaries to the flow (Fig. 8(d)). All the profiles are reproduced with high accuracy by the continuum model.

Fig. 8 Measured (symbols) and predicted (lines) profiles of (a) mean velocity, (b) solid volume fraction, (c) granular temperature and (d) normalized fluctuation energy flux, Q/p, for e_n = 0.5, p = 2.9 × 10⁻⁵, h = 7.94, and different values of the boundary particle diameter d_w (see the legend in (a)).

The ratio of the shear stress to the pressure, s/p, is the macroscopic friction of the system and measures the capability of the frictionless moving spheres to act as a lubricant fluid with respect to the two rigid, bumpy boundaries. The dependence of the macroscopic friction on the coefficient of restitution, the imposed pressure and the bumpiness is illustrated in Fig. 9(a), where the symbols are the measurements in the DEM simulations and the lines the predictions of the continuum model. The continuum model and the DEM simulations are in remarkable agreement. Similarly, the continuum model correctly predict also the flow height H, that is the gap between the rigid, bumpy planes (Fig. 9(b)).

Fig. 9 Measured (symbols) and predicted (lines) (a) macroscopic friction and (b) flow height as functions of the coefficient of restitution for h = 7.94 and: d_w = 1 and p = 2.9 × 10⁻⁵ (red circles); d_w = 1 and p = 2.9 × 10⁻⁴ (green triangles); d_w = 1 and p = 2.9 × 10⁻³ (orange squares); d_w = 1/2 and p = 2.9 × 10⁻⁵ (purple diamonds); d_w = 1/4 and p = 2.9 × 10⁻⁵ (blue lower triangles).

Counter-intuitively, both the macroscopic friction and the flow height decrease as the particles become more dissipative and the rigid boundaries less bumpy, i.e., e_n and d_w decrease. For instance, we observe a 40% reduction in the macroscopic friction as d_w diminishes from 1 to 1/4. Crucially, low values of s/p and H are associated with a strong spatial heterogeneity of the solid volume fraction, and, in particular, with the tendency of the particles to accumulate in a dense, slow-moving core squeezed in between two regions of high shear and high agitation near the bumpy planes.

The imposed pressure plays a minor role in determining the macroscopic friction (Fig. 9(a)). However, the gap between the bumpy planes decreases as the pressure increases by two orders of magnitude (Fig. 9(b)). In this situations, as previously mentioned, the spatial gradients in both ν and T reduce (Fig. 7), and we observe an increase of the macroscopic friction of about 20% (Fig. 9(a)).

At large pressure (p = 2.9 × 10⁻³) and small restitution coefficient (e_n < 0.5), a persisting steady state cannot be attained even if the total simulation time is increased by one order of magnitude. In such conditions, we observe the formation of regular crystalline structures inside the flowing particles (Fig. 10(a)), that is a dynamic phase-transition from a random to an ordered or partially ordered state. Similar non-steady crystalline structures have been observed in 2D constant-volume shearing simulations, in the absence of physical boundaries, at very large solid volume fractions and high restitution coefficients.³³

Fig. 10 (a) Example of a snapshot taken from the simulation with e_n = 0.2, d_w = 1 and p = 2.9 × 10⁻³ in the steady state, characterized by the presence of crystalline structures inside the flowing particles. (b) Sketch of face centered cubic (FCC) and hexagonal close pack (HCP) structures.

Measurements of the temporal evolution of the fraction of flowing particles that are arranged in either a face centered cubic (FCC) or a hexagonal close pack (HCP) structure (Fig. 10(b)) during shearing permit to quantitatively identify this phase-transition. Indeed, all the simulations with p = 2.9 × 10⁻⁵, p = 2.9 × 10⁻⁴, and those with p = 2.9 × 10⁻³ and e_n ≥ 0.5, reach a steady state with no flowing particles arranged in regular structures, and such stationary configuration persists for more than 20.000 saving time steps. On the other hand, in simulations with p = 2.9 × 10⁻³ and e_n < 0.5, we observe a sudden rise of the fraction of flowing particles involved in regular structures, which rapidly reaches a value of about 30%. To identify particles in FCC and HCP structures, we use the adaptive common neighbor analysis as described in ref. 54, where the standard method^55,56 is extended to multi-phase systems, and implemented in the visualization tool for particle-based systems OVITO.‡⁵⁷

We further investigate the conditions for the phase-transition by performing simulations at imposed pressure larger than 2.9 × 10⁻³, with various coefficients of restitution, while keeping d_w = 1 and h = 7.94, and measuring the fraction of flowing particles arranged in FCC or HCP structures. We point out that the presence of crystals can be either localized in a small sub-region or span the entire domain. Nevertheless, the fraction of moving particles arranged in regular structures is always larger than 10% whenever crystallization occurs, and even reach 80% at large p and small e_n.

We can build, then, the phase-diagram depicted in Fig. 11(a) which indicates that crystallization occurs when the pressure is larger than a critical value, p_c, and that this critical value is an increasing function of the coefficient of restitution. Perhaps intuitively, more dissipative particles crystallize under lower pressures given that the granular temperature is reduced at smaller e_n (Fig. 5(c)), and the randomizing effect of agitation is suppressed.

Fig. 11 (a) Phase-diagram for the onset of crystallization for d_w = 1 and h = 7.94, representing DEM simulations in the presence (filled circles) and in the absence (open circles) of particles arranged in regular structures. The solid line represents the critical pressure predicted from the continuum model. (b) Critical pressure for the onset of crystallization as a function of the coefficient of restitution predicted by the continuum model for h = 7.94 and different values of the diameter of the boundary particles.

Obviously, the continuum model based on the kinetic theory of granular gases, which assumes randomness, cannot apply to crystallized systems. For example, the radial distribution function at contact is singular at the random close packing, while granular solid crystals exist even at larger solid volume fractions. Different expressions of the radial distribution function would permit to overcome this limitation,⁴² but the modelling of shearing crystallized systems is beyond the scope of the present work. Nonetheless, a phase-transition to an ordered, crystalline state is firstly possible when the solid volume fraction exceeds the freezing point, ν = 0.49.⁴¹ If we assume that the crystal nucleation in the flowing particles is induced by the regular geometry of the bumpy boundaries, we can then estimate that the critical pressure at which crystallization occurs corresponds to the pressure at which the solid volume fraction at the rigid boundaries, ν_w, obtained by solving numerically the differential equations of the continuum model, is equal to 0.49. Operatively, we add the boundary condition ν (z = 0) = ν_w = 0.49 to those of (12)–(17), and let p = p_c to be an unknown determined as part of the solution of the system of differential eqn (8)–(11). As shown in Fig. 11(a), the critical pressure evaluated from the continuum model on the basis of the condition ν_w = 0.49 permits to satisfactorily predict the onset of crystallization in the DEM simulations.

Interestingly, the critical pressure predicted from the continuum model shows a non-monotonic dependence on the coefficient of restitution for d_w larger than 1, but basically decreases as d_w decreases (Fig. 11(b)).

5 Conclusions

We have performed discrete numerical simulations of steady flows of identical, frictionless, inelastic spheres sheared between parallel, bumpy planes, under pressure-imposed conditions, by systematically varying the coefficient of restitution, the imposed pressure and the bumpiness of the boundaries. We have used the results of our simulations to test the predictions of the extended kinetic theory of granular gases.

In particular, we have first validated the constitutive relations of the pressure, the shear viscosity, the rate of collisional dissipation of the fluctuation kinetic energy and the pseudo-thermal diffusivity, for solid volume fractions smaller than the random close packing. These relations do not involve fitting parameters, but only microscopic contact properties of the particles and material parameters that have a clear physical meaning: the random close packing and the solid volume fraction at the freezing point (0.64 and 0.49 for frictionless spheres, respectively).

Then, we have phrased and numerically solved a two-point boundary-valued problem based on the constitutive relations of kinetic theory and the balance of fluctuation kinetic energy, by employing appropriate boundary conditions for the rates of supply of momentum and fluctuation energy into the flow by the bumpy planes. Such boundary conditions have been originally derived by Richman²⁴ for nearly elastic spheres, by means of statistical mechanics arguments, and later empirically modified by Berzi and Vescovi²⁵ on the base of numerical simulations. Here we have slightly revised the rate of supply of fluctuation energy to account for the more dissipative particles employed in the present DEM simulations.

We have compared the predicted spatial distributions of mean velocity, solid volume fraction, granular temperature and fluctuation energy flux against the results of the DEM simulations, and have shown remarkable qualitative and quantitative agreement, for all the combinations of input parameters investigated. The dependence of the macroscopic friction, i.e., the ratio of the shear stress to the pressure, on the restitution coefficient, the imposed pressure and the bumpiness of the planes is also very well captured by the extended kinetic theory. The small discrepancies between the theoretical predictions and the data are likely due to the empirical nature of the two key physical quantities which affect the theory: the radial distribution function at contact and the correlation length. Direct measurements of such quantities by means of numerical simulations would allow to improve these expressions and, as a consequence, to provide more accurate predictions.

We have observed that the shear resistance of the flow decreases as the particles become more dissipative and the boundaries less bumpy. As the coefficient of restitution decreases, spatial heterogeneities in the solid volume fraction and granular temperature intensify: the particles tend to accumulate in a dense, slow-moving central core sandwiched in between two more dilute regions of high shear and high agitation near the bumpy planes.

Finally, our simulations have shown that at large imposed pressures, the granular system becomes very dense and crystallizes. The phase-transition from a disordered/fluid to an ordered/crystalline state takes place at a critical value of the imposed pressure, which depends on the restitution coefficient. The extended kinetic theory can predict the critical pressure at which crystallization occurs by assuming that it corresponds to the pressure at which the solid volume fraction at the bumpy boundaries reaches the freezing point. The dependence of the critical pressure on the friction coefficient and a detailed characterization of the behaviour of the granular materials in the sheared crystalline state will be the subject of future works.

Author contributions

All authors designed research. DV performed simulations and evaluated the data. DV and DB performed the analytical calculations and wrote the manuscript. All authors discussed the results, read, revised and approved the final manuscript.

Data availability

Data for this article are available at Zenodo at: https://doi.org/10.5281/zenodo.13785775.

Conflicts of interest

There are no conflicts to declare.

Appendix A: coarse-graining procedure

Considering each particle i to be a sphere of volume W_i, position rⁱ = (rⁱ_x,rⁱ_y,rⁱ_z) and velocity Vⁱ = (Vⁱ_x,Vⁱ_y,Vⁱ_z), then volume fraction, ν, flow velocity, v = (u,v,w), and granular temperature, T, at the center of the local, averaging volume W, are computed aswhere the summation over i extends over all particles which are wholly contained in W. We have introduced the fluctuation velocity of particle i with respect to the (local) mean velocity: Ṽⁱ = (Ṽⁱ_x,Ṽⁱ_y,Ṽⁱ_z) = Vⁱ − v.

Pressure, shear stress and fluctuating energy flux are calculated as the sum of a “kinetic” and a “contact” contribution:

p = p^K + p^C, (24)

s = s^K + s^C, (25)

Q = Q^K + Q^C. (26)

The kinetic contributions to p, s and Q are given aswhereas contact contributions areIn eqn (30)–(32), the summation over i and then over j > i extends over all pairs of particles which are in contact and such that the centroids of the contact are located within W, counting each pair of particles only once, whereas b^ij = (bⁱ_x,bⁱ_y,bⁱ_z) = rⁱ − r^j and f^ij = (f^ij_x,f^ij_y,f^ij_z) are the branch vector and the interaction force between particles i and j at contact, respectively.

Finally, the rate of collisional dissipation of the fluctuation kinetic energy is computed as

We emphasize that the above described coarse graining is performed inside the domain at a distance of at least one particle diameter from the edges of the boundary particles. This permits to avoid complications arising from considering different particle sizes, and, as a consequence, W_i = π/6 for all particles i.

Acknowledgements

The authors would like to thank Dr Eivind Bering for writing the first version of the LAMMPS script. This project has received funding from the HORIZON-EIC-2021-PATHFINDEROPEN-01 N. 101046693, SSLiP project. Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or European Commission. Neither the European Union nor the granting authority can be held responsible for them.

Notes and references

1. P. Cundall and O. Strack, Geotechnique, 1979, 29, 47–65.
2. G. D. R. MiDi, Eur. Phys. J. E: Soft Matter Biol. Phys., 2004, 14, 341–365.
3. F. da Cruz, S. Emam, M. Prochnow, J. Roux and F. Chevoir, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 021309.
4. K. Kamrin and G. Koval, Phys. Rev. Lett., 2012, 108, 1–5.
5. S. Chialvo and S. Sundaresan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021305.
6. A. Singh, V. Magnanimo, K. Saitoh and S. Luding, New J. Phys., 2015, 17, 043028.
7. D. Berzi and J. Jenkins, Soft Matter, 2015, 11, 4799–4808.
8. D. Vescovi and S. Luding, Soft Matter, 2016, 12, 8616–8628.
9. H. Shi, S. Roy, T. Weinhart, V. Magnanimo and S. Luding, Granular Matter, 2020, 22, 1–20.
10. D. Berzi, Phys. Rev. Fluids, 2024, 9, 034304.
11. J. Jenkins and S. Savage, J. Fluid Mech., 1983, 130, 187–202.
12. C. Lun, J. Fluid Mech., 1991, 223, 539–559.
13. V. Garzó and J. Dufty, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 59, 5895–5911.
14. I. Goldhirsch, Annu. Rev. Fluid Mech., 2003, 35, 267–293.
15. N. Mitarai and H. Nakanishi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 031305.
16. J. Jenkins, Phys. Fluids, 2006, 18, 103307.
17. J. Jenkins, Granular Matter, 2007, 10, 47–52.
18. J. Jenkins and D. Berzi, Granular Matter, 2010, 12, 151–158.
19. D. Berzi, Acta Mech., 2014, 225, 2191–2198.
20. D. Gollin, D. Berzi and E. Bowman, Granular Matter, 2017, 19, 56.
21. J. Jenkins and M. Larcher, Phys. Rev. Fluids, 2023, 8, 024303.
22. M. Islam, J. Jenkins and S. Das, J. Fluid Mech., 2022, 950, A13.
23. D. Vescovi, D. Berzi, P. Richard and N. Brodu, Phys. Fluids, 2014, 26, 053305.
24. M. Richman, Acta Mech., 1988, 75, 227–240.
25. D. Berzi and D. Vescovi, Comput. Part. Mech., 2017, 4, 373–377.
26. D. Risso and P. Cordero, J. Stat. Phys., 1996, 82, 1453–1466.
27. R. Garcia-Rojo, S. Luding and J. Brey, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 061305.
28. E. Khain, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 051310.
29. E. Khain, Phys. Rev. E, 2018, 98, 012903.
30. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in‘t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
31. C. Campbell and C. Brennen, J. Fluid Mech., 1985, 151, 167–188.
32. M. Louge, Phys. Fluids, 1994, 6, 2253–2269.
33. M. Alam and S. Luding, Phys. Fluids, 2003, 15, 2298–2312.
34. M. Babic, Int. J. Eng. Sci., 1997, 35, 523–548.
35. M. Babic, Phys. Fluids, 1997, 9, 2486–2505.
36. T. Weinhart, A. Thornton, S. Luding and O. Bokhove, Granular Matter, 2012, 14, 289–294.
37. T. Weinhart, R. Hartkamp, A. Thornton and S. Luding, Phys. Fluids, 2013, 25, 070605.
38. B. Miller, C. O'Hern and R. Behringer, Phys. Rev. Lett., 1996, 77, 3110–3113.
39. P. Zamankhan, T. Tynjala, W. Polashenski, P. Zamankhan and P. Sarkomaa, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 7149–7156.
40. N. Carnahan and K. Starling, J. Chem. Phys., 1969, 51, 635–636.
41. S. Torquato, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, 3170–3182.
42. S. Luding, Nonlinearity, 2009, 22, R101–R146.
43. S. Chialvo, J. Sun and S. Sundaresan, Phys. Fluids, 2013, 25, 070603.
44. D. Berzi, J. Jenkins and P. Richard, J. Fluid Mech., 2020, 885, A27.
45. D. Berzi and D. Vescovi, Phys. Fluids, 2015, 27, 013302.
46. J. Jenkins, S. Myagchilov and H. Xu, Unpublished.
47. E. Lim, Eur. Phys. J. E: Soft Matter Biol. Phys., 2010, 32, 365–375.
48. M. Alam and P. Nott, J. Fluid Mech., 1998, 377, 99–136.
49. P. Nott, M. Alam, K. Agrawal, R. Jackson and S. Sundaresan, J. Fluid Mech., 1999, 397, 203–229.
50. M. Alam, V. Arakeri, P. Nott, J. Goddard and H. Herrmann, J. Fluid Mech., 2005, 523, 277–306.
51. M. Alam, P. Shukla and S. Luding, J. Fluid Mech., 2008, 615, 293–321.
52. P. Shukla and M. Alam, Phys. Rev. Lett., 2009, 103, 068001.
53. P. Shukla and M. Alam, J. Fluid Mech., 2011, 666, 204–253.
54. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2012, 20, 045021.
55. J. Honeycutt and H. Andersen, J. Phys. Chem., 1987, 91, 4950–4963.
56. D. Faken and H. Jonsson, Comput. Mater. Sci., 1994, 2, 279–286.
57. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.

---

Footnotes

† https://www.lammps.org.

‡ https://www.ovito.org.

---