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 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 by renormalizing the coefficient of roughness. For a relatively dense granular flow (), 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.


granular flow, constitutive equation, micropolar fluid, kinetic equation \pacs45.70.Mg, 47.57.Gc, 47.50.-d


Для того, щоб виокремити матерiальнi рiвняння для системи були здiйсненi числовi симуляцiї двовимiрних гранульованих потокiв пiд дiєю однорiдного зсуву i зовнiшнього крутильного моменту. Результат чисельних симуляцiй проаналiзовано на основi моделi мiкрополярного плину. В симуляцiях реалiзується поле однорiдного середнього зсуву, яке не є пiдпорядковане полю вихоровостi. Оцiнки напружеь, зробленi на основi кiнетичної теорiї Люна [Lun, J. Fluid Mech. 233(1991) 539], добре узгоджуються з результатами симуляцiй в областi низьких часток , але узгодження погiршується, якщо ця величина зростає. Проте, оцiнки, зробленi в кiнетичнiй теорiї можуть бути пiдiгнанi до результатiв симуляцiй аж до шляхом ренормалiзацiї коефiцiнта шорсткостi. Для вiдносно густого гранульованого потоку (), результати симуляцiй також порiвнюються з теорiєю Канатанi [Kanatani, Int. J. Eng. Sci 17(1979) 419]. Знайдено, що дисипативна функцiя i її декомпозицiя в матерiальнi рiвняння в теорiї Канатанi не узгоджується з результатами симуляцiй. \keywordsгранульований потiк, матерiальнi рiвняння, мiкрополярний плин, кiнетичне рiвняння


201114113401 \doinumber10.5488/CMP.14.13401 Constitutive equations for granular flow]Constitutive equations for granular flow with uniform mean shear and spin fields \addressGraduate School of Pure and Applied Sciences, University of Tsukuba,
1–1–1 Tennodai, Ibaraki 305–8571, Japan

1 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 -particles with 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 , and the spatial derivative of the spin field, which will be denoted by , 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 and . Note that 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. ) with are approximately realized in the region. The situation and 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 and the magnitude of can be controlled (see section 6). Thus, we concentrate on the dependence of the constitutive equations with fixed to . The study of the 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.

2 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 and the same moment of inertia . Let and be, respectively, the velocity, the spin and the position of the particle at time where and are coordinate indices of -dimensional space. Here, can formally be an arbitrary positive integer larger than . In this paper, we use the convention that the spin is expressed by a skew-symmetric tensor whose -th component gives the angular velocity in the coordinate plane. Let be the probability density function in the phase space of -particles system satisfying Liouville equation. Here, the bold letters , and denote vector or tensor, the superscript on and is the index of the particle and is symmetrized with respect to interchanges of the particles. The -particles set distribution function is given by


and the number density of the -particles sets is given by


Macroscopic fields such as the mass density field , the moment of inertia density field , the velocity field and the spin field are introduced as




for an arbitrary function of and . 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 , is the stress tensor, the couple-stress tensor, the external body force and the external body torque. Here and hereafter, summations are taken over repeated subscripts and we employ the notations,


for an arbitrary tensor . Equations (6)–(9) correspond, respectively, to the conservation laws of mass, moment of inertia, momentum and angular momentum. Note that and are mutually independent degrees of freedom. When there is no spin field , (6)–(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 and are known in terms of and . The equations which yield such relations are called constitutive equations.

Let and be the fluctuations of velocity and spin of individual particles around the macroscopic fields, i.e.,


The ‘‘internal energies’’ per unit mass and associated with the fluctuations and , respectively, are given by


Here, the subscripts t and r indicate the translational and rotational motions, respectively. Two kinds of ‘‘temperature’’, and , are introduced by the relations and . The total ‘‘internal energy’’ is given by . 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 (6)–(9) that the kinetic energy of macroscopic fields per unit mass,






where we have introduced notations,


and is the pressure. Note that is the work done on a fluid element per unit volume per unit time by the stress , the couple stress , the external body force and the external body torque . From (15) and the energy budget equation,


we obtain


where is input of the ‘‘internal energy’’ per unit volume per unit time, is the ‘‘internal energy’’ flux. The quantity gives the energy transferred from the kinetic energy to the ‘‘internal energy’’ per unit volume per unit time. When macroscopic fields and are constant in time and space, we have from (6), and (20) reduces to


The equation (21) implies that the energy transfer rate from the kinetic energy to the ‘‘internal energy’’ balances with the dissipation rate of ‘‘internal energy’’ for macroscopically steady states. is called the dissipation function since it gives the energy going out from the kinetic energy per unit time per unit volume in the macroscopically steady state.

3 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.

Figure 1: Configuration of two contacting granular particles.

The dimension of the configuration space is and the particles are assumed to be -dimensional spheres with the same radius . The mass and the moment of inertia are, respectively denoted by and as in section 2. Let the colliding two particles be labeled as and , and a unit vector be defined by


(see figure 1). In this section, we employ the notation that the quantities without (with) check  denote those just before (after) the contact. Let be the impulse of the force exerted on the particle by the particle through the contact point. We have


Let us assume that the velocity difference between the two surfaces at the contact point,




changes as,


after the collision. Here, the coefficient of restitution and the coefficient of roughness are assumed to be constants satisfying and . From (23), (24), (27), (28), we obtain


where , and .

Assuming the binary collisions, the equation of motion for , where is an arbitrary function of and , can be written as




By substituting and in (31), one finds that the stress tensor can be written as




where denotes the kinetic contribution to the stress due to the particles that cross the plane perpendicular to -axis, and denotes the collisional contribution due to the collisions of the two particles in different sides of the plane.

The particle-pair distribution function is assumed to be approximated by the form


where is a radial distribution given by


with an assumption that it is insensitive to the direction of a unit vector . Assuming that is sufficiently smaller than the typical spatial scale in which the amplitude of varies, we have


Let be written as


where is the distribution function at a local equilibrium state and represents the perturbation. The function is given by


Note that different temperatures and 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 , i.e., , 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 , and the translational and rotational kinetic energy fluxes and given by


In what follows, we restrict ourselves to the case of macroscopic fields and being constant in space and time. In such a case, there are no energy fluxes and and the perturbation depends solely on . Let us assume the form


where the factor is determined from the consistency condition that (35) is satisfied. Note that the total is given by


where there is no anti-symmetric part in the definition (35). By substituting into (31), we obtain


where and are given by (35) and (36), respectively, and by (33). and reduce to functions of , and , after performing the integrations with respect to and in (32) and (33). Substitute the functions into (46) and assume that the magnitudes of and are so small that the second or higher order terms in and can be neglected. Then we arrive at the equations,




and is the volume fraction which may be related to the number density of the particle as


Now, the constitutive equation for the total stress tensor is given by


with (47)–(49). From symmetry, we have  1. The obtained constitutive equations are a restricted version of those in [10] in the sense that constants , , and are assumed, and they are a generalized version in the sense that is arbitrary whereas in [10].

4 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 , mass and moment of inertia is considered. Let the two contacting spherical particles be labeled as 1 and 2. The unit vector is the same as (22) (see figure 1). From (12), the quantities associated with particles, and , can be related to the macroscopic fields as


where and denote the fluctuations around the macroscopic fields. Since is small compared to the typical scale in which the amplitudes of fields vary, and can be approximated by its Taylor series about up to the first-order terms, i.e.,


where we have omitted writing the argument in and . From (25), (53) and (54), the velocity difference between the two surfaces at the contacting point is now written in terms of , , , , and . Assume that the fluctuations and can be neglected and that is isotropically distributed. Then, we have




and denote the average over . Here, we have used


which can be derived using the general expressions of isotropic tensors and . Suppose that the average energy dissipation per unit time at a contacting point is estimated by , where is the kinetic friction coefficient and is the average amplitude of the force applied in the direction at the contacting point. Let be a number of contacting points per a particle. Then, the energy dissipation per unit volume per unit time is given by


where the factor is introduced to cancel the double counting of the contacting points. Note that the pressure due to the contacting force can be estimated by


where is the volume fraction and is the surface area of the particle. From (55), (59) and (60), one arrives at


The estimate of is given as a homogeneous function of , and , i.e., 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 (21) and (61), one obtains


which is a homogeneous function of , and , i.e., with . With the help of Euler’s homogeneous function theorem,


the choice


is made for constitutive equations which are consistent with (17).

5 Comments on the kinetic theory and Kanatani’s theory

In the kinetic theory by Lun, the binary collision is assumed, i.e., -particle collisions for are neglected, and the particle-pair distribution is assumed to be approximated by a product of and as (37). These assumptions seem to be valid for small volume fraction, . This is because, in such a case, -particle collisions for would be rare and the velocities of the two particles before a binary collision would be statistically almost independent. Furthermore, macroscopic fields and 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 in (41), different temperatures and 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 and in the distribution function . Therefore, there can be a substantial deviation from (40) with (41) and (44) in . Since our aim is to extract information for the constitutive equations, full detailed information on 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 is not so large.

Note that it is suggested by a number of studies that depends on the angle of the oblique collision of the particles, where and . Walton and Braun [24] derived


where is the maximum value of the coefficient of roughness and 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 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 , a particular choice (64) is made. When the order of homogeneity is , the choice of constitutive equations (64) leads to the linear response satisfying Onsager’s reciprocal relations [25]. Thus, the choice is valid. However, when 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.

6 Numerical simulations

6.1 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 of the configuration space is . All the particles are disks of the same radius . The mass is uniformly distributed inside the particles and thus the mass and the moment of inertia of the particles are related as . As a mechanism of the contact, we assume Hooke’s law of elasticity in the direction parallel to of (22), that is, where is the overlap length and is the force constant. The kinetic friction force 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 and , inelasticity in the direction parallel to 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 , the mass , the force constant and the kinetic friction coefficient .

particles are put in a square domain with the length of sides and . The area fraction is given by . We fix for all the simulations except for the case in which . The time scale associated with the collision is given by