Dense , inhomogeneous shearing flows of spheres

We make use of recent extensions of kinetic theory of granular gases to include the role of particle stiffness in collisions to deal with pressure-imposed shearing flows between bumpy planes in relative motion, in which the solid volume fraction and the intensity of the velocity fluctuations are not uniformly distributed in the domain. As in previous numerical simulations on the flow of disks in an annular shear cell, we obtain an exponential velocity profile in the region where the volume fraction exceeds the critical value at which a rate-independent contribution to the stresses arises. We also show that the thickness of the inertial region, where the solid volume fraction is less than the critical value, and the shear stress at the moving boundary are determined functions of the relative velocity of the boundaries.


Introduction
Recently [1], it has been quantified how the particle stiffness influences the frequency of collisions responsible for the rate-dependent components of particle stresses, and induces the development of rateindependent components to the stresses associated with permanent deformations.That theory has been successfully tested against discrete numerical simulations of simple shearing, in which the velocity fluctuations and the solid volume fraction were uniformly distributed.It has been also briefly shown [1] how that theory can be applied to nonhomogeneous shearing flows, to explain the presence of an exponentially decaying velocity profile in very dense regions (creeping flow [2,3]).Here, we use a simplified version of the theory of ref. [1] to solve for the flow of identical spheres sheared between two parallel, bumpy planes, one moving at constant velocity and one at rest.To mimic the results of previous discrete numerical simulations of disks in an annular shear cell [4], we take the pressure constant in the flow and the shear stress inversely proportional to the square of the distance from the moving boundary.Suitable boundary conditions allow us to obtain solutions to the flow that are in qualitative agreement with the numerical simulations.

Theory
An assembly of identical spheres of mass density ρp and diameter d is sheared between two parallel planes, with the upper plane moving at constant velocity, in the absence of gravity.The two planes are made bumpy by gluing a layer of spheres (identical to those of the flow) in a regular, hexagonal array.The bumpiness ψ of the planes is, then, determined by the mean interparticle distance between the glued grains [5].We take x and y to be the directions parallel and perpendicular to the planes, respectively.The x-component of the particle velocity, the only one present, is u, while ν is the solid volume fraction.The distance, R, between the planes and the velocity, V, of the moving plane are the control variables of the system.The flow configuration is depicted in Fig. 1.To mimic the annular shear cell of Koval et al [4] with the above described configuration, we take the distributions of particle stresses to be those determined from the numerical simulations: a constant particle pressure p; and for the particle shear stress, ( ) ¼ , where S is the particle shear stress at y = R.As observed in the simulations, the flow domain can be in general divided into two sub-regions: an inertial region adjacent to the moving boundary of thickness δ and a quasi-static region of thickness R-δ (Fig. 1).Koval et al. [4] distinguished the two regions on the basis of the local value of the stress ratio s/p being larger or lower than the minimum (yield) value in simple shearing.Here, based on the analysis reported in ref. [1], we identify the inertial (quasi-static) region with the region where the solid volume fraction is less (larger) than the critical value, νc, necessary to develop elastic stresses associated with permanent deformations of the aggregate [6].

Inertial region
In the inertial region, we use the constitutive relations for the pressure and the shear stress of kinetic theory [7] in the dense limit [8], and where: e is the coefficient of normal restitution (the negative of the ratio of pre-to post-collisional relative velocity along the line of centres of two colliding grains); g0 is the radial distribution function at contact, here taken to be the dense limit of that suggested in ref. [9], ( ) where Q is the flux of fluctuation energy: with M = (1+e)/2+9π(1+e) 2 (2e-1)/[128-56(1-e)]; and Γ is the collisional rate of dissipation: ( ) where ε = e-3/2μ exp(-3μ) is an effective coefficient of restitution that takes into account the role of particle surface friction, μ [10].The quantity L in eq. ( 5) is the correlation length of extended kinetic theory [11,12] to take into account the breaking of molecular chaos at volume fractions larger than 0.49.Here, as in ref. [13], we use: Without loss of generality, we now take the particle density, diameter and pressure to be unity.As in ref. [8], we derive eq. ( 1) with respect to y, and use the expression of the radial distribution function at contact and eq.( 4) to obtain a differential equation for the solid volume fraction, where the subscript indicates the derivative with respect to ν and T is, from eq. ( 1), ( ) A differential equation for the particle velocity is obtained by inverting eq. ( 2), Finally, the differential equation for the energy flux is, from eq. ( 3) with eqs.( 5) and ( 8), ( ) where L is calculated from eqs. ( 6) and (8).The system of the 3 differential equations ( 7)-( 9) allows the determination of the distribution of ν, u and Q in the inertial region once appropriate boundary conditions are provided.At y = R, we use the boundary conditions derived by Richman [5] for the flows of spheres over bumpy planes: a slip velocity given by where B = 12π/2 1/2 [1+5/(8νg0)]; and an energy flux given as In eqs.( 10) and (11), TR and LR are the granular temperature and correlation length evaluated at y = R.At y = δ (the interface between the inertial and the quasi-static region), we use the boundary condition for the energy flux at an erodible boundary [14], modified with the introduction of the correlation length, ( ) where TR-δ and LR-δ are the granular temperature and correlation length evaluated at y = δ.We also fix the solid volume fraction there to be equal to a given value νE, less than νc.We cannot use νc as a boundary condition for the solid volume fraction, because the radial distribution function at contact is singular there, and so would be the stresses [15].Including the role of the particle stiffness in collisions would allow for the stresses to be non-singular at ν = νc [1].Taking νE less than νc is a simpler, but equivalent, alternative to that approach.Finally, we set the value uE of the particle velocity at y = δ, based on the analytical solution to the flow in the quasi-static region (see next Section).The five boundary conditions the solution (using the Matlab ® routine 'bvp4c') of the system of differential equations ( 7)-( 9) and the determination of the thickness δ and the value of the shear stress S at y = R as part of the solution.

Quasi-static region
In the quasi-static region, i.e., y δ, we assume, as in ref. [1], that the stresses have an elastic component, associated with permanent deformations of the particles, and an inertial component, associated with the momentum exchange due to the fluctuating motion of the particle centres of mass.We assume that the ratio of the shear stress to the pressure is equal to the ratio of eq. ( 2) to eq. ( 1) even in the quasi-static region (this allows the reproduction of the results of numerical simulations of simple shearing [1]).With p = d =1, this gives ( ) We also assume that the fluctuation energy balance in this quasi-static region reduces to a balance between the diffusion and the collisional dissipation of energy, i.e., that the energy production due to the work of the inertial component of the shear stress is negligible.The constitutive relations for the energy flux and the dissipation rate are ( ) respectively, where E is the Young's modulus of the particles.For simplicity, we take L in eq. ( 15) to be equal to the its value in simple shearing [1], ( ) Using eqs.( 14)-( 16) in the energy balance, and neglecting the work of the shear stress, gives where ( ) where, as already mentioned, TR-δ is the granular temperature at y = R-δ.Using eq. ( 18) in eq. ( 13), with the distribution of the shear stress with the distance from the planes measured in the numerical simulations, and integrating, with the boundary condition u = 0 when y = 0, gives When y = R-δ, eq. ( 19) provides the value of the velocity uE to be used as boundary condition in the numerical solution of the inertial region (see previous Section): ( ) ( )

Results and conclusions
We now obtain the results of the present theory for a representative set of parameters.We take: e = 0.7 and μ = 0.5, so that νc = 0.587 [6]; νE = 0.586 and ψ = π/6.We use different values of total cell thickness R and plane velocity V and make qualitative comparisons with the numerical simulations on disks in an annular shear cell [4].
Fig. 2 shows the particle velocity, scaled with the velocity of the moving plane, as a function of the distance from the moving plane.As in ref. [4], it is possible to appreciate the slip velocity at y = R and the exponential behaviour (eq.( 19)) in the quasi-static region; when R = 50 and V = 2.5, it is also possible to notice the change in the concavity of the velocity profile in the inertial region.The thickness of the inertial region, marked by the knee in the velocity profiles of Fig. 2, apparently increases with increasing R and V.
The behaviour of δ with V is actually non-monotonic, as shown in Fig. 3 for the case R = 50.The thickness δ is zero when V is less than 0.1, and has a maximum when V is near 2.5; for large values of V, it decreases slightly and seems to saturate for even larger values.This is qualitatively similar to the relation between the shear stress at y = R and the velocity of the moving plane depicted in Fig. 4, although in that case the decrease after the peak is much more pronounced.Koval et al. [4] obtained monotone increase of δ and S with V.However, they did not perform simulations for V larger than 2.5, and their figures 9a and 3b seem to confirm that the measurements are approaching a maximum at that value.Although the results of the present theory are in good qualitative agreement with 2D discrete numerical simulations, 3D discrete numerical simulations are required to fully test the capability of the theory developed in ref. [1] to reproduce dense granular flows in nonhomogeneous configurations.We postpone these comparisons to future work.

Fig. 1 .
Fig. 1.Sketch of the flow configuration with the frame of reference.
to fit the results of numerical simulations of simple shearing when e is less than 0.95; T is the granular temperature (one third of the mean square of the velocity fluctuations); and J = (1+e)/2+π(1+e) 2 (3e-1)/[96-24(1-e) 2 -20(1-e 2 )].Here, and in what follows, a prime indicates a derivative with respect to y.The fluctuation energy balance is

Fig. 3 .
Fig. 3. Thickness of the inertial region as a function of the velocity of the moving plane for R = 50.

Fig. 4 .
Fig. 4. Shear stress at y = R as a function of the velocity of the moving plane for R = 50.