Constitutive equations for granular flow with uniform mean shear and spin fields

Numerical simulations of two-dimensional granular flows under uniform shear and external body torque were performed in order to extract the constitutive equations for the system. The outcome of the numerical simulations is analyzed on the basis of the micropolar fluid model. Uniform mean shear field and mean spin field, which is not subordinate to the vorticity field, are realized in the simulations. The estimates of stresses based on kinetic theory by Lun [Lun, J. Fluid Mech., 1991, 233, 539] are in good agreement with the simulation results for a low area fraction $\nu=0.1$ but the agreement becomes weaker as the area fraction gets higher. However, the estimates in the kinetic theory can be fitted to the simulation results up to $\nu=0.7$ by renormalizing the coefficient of roughness. For a relatively dense granular flow ($\nu=0.8$), the simulation results are also compared with Kanatani's theory [Kanatani, Int. J. Eng. Sci., 1979, 17, 419]. It is found that the dissipation function and its decomposition into the constitutive equations in Kanatani's theory are not consistent with the simulation results.


Introduction
Collective motions of granular materials behave like fluid motions under appropriate conditions. However, unlike the Newtonian fluids, the basic equations for the collective motions of granular materials have not been well established yet. One of the difficulties of the problem lies in the fact that the scale of macroscopic collective motions is not well separated from the microscopic scale of the system such as the radius of the granular particles. Thus, applicability of arguments based on the scale separation would be limited. Many detailed properties of the individual particles would directly affect the behavior of the macroscopic flow.
One possible way to treat such granular flows is to model them as flows of a micropolar fluid, a fluid with polar micro-structures such as spin [1,2]. By applying micropolar fluid mechanics to granular flows, the spin of the granular particles can be coupled to the dynamics of the macroscopic collective motions of the granular particles. The microscopic properties of the granular particles are reflected in the equations of motion for the macroscopic fields through the constitutive equations, i.e., the relations between strains and stresses (see section 2).
For sparse and rapid granular flows, the equations of motion as a micropolar fluid can be derived within the framework of a kinetic theory. Firstly, the kinetic theory was developed without introducing the frictional interactions between particles and the degrees of freedom for spin by Savage and Jeffrey [3], and Jenkins and Savage [4]. Although Jenkins and Richman [5], Jenkins and Zhang [6] and Yoon and Jenkins [7] introduced frictional interaction between particles to the kinetic theory, they eliminated the macroscopic degrees of freedom of the spin field by assuming that the macroscopic spin field is subordinate to the vorticity field. Such an assumption may be justified, for example, for steady flows far from the boundary. In the kinetic theories, the effect of frictional interactions can be absorbed into the renormalized restitution coefficient. Saitoh and Hayakawa [8] performed numerical simulations of two-dimensional granular flow under a plane shear and confirmed that the hydrodynamic equations derived from the kinetic theories agree with the simulation results. In some cases, such as flows near boundaries, discrepancy between the spin field and the vorticity field is not negligible. The kinetic theories retaining the spin field as independent macroscopic degrees of freedom were developed by Lun and Savage [9], Lun [10] and Goldshtein and Shapiro [11]. Mitarai et al. [12] performed numerical simulation of a collisional granular flow on a slope and showed that the velocity and spin field profiles are in agreement with the micropolar fluid equations based on constitutive equations which are consistent with that in [10].
When the granular particles become dense enough and the volume fraction exceeds the critical value, the collective motions of particles stop to behave like a fluid in a sense that a finite shear stress is required to create an infinitesimal strain. Such a phase is called the jammed phase. The phase transition between unjammed phase and jammed phase is called the jamming transition. Scaling laws near the critical point of the jamming transition have been suggested and verified in the numerical simulations by Hatano [13], and Otsuki and Hayakawa [14,15]. The frictional interactions among particles were not considered in their studies. Recently, a number of results on the jamming transition based on numerical simulations including the frictional interactions have been reported (e.g., Silbert et al. [16], Zhang and Makse [17] and Shundyak, Hecke and Saarloos [18].) There can be a substantial intermediate regime of the volume fraction between the kinetic region with low volume fraction and the critical region near the jamming transition. In this regime, interaction of n-particles with n > 2 would become important. Kanatani [19] developed a micropolar fluid theory for relatively dense granular flows in which particles are almost regularly in contact with the other particles. The regime where Kanatani's theory is applicable is possibly located in this intermediate regime. Kano et al. [20] showed that numerical simulation of a granular flow on an inclined trough is in qualitative agreement with the micropolar fluid equation based on Kanatani's theory regarding the velocity profile.
In this paper, we focus on the constitutive equations for granular flows. As an intrinsic nature of the granular flows, the spin field associated with granular particles is not subordinate to the velocity field of their mean flow. This situation is analogous to the case of micropolar fluids in which the spin field is regarded as an independent degree of freedom. As we will see in section 2, both the difference between the vorticity and spin fields, which will be denoted by R ji , and the spatial derivative of the spin field, which will be denoted by Ω kji , contribute to the constitutive equations. From a theoretical point of view, it is desirable to analyze them separately. Therefore, let us consider the case of R ji = 0 and Ω kji = 0. Note that R ji = 0 near the boundary. When sufficient numbers of particles are contained in a region near the boundary with the length scale smaller than the typical length scale in which the shear and spin fields changes, we may consider that uniform shear and spin fields (i.e. Ω kji = 0) with R ji = 0 are approximately realized in the region. The situation R ji = 0 and Ω kji = 0 would be also obtained by applying external torque to each particle. Note that, provided that the micropolar fluid picture is appropriate for the granular flows, the constitutive equations depend solely on velocity, spin fields and their spatial derivatives at the local point under consideration and independent of driving forces that generate the fields. A possible way to apply the external torque to each particle in experiments is to embed the source of angular momentum inside each particle. That is, the particle is supposed to be a kind of micro-machine composed of an outer shell and an inner sphere with the friction between them being small. Initially, the inner sphere is made to rotate with a high angular velocity by some means while the outer shell is not rotating. Then, the angular momentum of the inner sphere is continuously supplied to the outer shell through the friction until the inner sphere loses its substantial angular momentum. By virtue of small friction, one can realize a longer period for the experiment. By considering the inner sphere as an exterior system, the situation implies that the external torque is continuously applied to the particle (the outer shell). The inner sphere can be replaced by a liquid with low viscosity such as a super fluid. The actual setting of the above system for the experiment may be quite difficult. However, in numerical experiments, it is quite easy to apply the external torque to each particle.
Taking the above into consideration, we performed numerical simulations of two-dimensional granular flows under uniform shear and uniform external torque field. By virtue of the external torque field and the applied boundary conditions, macroscopically uniform vorticity and spin fields are realized and their magnitudes are controlled independently, which means that Ω kji = 0 and the magnitude of R ji can be controlled (see section 6). Thus, we concentrate on the R ji dependence of the constitutive equations with Ω kji fixed to 0. The study of the Ω kji dependence of the constitutive equations will be the next step and will not be referred to in this paper. Unlike the preceding numerical studies such as [20] and [12], we are able to obtain not only the velocity and spin field profiles, which are the results of the constitutive equations, but also the constitutive equations directly. Since the subject of this paper is the micropolar fluid aspect of the granular flows, the value of area fraction is varied within the unjammed region. We compare the results from the numerical simulations with those from the kinetic theory by Lun [10], which is capable of treating cases that the spin field is not subordinate to the vorticity field. For the intermediate regime noted above, we also compared the simulation results with Kanatani's theory [19].
This paper is organized as follows. In section 2, a brief review of the micropolar fluid theory is given. In sections 3 and 4, the kinetic theory by Lun and Kanatani's theory are reviewed, respectively. In section 5, comments on the two theories are given. In section 6, the results of the numerical simulations are shown and they are compared with the theories. In section 7, discussion is presented.

Equations for micropolar fluid
In this paper, we consider collective motions of particles. For simplicity, we assume that the particles are of the same mass m and the same moment of inertia I. Let c i (t), w ji (t) and r i (t) be, respectively, the velocity, the spin and the position of the particle at time t where i and j are coordinate indices of d-dimensional space. Here, d can formally be an arbitrary positive integer larger than 1. In this paper, we use the convention that the spin is expressed by a skew-symmetric tensor w ji whose (j, i)-th component gives the angular velocity in the (j, i) coordinate plane. Let F (N ) (c (1) , w (1) , r (1) ; · · · ; c (N ) , w (N ) , r (N ) ; t) be the probability density function in the phase space of N -particles system satisfying Liouville equation. Here, the bold letters c, w and r denote vector or tensor, the superscript (α) on c (α) , w (α) and r (α) is the index of the particle and F (N ) is symmetrized with respect to interchanges of the particles. The s-particles set distribution function f (s) (s N ) is given by (1) , w (1) , r (1) ; · · · ; c (N ) , w (N ) , r (N ) ; t), (2.1) and the number density of the s-particles sets n (s) is given by Macroscopic fields such as the mass density field ρ(r, t), the moment of inertia density field ρ I (r, t), the velocity field v(r, t) and the spin field ω(r, t) are introduced as ρ(r, t) := mn (1) (r, t), ρ I (r, t) := In (1) (r, t), for an arbitrary function ψ of c and w. Hereafter, indices or subscripts of spatial or time coordinates will be suppressed unless we need to emphasize them.
These macroscopic fields satisfy the following equations, where D/Dt = ∂/∂t+v i ∂ i , σ ji is the stress tensor, λ kji the couple-stress tensor, b i the external body force and τ ji the external body torque. Here and hereafter, summations are taken over repeated subscripts and we employ the notations, for an arbitrary tensor T ji . Equations (2.6)-(2.9) correspond, respectively, to the conservation laws of mass, moment of inertia, momentum and angular momentum. Note that v i and ω ji are mutually independent degrees of freedom. When there is no spin field ω ji , (2.6)-(2.8) reduce to the equations of ordinary fluids. A fluid with the degree of freedom of the spin field is called micropolar fluid. The equations of motions will be closed when the stresses σ ji and λ kji are known in terms of ∂ j v i and Ω kji := ∂ k ω ji . The equations which yield such relations are called constitutive equations. Let C and W be the fluctuations of velocity c and spin w of individual particles around the macroscopic fields, i.e., The "internal energies" per unit mass ǫ t and ǫ r associated with the fluctuations C and W , respectively, are given by Here, the subscripts t and r indicate the translational and rotational motions, respectively. Two kinds of "temperature", T t and T r , are introduced by the relations ǫ t = (d/2)T t and ǫ r = (d(d − 1)/4)T r . The total "internal energy" is given by ǫ U = ǫ t + ǫ r . Here, we have assumed that the particles are rigid bodies and that there is no potential force acting among particles, that is, there is neither contribution of elastic energy nor of potential energy to the "internal energy". One can show from (2.6)-(2.9) that the kinetic energy of macroscopic fields per unit mass, 14) with where we have introduced notations, and p = −σ ii /d is the pressure. Note that Ψ is the work done on a fluid element per unit volume per unit time by the stress σ ji , the couple stress λ kji , the external body force b i and the external body torque τ ji . From (2.15) and the energy budget equation, where q is input of the "internal energy" per unit volume per unit time, h j is the "internal energy" flux. The quantity Φ − p∂ i v i gives the energy transferred from the kinetic energy to the "internal energy" per unit volume per unit time. When macroscopic fields ρ, ∂ j v i , ǫ U and h j are constant in time and space, we have ∂ i v i = 0 from (2.6), and (2.20) reduces to The equation (2.21) implies that the energy transfer rate Φ from the kinetic energy to the "internal energy" balances with the dissipation rate of "internal energy" −q for macroscopically steady states. Φ is called the dissipation function since it gives the energy going out from the kinetic energy ǫ K per unit time per unit volume in the macroscopically steady state.

Kinetic theory for collisional granular flow
A Kinetic theory for collisional granular flow is developed by Savage and Jeffrey [3], and Jenkins and Savage [4]. It is extended to include the effects of surface friction and inertial moment of particles by Lun and Savage [9], and Lun [10]. From assumptions on the collisional process of two particles, the constitutive equations for the granular material as a micropolar fluid are derived. We briefly review the theory by Lun [10] in the following.
The dimension of the configuration space is d and the particles are assumed to be d-dimensional spheres with the same radius a. The mass and the moment of inertia are, respectively denoted by m and I as in section 2. Let the colliding two particles be labeled as 1 and 2, and a unit vector k be defined by figure 1). In this section, we employ the notation that the quantities without (with) checkď enote those just before (after) the contact. Let J be the impulse of the force exerted on the particle 2 by the particle 1 through the contact point. We have

Iw
(1) Let us assume that the velocity difference ξ i between the two surfaces at the contact point, changes as, after the collision. Here, the coefficient of restitution e and the coefficient of roughness β are assumed to be constants satisfying 0 e 1 and −1 β 1.
Assuming the binary collisions, the equation of motion for ψ , where ψ is an arbitrary function of c and w, can be written as By substituting ψ = m and ψ = mc i in (3.10), one finds that the stress tensor σ ji can be written as ji denotes the kinetic contribution to the stress σ ji due to the particles that cross the plane perpendicular to j-axis, and σ (c) ji denotes the collisional contribution due to the collisions of the two particles in different sides of the plane.
The particle-pair distribution function f (2) is assumed to be approximated by the form where g 0 (r ′ ; r, t) is a radial distribution given by with an assumption that it is insensitive to the direction of a unit vector e. Assuming that a is sufficiently smaller than the typical spatial scale in which the amplitude of f (1) varies, we have Let f (1) be written as where f (1) 0 is the distribution function at a local equilibrium state and φ represents the perturbation.
is given by Note that different temperatures T t and T r are assigned to transitional and rotational degrees of freedom, respectively. Assuming that the mean is much smaller than the typical magnitude of the fluctuation for the distribution of w ji , i.e., |ω ji | ≪ (m/I)T r , we have It is assumed that the function φ is approximated by a linear function of degrees of nonequilibrium such as the symmetric and traceless kinetic stress tensor σ In what follows, we restrict ourselves to the case of macroscopic fields ∂ j v i , ω ji , T t and T r being constant in space and time. In such a case, there are no energy fluxes h where the factor −1/2ρT t is determined from the consistency condition that (3.14) is satisfied. Note that the total σ (k) is given by where there is no anti-symmetric part σ (k) [ji] in the definition (3.14). By substituting ψ = mC i C j /2 into (3.10), we obtain where σ (k) and σ (c) are given by (3.14) and (3.15), respectively, and χ(·) by (3.12). σ (c) and χ(mC i C j /2) reduce to functions of ∂ j v i , ω ji and σ (k) , after performing the integrations with respect to k, c (i) and w (i) (i = 1, 2) in (3.11) and (3.12). Substitute the functions into (3.25) and assume that the magnitudes of ∂ j v i and ω ji are so small that the second or higher order terms in ∂ j v i and ω ji can be neglected. Then we arrive at the equations, and ν is the volume fraction which may be related to the number density n of the particle as (3.30) Now, the constitutive equation for the total stress tensor σ is given by with (3.26)-(3.28). From symmetry, we have µ kji (E ji , R ji , Ω kji = 0) = 0 1 . The obtained constitutive equations are a restricted version of those in [10] in the sense that constants ∂ j v i , ω ji , T t and T r are assumed, and they are a generalized version in the sense that d is arbitrary whereas d = 3 in [10].

Kanatani's theory
A micropolar fluid theory for relatively dense granular flows was developed by Kanatani [19]. The theory can be outlined as follows. As in section 3, a system of the particles with the same radius a, mass m and moment of inertia I is considered. Let the two contacting spherical particles be labeled as 1 and 2. The unit vector k is the same as (3.1) (see figure 1). From (2.12), the quantities associated with particles, c (α) i and w (α) ji (α = 1, 2), can be related to the macroscopic fields as c ji denote the fluctuations around the macroscopic fields. Since |r (1) −r (2) | = 2a is small compared to the typical scale in which the amplitudes of fields vary, v i (r (2) ) and ω ji (r (2) ) can be approximated by its Taylor series about x (1) up to the first-order terms, i.e., v i (r (2) where we have omitted writing the argument and Ω kji (x (1) ). From (3.4), (4.1) and (4.2), the velocity difference ξ i between the two surfaces at the contacting ji can be neglected and that k i is isotropically distributed. Then, we have and · · · denote the average over k i . Here, we have used which can be derived using the general expressions of isotropic tensors and k i k i = 1. Suppose that the average energy dissipation per unit time at a contacting point is estimated by µf ξ, where µ is the kinetic friction coefficient and f is the average amplitude of the force applied in the direction k at the contacting point. Let N c be a number of contacting points per a particle. Then, the energy dissipation per unit volume per unit time −q is given by where the factor 1/2 is introduced to cancel the double counting of the contacting points. Note that the pressure p c due to the contacting force can be estimated by where ν is the volume fraction and S is the surface area of the particle. From The estimate of p c is given as a homogeneous function of E ji , R ji and Ω kji , i.e., p c (aE ji , aR ji , aΩ kji ) = a ζ ′ p c (E ji , R ji , Ω kji ) with a constant ζ ′ . The form of the function depends on whether the flow is slow or fast, and here we omit the review. See the original paper [19] for the details. Let us restrict ourselves to macroscopically uniform and steady states. From (2.21) and (4.9), one obtains which is a homogeneous function of E ji , R ji and Ω kji , i.e., Φ(aE ji , aR ji , aΩ kji ) = a ζ Φ(E ji , R ji , Ω kji ) with ζ = ζ ′ + 1. With the help of Euler's homogeneous function theorem, the choice is made for constitutive equations which are consistent with (2.17).

Comments on the kinetic theory and Kanatani's theory
In the kinetic theory by Lun, the binary collision is assumed, i.e., n-particle collisions for n > 2 are neglected, and the particle-pair distribution f (2) is assumed to be approximated by a product of g 0 and f (1) as (3.16). These assumptions seem to be valid for small volume fraction, ν ≪ 1. This is because, in such a case, n-particle collisions for n > 2 would be rare and the velocities of the two particles before a binary collision would be statistically almost independent. Furthermore, macroscopic fields ∂ j v i and ω ji are assumed to be small in comparison to their microscopic fluctuations. All of these assumptions would become inappropriate with the increase of ν.
Even when ν is small, there are a few points one should consider carefully in the kinetic theory by Lun. Regarding the distribution function f in (3.20), different temperatures T t and T r are assigned to transitional and rotational degrees of freedom. This is not a redundant setting since the violation of equipartition of energy between transitional and rotational degrees of freedom indeed occurs. For example, Huthmann and Zippelius [21] showed that the equipartition is immediately destroyed even if it is initially satisfied. Some studies (e.g., Goldhirsh, Noskowicz and Bar-Lev [22], Kranz et al. [23]) suggest that there are substantial deviations from Gaussian and correlation between C i and W ji in the distribution function f (1) . Therefore, there can be a substantial deviation from (3.19) with (3.20) and (3.23) in f (1) . Since our aim is to extract information for the constitutive equations, full detailed information on f (1) is not necessary. However, we do not know how much information is sufficient for our purpose a priori. Here, it is assumed that the effect of the corrections on moments higher than two in f (1) 0 is not so large. Note that it is suggested by a number of studies that β depends on the angle ϑ := arctan(ξ t /ξ n ) of the oblique collision of the particles, where ξ n = k i ξ i and ξ t = ξ i ξ i − ξ 2 n . Walton and Braun [24] where β 0 is the maximum value of the coefficient of roughness and ϑ 0 is the critical angle. The form of β is verified in the numerical simulation by Saitoh and Hayakawa [8]. Therefore, the constant β assumed in the kinetic theory by Lun, should be regarded as an approximation. Again, we do not know a priori how much detailed information on β is required to obtain constitutive equations. It should be verified whether the simplified model of constant β is adequate herein. In Kanatani's theory, a relatively large ν is considered. Every particle is in contact with some other particles during almost all the time. The estimate of the energy dissipation µf ξ per unit time per a contact point implies that the velocity difference ξ is maintained to the order of the macroscopic time scale. However, it is not clear whether this is the general case for a large ν. The velocity difference ξ possibly decays in the order of microscopic time scale during the contact due to friction. Furthermore, among many choices of constitutive equations that are consistent with a given dissipation function Φ(E ji , R ji , Ω kji ), a particular choice (4.12) is made. When the order of homogeneity ζ is 2, the choice of constitutive equations (4.12) leads to the linear response satisfying Onsager's reciprocal relations [25]. Thus, the choice is valid. However, when ζ = 2 and the constitutive equations are nonlinear, such a validation is not possible.
Some of the points given above will be examined by means of numerical simulations in the next section.

Setting
We performed numerical simulations of a system of granular particles using a distinct element method (DEM). In DEM, the interaction forces for every contacting pair of particles are calculated at each time step, and the equations of the motion are solved using a difference method. The dimension d of the configuration space is 2. All the particles are disks of the same radius a. The mass 13401-10 is uniformly distributed inside the particles and thus the mass m and the moment of inertia I of the particles are related as K = I/ma 2 = 1/2. As a mechanism of the contact, we assume Hooke's law of elasticity in the direction parallel to k of (3.1), that is, F = κl where l = max(2a−|x (β) −x (α) |, 0) is the overlap length and κ is the force constant. The kinetic friction force µF is applied at the contacting point in the direction to reduce the velocity difference ξ between the surfaces of the two particles (see section 4). We have chosen the above model of contact since it is one of the simplest models that induces both energy dissipation and rotational degree of freedom of particles, that are essential elements to consider the micropolar nature of the granular flow. For simplicity, we have neglected nonlinearity in the relation between l and F , inelasticity in the direction parallel to k and the static friction in the direction of the velocity difference ξ. By virtue of simplicity of the model, the microscopic characteristics of granular particles are completely determined by four parameters, the radius a, the mass m, the force constant κ and the kinetic friction coefficient µ.
N particles are put in a square domain with the length of sides L x and L y . The area fraction ν is given by ν = N πa 2 /L x L y . We fix L x /a = L y /a = 200 for all the simulations except for the case ν = 0.1 in which L x /a = L y /a = 600. The time scale associated with the collision is given by t col = (m/κ) 1/2 . The time step ∆t of DEM is set to 10 −2 t col , which is sufficiently smaller than t col so that the collision process is resolved in the simulation. Since there is no damping in the direction parallel to k, the time scale associated with the damping is t damp = ∞.
For the kinetic friction coefficient, we used µ = 0.3 and 0.8. Note that kinetic theories which eliminate the spin field as an independent field are suggested [6,7] for small kinetic friction coefficient. Therefore, the aspect of granular flows as a micropolar fluid would become significant for large µ. In order to investigate this aspect, we have chosen somewhat larger µ compared to those in some measurements, e.g. µ < 0.2 in [26].
We employ the Lees-Edwards periodic boundary conditions [27] for the velocity c of particles, whereγ is a constant. The external torqueτ ji =τ (δ jy δ ix − δ jx δ iy ), constant in time, is applied to every particle. The macroscopic fields ρ(x, y), v i (x, y) and ω ji (x, y) are defined by the spatial averages, where c (α) and w (α) are, respectively, the velocity and the spin of the particle labeled by α, N ∆ := α 1 and α denotes the summation over the particles such that x − ∆x/2 x (α) < x + ∆x/2 and y − ∆y/2 y (α) < y + ∆y/2. The notations E ji (x, y), R ji (x, y) and Ω kji (x, y) are introduced similarly to (2.18). The simulations were performed up to the time that the velocity field relaxes to a quasi-steady state. For all the sets of parameters that we investigate in this paper, the velocity field relaxed to a uniform simple shear profile with the shear rateγ, i.e., v i (y) = (γy + A)δ ix where A is a constant. Here, we have chosen ∆x = L x so that v i is independent of x. In terms of E ji (y), we have It was found in the simulations that, when E ji (y) relaxed to a uniform field, ω ji (y) also relaxed to a uniform field. Consequently, we have The shear rateγ is directly controlled by the boundary conditions (6.1) and (6.2). We observed R = 0, i.e., the spin ω ji is subordinate to ∂ [j v i] , when there is no external torqueτ . The value of R is controlled by changing the value ofτ . In the present study, we variedτ in the rangeτ 0. In such cases, we observed R 0, i.e., ω yx ∂ [y v x] , at the quasi-steady states. Let us consider a line x i = h where h is a constant. In the simulation, the kinetic contribution to the stress tensor, σ

13401-11
where the summation (h) α is taken over the particles α which cross the line x j = h during the time step ∆t. The contribution of contacts to the stress, σ where F α→β i is the contacting force applied to the particle β by the particle α and the summation α,β is taken over the contacting pairs {α, β} such that x j . The mean stress σ ji (h) over a line x j = h can be averaged over h to give the mean stress σ ji of the whole domain,  such a case, σ + and σ − do not depend on t col . Then, a dimensional consideration yields similarity forms, whereσ ± are dimensionless functions, which may depend on the kinetic friction coefficient µ and the area fraction ν. The origin of the function sgn(γ) is a requirement of the symmetry σ ± (−γ, −R) = −σ ± (γ, R). Note that, in the case of R = 0, (6.12) reduces to the well known Bagnold scaling [28], σ + ∝γ 2 . In figure 2, σ ± /mγ 2 from the simulations are plotted as a function of R/γ forγt col = 0.5 × 10 −3 , 1.0 × 10 −3 and 2.0 × 10 −3 (µ and ν are fixed). The collapse of the data implies that the similarity forms (6.12) are valid within the parameter range of the simulations. By virtue of (6.12), we can fixγ and vary only R to obtain the constitutive equations. In the following simulations, we putγt col = 1.0 × 10 −3 . The value of the radial distribution function at r = 2a, g 0 = g 0 (2a), in the simulations is given in figure 3(a). The data are given for R = 0, but they are almost independent of R. Note that, for the elastic system (e = 1) without friction (β = −1) and shearing (γ = 0), the low density expansion of g 0 at a thermal equilibrium state can be given up to the first order of ν as g 0 = 1+ν[(8/3)−( √ 3/π)] (see, for example, [29]), and it is also shown in a dashed line in figure 3 for comparison. One can see from figure 3 that g 0 of the simulation is approximated by the low density expansion for relatively low area fraction ν = 0.1, but it largely deviates from the low density expansion when ν 0.6. It is known that, when ν exceeds a critical value ν J , which is called a jamming point, the system enters the jammed phase, in which the shear stress remains nonzero in the limit of zero strain, i.e., yields stress. When ν approaches ν J , quantities such as the pressure p and g 0 are expected to obey certain scaling laws. It is found that the present system suffers the phase separation between crystallized region and fluid region for ν > 0.815, before entering the jammed phase. Since it is impossible to estimate ν J in the present system, we borrowed the value ν J = 0.84, which is a round-off value of ν J reported in [14] for two-dimensional poly-disperse frictionless granular system, as a reference. In figure 3(b), p and g 0 are plotted as functions of |δν| = |ν − ν J | for various ν and R = 0. Since the crystallized phase is out of the scope of the present paper, data for ν > 0.815 are omitted in the figure. The scaling laws p ∼ |δν| −4 and g 0 ∼ |δν| 4 predicted by Otsuki and Hayakawa [14] for the particles with frictionless, linear spring interactions, are indicated by lines in the figures as a reference. Since there are linear spring and frictional interactions in the present model, the scaling law could be modified. However, since the present system experiences crystallization, the scaling range is very narrow, in case it exists. For each particle-pair collision, we can define the interfering number n int , which is the maximum number of other particles contacting with the particle-pair simultaneously during the collisions of the particle-pair. n int = 0 implies that the particle-pair collision is not interfered by the other particles. In figure 4, the probability P (n int ) is given for ν = 0.1, 0.7 and 0.8 in the case of µ = 0.3. The data for µ = 0.8 are omitted since they almost collapse with those of µ = 0.3. The fact that P (n int ) has a peak at n int = 0 for ν = 0.1 and 0.7 is consistent with the kinetic theory. However, the peak shifts to n int = 2 for ν = 0.8. This suggests the failure of the applicability of the kinetic theory based on binary collisions. Note that the coordination number Z at ν = 0.8 is about 0.6−0.7 and it is still much smaller than the value Z = d + 1 = 3 for the isostatic condition [30].

13401-13
Summarizing the above observations, we can consider that 0.7 ν 0.8 is located in an intermediate regime where ν is neither low enough so that the kinetic theory is applicable nor it is high enough so that the system is quite near the jamming point ν J .
The above conclusion might seem to be inconsistent with the results of the numerical simulations of frictional granular shear flows by Otsuki and Hayakawa [31]. According to [31], when friction is introduced to the system, the critical point ν J of the area fraction splits into two points ν S and ν L , where the former belongs to the solid branch and the latter liquid branch. The point ν = 0.80 for µ = 0.8 is located in the region ν S < ν < ν L , which implies that the selection of the solid or liquid phase depends onγt col . Although the figures forγt col = 10 −3 , which is the value in the present study, are not given in [31], it can be located in the solid (jammed) phase. Such a discrepancy between the results of the present simulations and those in [31] can occur since the details of the system are different in the present study and [31]. For example, the radius of particles is unique in the present study whereas particles with four different radii are present in [31]; the normal viscosity is not introduced in the present study whereas the normal viscous constant is presumably non-zero (although the actual value is not given) in [31]; and contacting particles are always slipping in the tangential direction in the present study whereas contacting particles can stick when the normal contact force is strong enough in [31]. Especially, the latter fact that the particles are always slipping in the present study would act to raise the value of ν J (ν S or ν L ) from those in [31].

Comparison with the kinetic theory
In order to compare the simulation results with the kinetic theory by Lun, we need to determine the coefficient of roughness β in the simulations. In the simulations, we define β by  whereǩ is the unit vector parallel toř (2) −ř (1) . Since a collision has a finite duration in the simulations, k andǩ are distinguished here. The particle-pair possibly contacts with other particles during their own collision, especially when the area fraction ν is high and the duration of the collision is long.β(ϑ) is defined as the average of β over the particle pairs with the oblique collision angle ϑ. In figure 5,β(ϑ) is plotted with the probability density function P (ϑ) of ϑ. The data is for the simulation without the external torque and the dependency ofβ(ϑ) on the external torque is weak (the figure omitted). For all (ν, µ) we investigated, we have −1 β (ϑ) 1 which is consistent with the assumption −1 β 1 in the kinetic theory. When ν = 0.1,β is in agreement with the theoretical estimate (5.1) with β 0 = 0. As ν increases,β(ϑ) deviates from the estimate (5.1). The deviation is in the direction of decreasing the ϑ dependence ofβ(ϑ). This deviation can be understood as the effect of the collision of the particle-pair with the other particles during the collision within the particle-pair. The other particle randomly changes the velocity differenceξ just after the particle-pair collision. Therefore, β becomes a random variable for fixed ϑ. For reference, the standard deviations ∆β at ϑ = 1 are 0.004(ν = 0.1), 0.05(ν = 0.7), 0.31(ν = 0.8) for µ = 0.3 and 0.0008(ν = 0.1), 0.05(ν = 0.7), 0.14(ν = 0.8) for µ = 0.8. As expected, ∆β increases with ν.
In figure 5, the dashed line shows cos ϑ, the form for the random collisions without spin. Particles tend to collide with large angle ϑ when spin and frictional interaction are introduced. This can be understood by a simple model of a particle-pair with the same spin w colliding with each other  with the relative velocity c. Let b be the impact parameter and Θ(−π/2 Θ π/2) the angle between c and k. Then, we have b = 2a sin Θ, tan ϑ = tan Θ − ς cos Θ , (6.14) where ς := 2aw/c. Assuming the uniform distribution of b in −a b a, we have for the probability densityP (ϑ) of ϑ (−π/2 ϑ π/2) as In figure 5, the symmetrized probability density function P (ϑ) :=P (ϑ) +P (−ϑ)(0 ϑ π/2) with ς = 0.5 is given in a dot-dashed line for reference. The model probability density function roughly fits the corresponding simulation results for ν = 0.1 and 0.7. In this simple model, c and w are represented by constants and their probability density functions are not considered. If we take c and w as the root mean square of the fluctuations in the simulations, we have ς ∼ 0.7 and it is not far from ς = 0.5. Thus, P (ϑ) for ν 0.7 can be roughly explained by the effect of spin. However, it is found that P (ϑ) of ν = 0.8 is hard to be fitted by that of the model. Especially, the bump around ϑ = 0 is not explained by the model. It seems that n-particles interaction with n > 2 should be taken into account in order to understand P (ϑ) for ν = 0.8.
In the kinetic theory by Lun, β is treated as a constant. Whether we can, or cannot, approximate ϑ dependentβ(θ) by a constant is not obvious. Here, we take a simple average, β := π/2 0 dϑβ(ϑ)P (ϑ). The values ofβ for various ν and µ are given in figure 6. The stresses σ Here, the values obtained in the simulations are used for the "temperatures" T t and T r , and the value of the radial distribution function at r = 2a, g 0 = g 0 (2a). The coefficient of restitution e is set to 1 following the situation of the simulation. β is set to beβ.
When the area fraction is as small as ν = 0.1, the kinetic theory is in good agreement with the simulation results with regard to σ + . The agreement between the kinetic theory and the simulation for σ − is not as good as that of σ + . The relative deviations are about 40% for µ = 0.3 and 20% for µ = 0.8. This may be due to the approximation of the constant β in the theory. A better agreement for µ = 0.8 compared to that for µ = 0.3 can be explained by the fact that the range of ϑ in whichβ(ϑ) can be approximated by a constant is wider for µ = 0.8.
As ν increases, the ratio σ − . As expected from the fact that the peak of P (n int ) is shifted from n int = 0 for ν = 0.8 (see figure 4), the kinetic theory is no more appropriate for ν = 0.8 and the formal application yields the overestimates of σ (c) + by about 40% or higher. For σ − , the overestimate is larger by the factor of more than 2 (the figures omitted). Now, let us find the effective β eff that fits best with the constitutive relations obtained in the simulation. For every set of (µ, ν), we chose β = β eff where β eff is the value of β which minimizes the summation of squared relative errors, i.e., with the superscripts 'th' and 'sim' denoting the kinetic theory and the simulation respectively, and the subscript i the index of simulations with different values of R. We excluded i with R i = 0 in the summation for σ − since σ − (R = 0) is ideally 0 and thus it is not appropriate to use the relative error. In the vicinity of β = β eff , the function C(β) can be approximated by a quadratic function of β, i.e, C(β) ≃C(β) = α(β − β eff ) 2 + C min . We determined the error ∆β of β to satisfỹ C(β ± ∆β) = 2C min .
The obtained values of β eff and ∆β are given in figure 6 together withβ. Stresses σ (k) + , σ (c) + and σ − obtained from (3.26), (3.27) and (3.28) with β = β eff are given in figures 7 and 8. One can see from the figures that σ − can be well fitted to the simulation data by adjusting β. On the other hand, σ + is less sensitive to β and so the deviation from the simulation data is not effectively reduced by adjusting β. The values of β eff are considerably smaller thanβ for all ν and µ that we have investigated. The discrepancy betweenβ and β eff becomes larger as ν becomes larger.

Comparison with Kanatani's theory
As we have confirmed in the previous subsection, the estimates from the kinetic theory are in fairly good agreement with the simulation results up to ν = 0.7. When ν = 0.8, the peak of P (n int ) is not n int = 0, that is, most pair collisions are interfered by the other particles. Consequently, the kinetic theory is not applicable. In this subsection, we examine the applicability of Kanatani's theory for relatively dense granular flows to the results of the simulation for ν = 0.8.   Under the conditions (6.4) and (6.5), the dissipation function Φ of (2.17) reduces to Φ(γ, R) = σ + (γ, R)γ + 2σ − (γ, R)R.

13401-17
(6.18) In figure 9, Φ obtained from the simulation is plotted with the estimate (4.10) of Φ in Kanatani's theory. In the estimate, the values from the simulation are used for p c . One can see from the figure that they are not in agreement. Note that the estimate (4.10) depends linearly on µ besides the possible µ dependence in p c . However, Φ from the simulations depends on µ less sensitively. It is suggested from the disagreement of the µ dependence that the process of energy dissipation assumed in Kanatani's theory is not valid in the present situation of the simulation. Although, at ν = 0.8, the effect of a multi-contact of the particle may be significant, the velocity difference at the contact point would not be completely sustained during the contact. Now, let us examine the choice (4.12) of the constitutive equations in Kanatani's theory. In order to focus on the examination of the choice (4.12), let us use Φ obtained in the simulation rather than Φ of the estimates in Kanatani's theory. The form (6.18) with (6.12) implies that the degree of homogeneity ζ is 3. The choice (4.12) reads, in the present context, as follows: by fitting the function of the form (6.23) directly to the simulation data. One can see that they disagree. Especially, in the decomposition based on (6.19), we have B + = A − according to (6.24). However, this is not the case in the simulation. The disagreement can be also checked in figure 10; the normalized stressesσ ± = σ ± /mγ 2 of (6.23) with (6.25) and (6.26) are plotted together with the normalized stresses directly measured in the numerical simulation and their fits (6.23) with (6.27) and (6.28). Thus, the choice (4.12) of the constitutive equations made by Kanatani is not appropriate in the present situation of the simulation.

Discussion
We confirmed in section 6.2 that the constitutive equation for σ + in the kinetic theory by Lun is in good agreement with the simulation results for a small area fraction ν = 0.1. But the agreement for σ − is not as good as σ + . Since σ − is sensitive to β, this may be due to the oversimplification of modeling β by a constant in the Lun's theory. As ν increases, discrepancy between the kinetic theory and the simulation results on both σ + and σ − increases. However, σ − in the kinetic theory can be formally fitted to the simulation results by introducing the effective coefficient of roughness β eff which is substantially smaller than the actual averaged betaβ in the simulation. This implies that, although the particle-pairs almost get stuck (β ∼ 0) in the tangential direction during the collision, the macroscopic field feels as if the particles were slipping (β eff ∼ −1). A possible scenario is that the collective motions of some stuck particles affect the macroscopic field as motions of virtual particles with renormalized mass m, radius a, coefficient of restitution e, and β. In the present study, we renormalized β with fixed m, a and e. Then β is renormalized to reduce its value.
We have seen in section 6.3 that Kanatani's theory is not applicable to the present situation of simulation (ν = 0.8) in the following sense; (i) the kinetic friction coefficient µ dependence of the dissipation function Φ and (ii) the choice (4.12) made for determining the constitutive equations from Φ, are in disagreement with those in the simulations.
To summarize, there is a regime of relatively dense (ν ∼ 0.8 in the present simulation) granular flows in the form of micropolar fluid, where physical pictures based on neither kinetic theories nor Kanatani's theory are adequate. In this regime, a new microscopic description of particle interaction, that is different from mutually independent short-time collisions in the kinetic theory or long-time contacts with sustained velocity difference in Kanatani's theory, is necessary. As we have seen in sections 6.1 and 6.2, there are considerable changes in the interfering number n int ,β(ϑ) and P (ϑ) when ν increases from 0.7 to 0.8. These changes suggest the significance of the effect of n-particle interactions with n > 2 in this regime. It would be a future study to see whether the effect is included in the kinetic theory with renormalized parameters as we have preliminary analyzed or in some variant of Kanatani's theory with an appropriate dissipation function and its decomposition.