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

Повний опис

Збережено в:

Бібліографічні деталі
Дата:	2011
Автори:	Takechi, K., Yoshida, K., Arimitsu, T.
Формат:	Стаття
Мова:	English
Опубліковано:	Інститут фізики конденсованих систем НАН України 2011
Назва видання:	Condensed Matter Physics
Онлайн доступ:	http://dspace.nbuv.gov.ua/handle/123456789/119938
Теги:	Додати тег Немає тегів, Будьте першим, хто поставить тег для цього запису!
Назва журналу:	Digital Library of Periodicals of National Academy of Sciences of Ukraine
Цитувати:	Constitutive equations for granular flow with uniform mean shear and spin fields / K. Takechi, K. Yoshida, T. Arimitsu // Condensed Matter Physics. — 2011. — Т. 14, № 1. — С.13401: 1-22. — Бібліогр.: 31 назв. — англ.

Репозитарії

Digital Library of Periodicals of National Academy of Sciences of Ukraine

id	irk-123456789-119938
record_format	dspace
spelling	irk-123456789-1199382017-06-11T03:04:38Z Constitutive equations for granular flow with uniform mean shear and spin fields Takechi, K. Yoshida, K. Arimitsu, T. 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 ν = 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 ν = 0.7 by renormalizing the coefficient of roughness. For a relatively dense granular flow (ν= 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. Для того, щоб виокремити матер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 низьких часток ν = 0.1, але узгодження погiршується, якщо ця величина зростає. Проте, оцiнки, зробленiв кiнетичнiй теорiї можуть бути пiдiгнанiдо результатiв симуляцiй аж до ν = 0.7 шляхом ренормалiзацiї коефiцiнта шорсткостi. Для вiдносно густого гранульованого потоку(ν= 0.8), результати симуляцiй також порiвнюються з теорiєю Канатанi [Kanatani, Int. J. Eng. Sci17(1979) 419]. Знайдено, що дисипативна функцiя i її декомпозицiя в матерiальнi рiвняння в теорiї Канатанiне узгоджується з результатами симуляцiй. 2011 Article Constitutive equations for granular flow with uniform mean shear and spin fields / K. Takechi, K. Yoshida, T. Arimitsu // Condensed Matter Physics. — 2011. — Т. 14, № 1. — С.13401: 1-22. — Бібліогр.: 31 назв. — англ. 1607-324X PACS: 45.70.Mg, 47.57.Gc, 47.50.-d DOI:10.5488/CMP.14.13401 arXiv:1106.5306 http://dspace.nbuv.gov.ua/handle/123456789/119938 en Condensed Matter Physics Інститут фізики конденсованих систем НАН України
institution	Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection	DSpace DC
language	English
description	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 ν = 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 ν = 0.7 by renormalizing the coefficient of roughness. For a relatively dense granular flow (ν= 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.
format	Article
author	Takechi, K. Yoshida, K. Arimitsu, T.
spellingShingle	Takechi, K. Yoshida, K. Arimitsu, T. Constitutive equations for granular flow with uniform mean shear and spin fields Condensed Matter Physics
author_facet	Takechi, K. Yoshida, K. Arimitsu, T.
author_sort	Takechi, K.
title	Constitutive equations for granular flow with uniform mean shear and spin fields
title_short	Constitutive equations for granular flow with uniform mean shear and spin fields
title_full	Constitutive equations for granular flow with uniform mean shear and spin fields
title_fullStr	Constitutive equations for granular flow with uniform mean shear and spin fields
title_full_unstemmed	Constitutive equations for granular flow with uniform mean shear and spin fields
title_sort	constitutive equations for granular flow with uniform mean shear and spin fields
publisher	Інститут фізики конденсованих систем НАН України
publishDate	2011
url	http://dspace.nbuv.gov.ua/handle/123456789/119938
citation_txt	Constitutive equations for granular flow with uniform mean shear and spin fields / K. Takechi, K. Yoshida, T. Arimitsu // Condensed Matter Physics. — 2011. — Т. 14, № 1. — С.13401: 1-22. — Бібліогр.: 31 назв. — англ.
series	Condensed Matter Physics
work_keys_str_mv	AT takechik constitutiveequationsforgranularflowwithuniformmeanshearandspinfields AT yoshidak constitutiveequationsforgranularflowwithuniformmeanshearandspinfields AT arimitsut constitutiveequationsforgranularflowwithuniformmeanshearandspinfields
first_indexed	2025-07-08T16:56:32Z
last_indexed	2025-07-08T16:56:32Z
_version_	1837098644517421056
fulltext	Condensed Matter Physics 2011, Vol. 14, No 1, 13401: 1–22 DOI:10.5488/CMP.14.13401 http://www.icmp.lviv.ua/journal Constitutive equations for granular flow with uniform mean shear and spin fields K. Takechi, K. Yoshida, T. Arimitsu Graduate School of Pure and Applied Sciences, University of Tsukuba, 1–1–1 Tennodai, Ibaraki 305–8571, Japan Received May 12, 2010, in final form December 22, 2010 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 simu- lations 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 ν = 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 ν = 0.7 by renormaliz- ing the coefficient of roughness. For a relatively dense granular flow (ν = 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. Key words: granular flow, constitutive equation, micropolar fluid, kinetic equation PACS: 45.70.Mg, 47.57.Gc, 47.50.-d 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 c© K. Takechi, K. Yoshida, T. Arimitsu 13401-1 http://dx.doi.org/10.5488/CMP.14.13401 http://www.icmp.lviv.ua/journal K. Takechi, K. Yoshida, T. Arimitsu 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 re- gion 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 micropo- lar 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 Rji, 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 Rji 6= 0 and Ωkji = 0. Note that Rji 6= 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 Rji 6= 0 are approximately realized in the region. The situation Rji 6= 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 13401-2 Constitutive equations for granular flow 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 Rji can be controlled (see section 6). Thus, we concentrate on the Rji 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. 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 m and the same moment of inertia I. Let ci(t), wji(t) and ri(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 wji 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 6 N) is given by f (s)(c(1),w(1), r(1); · · · ; c(s),w(s), r(s); t) = N ! (N − s)! N ∏ α=s+1 ( ∫ dc(α) ∫ dw(α) ∫ dr(α) ) × F (N)(c(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 n(s)(r(1), · · · , r(s); t) = s ∏ α=1 ( ∫ dc(α) ∫ dw(α) ) × f (s)(c(1),w(1), r(1); · · · ; c(s),w(s), r(s); t). (2.2) 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), (2.3) v(r, t) := 〈c〉r,t , ω(r, t) := 〈w〉r,t , (2.4) where 〈ψ(c,w)〉r,t := 1 n(1)(r, t) ∫ dc ∫ dwψ(c,w)f (1)(c,w, r; t), (2.5) for an arbitrary function ψ of c andw. Hereafter, indices or subscripts of spatial or time coordinates will be suppressed unless we need to emphasize them. 13401-3 K. Takechi, K. Yoshida, T. Arimitsu These macroscopic fields satisfy the following equations, Dρ Dt + ρ∂ivi = 0, (2.6) DρI Dt + ρI∂ivi = 0, (2.7) ρ Dvi Dt = ∂jσji + ρbi , (2.8) ρI Dωji Dt = 2σ[ji] + ∂kλkji + ρIτji , (2.9) whereD/Dt = ∂/∂t+vi∂i, σji is the stress tensor, λkji the couple-stress tensor, bi the external body force and τji the external body torque. Here and hereafter, summations are taken over repeated subscripts and we employ the notations, T(ji) = 1 2 (Tji + Tij) , T[ji] = 1 2 (Tji − Tij) , (2.10) T{ji} = T(ji) − 1 d Tkkδji , (2.11) for an arbitrary tensor Tji. Equations (2.6)–(2.9) correspond, respectively, to the conservation laws of mass, moment of inertia, momentum and angular momentum. Note that vi 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 ∂jvi 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., C = c− v, W = w − ω. (2.12) The “internal energies” per unit mass ǫt and ǫr associated with the fluctuations C and W , respec- tively, are given by ǫt := 1 2 〈 CiCi 〉 , ǫr := ρI 4ρ 〈 WjiWji 〉 . (2.13) Here, the subscripts t and r indicate the translational and rotational motions, respectively. Two kinds of “temperature”, Tt and Tr, are introduced by the relations ǫt = (d/2)Tt and ǫr = (d(d − 1)/4)Tr. 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, ǫK = 1 2 vivi + ρI 4ρ ωjiωji , (2.14) satisfies ρ DǫK Dt = Ψ+ p∂ivi − Φ, (2.15) with Ψ := ∂k (σkivi) + 1 2 ∂k (λkjiωji) + ρbivi + 1 2 ρIτjiωji , (2.16) Φ := σ{ji}Eji + σ[ji]Rji + 1 2 λkjiΩkji , (2.17) where we have introduced notations, Eji := ∂{jvi}, Rji := ∂[jvi] − ωji , (2.18) 13401-4 Constitutive equations for granular flow Figure 1. Configuration of two contacting granular particles. 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 bi and the external body torque τji. From (2.15) and the energy budget equation, ρ D Dt (ǫK + ǫU) = Ψ + q − ∂ihi , (2.19) we obtain ρ DǫU Dt = Φ− p∂ivi + q − ∂ihi , (2.20) where q is input of the “internal energy” per unit volume per unit time, hj is the “internal energy” flux. The quantity Φ− p∂ivi gives the energy transferred from the kinetic energy to the “internal energy” per unit volume per unit time. When macroscopic fields ρ, ∂jvi, ǫU and hj are constant in time and space, we have ∂ivi = 0 from (2.6), and (2.20) reduces to Φ = −q. (2.21) 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. 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. 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 k := r(2) − r(1) \|r(2) − r(1)\| , (3.1) (see figure 1). In this section, we employ the notation that the quantities without (with) check ˇ denote 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 mč (1) i = mc (1) i − Ji , mč (2) i = mc (2) i + Ji , (3.2) Iw̌ (1) ji = Iw (1) ji − a(kjJi − kiJj), Iw̌ (2) ji = Iw (2) ji − a(kjJi − kiJj). (3.3) 13401-5 K. Takechi, K. Yoshida, T. Arimitsu Let us assume that the velocity difference ξi between the two surfaces at the contact point, ξi := (c−l + akjw + jl)(δli − klki), (3.4) where c± := c(1) ± c(2), w± := w(1) ± ω(2), (3.5) changes as, kiξ̌i = −ekiξi , (3.6) (δij − kikj)ξ̌j = −β(δij − kikj)ξj , (3.7) after the collision. Here, the coefficient of restitution e and the coefficient of roughness β are assumed to be constants satisfying 0 6 e 6 1 and −1 6 β 6 1. From (3.2), (3.3), (3.6), (3.7), we obtain č (1) i − c (1) i = −č(2)i + c (2) i = −η2v−i − (η1 − η2)kjc − j ki − η2akjw + ji , (3.8) w̌ (1) ji − w (1) ji = w̌ (2) ji − w (2) ji = − η2 Ka [ (kjc − i − kic − j + a(klkjw + li − klkiw + lj) ] , (3.9) where η1 = (1 + e)/2, η2 = (1 + β)K/2(K + 1) and K = I/ma2. Assuming the binary collisions, the equation of motion for 〈ψ〉, where ψ is an arbitrary function of c and w, can be written as ∂ ∂t ( n(1)〈ψ〉 ) = −∂i ( n(1)〈ciψ〉 ) − ∂iθi(ψ) + χ(ψ) + n(1) (〈 ∂ψ ∂ci 〉 bi + 1 2 〈 ∂ψ ∂wji 〉 τji ) , (3.10) where θi(ψ; r, t) :=− 1 2 ∫ dc(1) ∫ dc(2) ∫ dw(1) ∫ dw(2) ∫ k·c−>0 dk × k · c−(2a)dki(ψ̌(1) − ψ(1)) × f (2) ( c(1),w(1), r − ak; c(2),w(2), r + ak; t ) , (3.11) χ(ψ; r, t) := 1 2 ∫ dc(1) ∫ dc(2) ∫ dw(1) ∫ dw(2) ∫ k·c−>0 dk × k · c−(2a)d−1(ψ̌(1) − ψ(1) + ψ̌(2) − ψ(2)) × f (2) ( c(1),w(1), r − ak; c(2),w(2), r + ak; t ) . (3.12) By substituting ψ = m and ψ = mci in (3.10), one finds that the stress tensor σji can be written as σji = σ (k) ji + σ (c) ji , (3.13) with σ (k) ji :=− ρ〈CjCi〉, (3.14) σ (c) ji :=− θj(mCi), (3.15) where σ (k) 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 f (2)(c(1),w(1), r − ak, c(2),w(2), r + ak; t) ≃ g0(2a; r, t)f (1)(c(1),w(1), r − ak; t)f (1)(c(2),w(2), r + ak; t), (3.16) 13401-6 Constitutive equations for granular flow where g0(r ′; r, t) is a radial distribution given by g0(r ′; r, t) := n(2)(r − r′ 2 e, r + r′ 2 e; t) [n(1)(r, t)]2 , (3.17) 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 f (1)(c,w, r ± ak; t) ≃ f (1)(c,w; r, t) [ 1± ak ·∇ ln f (1)(c,w, r; t) ] . (3.18) Let f (1) be written as f (1)(c,w, r; t) = f (1) 0 (c,w, r; t)[1 + φ(c,w, r; t)], (3.19) where f (1) 0 is the distribution function at a local equilibrium state and φ represents the perturbation. The function f (1) 0 is given by f (1) 0 (c,w, r; t) = n(1)(r, t) (2πTt(r, t)) d 2 (2πmTr(r, t)/I) d(d−1) 4 × exp [ −Ci(r, t)Ci(r, t) 2Tt(r, t) − IWji(r, t)Wji(r, t) 4mTr(r, t) ] . (3.20) Note that different temperatures Tt and Tr 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 wji, i.e., \|ωji\| ≪ √ (m/I)Tr, we have exp { −IWji(r, t)Wji(r, t) 4mTr(r, t) } ≃ [ 1 + Iwjiωji(r, t) 2mTr(r, t) ] exp [ − Iwjiwji 4mTr(r, t) ] . (3.21) 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 σ (k) {ij}, and the translational and rotational kinetic energy fluxes h (k) t,i and h (k) r,i given by h (k) t,i (r, t) := ρ 2 〈CiCjCj〉r,t , h (k) r,i (r, t) := ρI 4 〈CiWljWlj〉r,t . (3.22) In what follows, we restrict ourselves to the case of macroscopic fields ∂jvi, ωji, Tt and Tr being constant in space and time. In such a case, there are no energy fluxes h (k) t,i and h (k) r,i and the perturbation φ depends solely on σ (k) {ij}. Let us assume the form φ = − 1 2ρTt σ (k) {ij}CiCj , (3.23) where the factor −1/2ρTt is determined from the consistency condition that (3.14) is satisfied. Note that the total σ(k) is given by σ (k) ji = −ρTtδji + σ (k) {ji} , (3.24) where there is no anti-symmetric part σ (k) [ji] in the definition (3.14). By substituting ψ = mCiCj/2 into (3.10), we obtain 0 = 1 2 ( σ (k) li ∂lvj + σ (k) lj ∂lvi ) + 1 2 ( σ (c) li ∂lvj + σ (c) lj ∂lvi ) + χ ( mCiCj 2 ) , (3.25) 13401-7 K. Takechi, K. Yoshida, T. Arimitsu where σ(k) and σ(c) are given by (3.14) and (3.15), respectively, and χ(·) by (3.12). σ(c) and χ(mCiCj/2) reduce to functions of ∂jvi, ω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 ∂jvi and ωji are so small that the second or higher order terms in ∂jvi and ωji can be neglected. Then we arrive at the equations, σ (k) {ji} = 2µkEji , (3.26) σ (c) {ji} = 2d d+ 2 g0ν [ (4η1 + dη2)ρa ( Tt π )1/2 + (2η1 + dη2)µk ] Eji , (3.27) σ (c) [ji] = 2dg0νρa ( Tt π )1/2 η2Rji , (3.28) where µk = 2(πTt) 1/2ρa νg0 {2−d−1(d+ 2) + νg0 4 [6η21 − 4η1 + 2dη1η2 − (d− 2)η2 − 2η22 ]− η2 2Tr 2KTt νg0} {2(3 + d)η1 + (d− 1)(d+ 3)η2 − 6η21 − 4dη1η2 − (d2 − 3)η22 + η2 2Tr KTt } , (3.29) and ν is the volume fraction which may be related to the number density n of the particle as ν = πd/2 Γ ( d+2 2 )adn. (3.30) Now, the constitutive equation for the total stress tensor σ is given by σji(Eji, Rji,Ωkji = 0) = −ρTtδji + σ (k) {ji} + σ (c) {ji} + σ (c) [ji] (3.31) with (3.26)–(3.28). From symmetry, we have µkji(Eji, Rji,Ωkji = 0) = 0 1. The obtained constitu- tive equations are a restricted version of those in [10] in the sense that constants ∂jvi, ωji, Tt and Tr are assumed, and they are a generalized version in the sense that d is arbitrary whereas d = 3 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 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 (α) i = vi(r (α)) + C (α) i , w (α) ji = ωji(r (α)) +W (α) ji , (4.1) where C (α) i andW (α) 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, vi(r (2)) and ωji(r (2)) can be approximated by its Taylor series about x(1) up to the first-order terms, i.e., vi(r (2)) = vi + 2akj∂jvi , ωji(r (2)) = ωji + 2aklΩlji , (4.2) where we have omitted writing the argument x(1) in vi(x (1)), ωji(x (1)), ∂jvi(x (1)) and Ωkji(x (1)). From (3.4), (4.1) and (4.2), the velocity difference ξi between the two surfaces at the contacting 1Expand µkji(Eji, Rji,Ωkji = 0) in the power series of Eji and Rji. Note that the coefficient tensors are of odd degree. Since there is no special direction in the granular system, the coefficient tensor must be isotropic. Isotropic tensors of odd degree are 0. 13401-8 Constitutive equations for granular flow point is now written in terms of ∂jvi, ωji, Ωkji, ki, C (α) i and W (α) ji . Assume that the fluctuations C (α) i and W (α) ji can be neglected and that ki is isotropically distributed. Then, we have ξ := 〈〈ξiξi〉〉1/2 = ( 23 d )1/2 aω̂, (4.3) where ω̂ := [ d 2(d+ 2) EjiEji + 1 2 RjiRji + a2 2(d+ 2) (ΩjjiΩlli +ΩkjiΩkji +ΩkjiΩjki) ]1/2 , (4.4) and 〈〈· · · 〉〉 denote the average over ki. Here, we have used 〈〈ki〉〉 = 0, 〈〈kjki〉〉 = 1 d δji , 〈〈klkjki〉〉 = 0, (4.5) 〈〈klkmkjki〉〉 = 1 d(d+ 2) (δlmδji + δljδmi + δliδmj), (4.6) which can be derived using the general expressions of isotropic tensors and kiki = 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 Nc be a number of contacting points per a particle. Then, the energy dissipation per unit volume per unit time −q is given by − q = Nc 2 ρ m µfξ, (4.7) where the factor 1/2 is introduced to cancel the double counting of the contacting points. Note that the pressure pc due to the contacting force can be estimated by pc = νNcf S , (4.8) where ν is the volume fraction and S is the surface area of the particle. From (4.3), (4.7) and (4.8), one arrives at − q = (2d)1/2µpcω̂. (4.9) The estimate of pc is given as a homogeneous function of Eji, Rji and Ωkji, i.e., pc(aEji, aRji, aΩkji) = aζ ′ pc(Eji, Rji,Ω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 Φ = (2d)1/2µpcω̂, (4.10) which is a homogeneous function of Eji, Rji and Ωkji, i.e., Φ(aEji, aRji, aΩkji) = aζΦ(Eji, Rji,Ωkji) with ζ = ζ′ + 1. With the help of Euler’s homogeneous function theorem, Φ = 1 ζ ( ∂Φ ∂Eji Eji + ∂Φ ∂Rji Rji + ∂Φ ∂Ωkji Ωkji ) , (4.11) the choice σ{ji} = 1 ζ ∂Φ ∂Eji , σ[ji] = 1 ζ ∂Φ ∂Rji , 1 2 µkji = 1 ζ ∂Φ ∂Ωkji , (4.12) is made for constitutive equations which are consistent with (2.17). 13401-9 K. Takechi, K. Yoshida, T. Arimitsu 5. 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 g0 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 ∂jvi 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 (1) 0 in (3.20), different temperatures Tt and Tr 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 Ci andWji 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 = kiξi and ξt = √ ξiξi − ξ2n. Walton and Braun [24] derived β ≃ { −1 + ( 1 + 1 K ) µ(1 + e) cotϑ (ϑ > ϑ0) β0 (ϑ < ϑ0) , (5.1) 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 Φ(Eji, Rji,Ω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 ζ 6= 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. 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 d of the configuration space is 2. All the particles are disks of the same radius a. The mass 13401-10 Constitutive equations for granular flow is uniformly distributed inside the particles and thus the massm and the moment of inertia I of the particles are related as K = I/ma2 = 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 Lx and Ly. The area fraction ν is given by ν = Nπa2/LxLy. We fix Lx/a = Ly/a = 200 for all the simulations except for the case ν = 0.1 in which Lx/a = Ly/a = 600. The time scale associated with the collision is given by tcol = (m/κ)1/2. The time step ∆t of DEM is set to 10−2tcol, which is sufficiently smaller than tcol 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 tdamp = ∞. 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 coef- ficient. 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, ci(x = L, y, t) = ci(x = 0, y, t), (6.1) ci(x, y = L, t) = ci ( (x − Lyγ̇t) mod Lx, y = 0, t ) + Lyγ̇δix , (6.2) 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), vi(x, y) and ωji(x, y) are defined by the spatial averages, ρ(x, y) = mN∆ ∆x∆y , vi(x, y) = 1 N∆ ∑ α c (α) i , ωji(x, y) = 1 N∆ ∑ α w (α) ji , (6.3) 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 6 x(α) < x +∆x/2 and y−∆y/2 6 y(α) < y+∆y/2. The notations Eji(x, y), Rji(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., vi(y) = (γ̇y+A)δix where A is a constant. Here, we have chosen ∆x = Lx so that vi is independent of x. In terms of Eji(y), we have Eji(y) = γ̇ 2 (δjyδix + δiyδjx) . (6.4) It was found in the simulations that, when Eji(y) relaxed to a uniform field, ωji(y) also relaxed to a uniform field. Consequently, we have Rji(y) = R (δjyδix − δiyδjx) , Ωkji(y) = 0. (6.5) 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 ∂[jvi], 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 6 0, i.e., ωyx > ∂[yvx], at the quasi-steady states. 13401-11 K. Takechi, K. Yoshida, T. Arimitsu Figure 2. Stresses σ± normalized by mγ̇2 as functions of R/γ̇ for the mean shear γ̇tcol = 0.5× 10−3, 1.0× 10−3 and 2.0× 10−3. The friction coefficient µ and the area fraction ν are fixed to µ = 0.3 and ν = 0.7. Let us consider a line xi = h where h is a constant. In the simulation, the kinetic contribution to the stress tensor, σ (k) ji , is estimated by σ (k) ji (h) = − 1 Lj ∑ α (h) sgn(v (α) j ) mv (α) i ∆t , (6.6) where the summation ∑(h) α is taken over the particles α which cross the line xj = h during the time step ∆t. The contribution of contacts to the stress, σ (c) ji , is estimated by σ (c) ji (h) = 1 Lj ∑ α,β (h) Fα→β i , (6.7) where Fα→β i is the contacting force applied to the particle β by the particle α and the summation ∑(h) α,β is taken over the contacting pairs {α, β} such that x (β) j 6 h < x (α) j . The mean stress σji(h) over a line xj = h can be averaged over h to give the mean stress σji of the whole domain, σji = σ (k) ji + σ (c) ji , (6.8) with σ (k) ji = − 1 LxLy ∑ α mv (α) j v (α) i , (6.9) σ (c) ji = 1 LxLy ∑ α,β Fα→β i ( x (α) j − x (β) j ) , (6.10) where the summations ∑ α and ∑ α,β are, respectively, taken over the particles and the contacting pairs. Note that the pairs {α, β} and {β, α} are not distinguished and are counted only once in the summation ∑ α,β. Under the conditions (6.4) and (6.5), the constitutive equations read σ± = σ±(γ̇, R), (6.11) where we have introduced the notations σ+ := σ{yx} and σ− := σ[yx]. There are three time scales involved in the dynamics of the present system. They are tγ̇ := \|γ̇\|−1, tR := \|R\|−1 and the time scale of collision tcol = (m/k)1/2. In the simulations, we set tγ̇ ∼ tR ≫ tcol. Let us assume that, in 13401-12 Constitutive equations for granular flow Figure 3. (a) The value of the radial distribution function at r = 2a, g0 = g(2a), in the simulations for various ν and R = 0. Dashed line indicates the value of the low density expansion to the first order. (b) The normalized pressure p/mγ̇2 and the g0 as functions of \|δν\| = \|ν − νJ\|, where νJ = 0.84. The slope of the lines indicate the scalings in theory in [15]. such a case, σ+ and σ− do not depend on tcol. Then, a dimensional consideration yields similarity forms, σ±(γ̇, R) = sgn(γ̇)m\|γ̇\|2σ̃± ( R γ̇ ) = m\|γ̇\|γ̇σ̃± ( R γ̇ ) , (6.12) 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 γ̇tcol = 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 γ̇tcol = 1.0× 10−3. The value of the radial distribution function at r = 2a, g0 = g0(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 g0 at a thermal equilibrium state can be given up to the first order of ν as g0 = 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 g0 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 g0 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 g0 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 g0 ∼ \|δν\|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. 13401-13 K. Takechi, K. Yoshida, T. Arimitsu Figure 4. Probability function P (nint) of the interfering number nint, the maximum number of particles contacting with the colliding particle-pair simultaneously. For each particle-pair collision, we can define the interfering number nint, which is the maximum number of other particles contacting with the particle-pair simultaneously during the collisions of the particle-pair. nint = 0 implies that the particle-pair collision is not interfered by the other particles. In figure 4, the probability P (nint) 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 (nint) has a peak at nint = 0 for ν = 0.1 and 0.7 is consistent with the kinetic theory. However, the peak shifts to nint = 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]. Summarizing the above observations, we can consider that 0.7 6 ν 6 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 γ̇tcol. Although the figures for γ̇tcol = 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]. 6.2. 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 (δij − ǩiǩj)ξ̌j = −β(δij − kikj)ξj , (6.13) 13401-14 Constitutive equations for granular flow Figure 5. Averaged coefficient of roughness β̄(ϑ) as a function of the oblique collision angle ϑ for various ν and µ. The solid line indicates the theoretical function (5.1). The dashed line shows P (ϑ) = cos ϑ, the form for the random collisions without spin. The dot-dashed line shows P (ϑ) of the modeled collision with ς = 0.5. Figure 6. Averaged coefficient of roughness β̄ obtained from the simulations (circle symbols) and the coefficient of roughness βeff fitted to the simulation data on the basis of the kinetic theory by Lun (square symbols) for various ν and µ. The error bars indicate ∆β. where ǩ is the unit vector parallel to ř(2) − ř(1). Since a collision has a finite duration in the simulations, k and ǩ 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 6 β̄(ϑ) 6 1 which is consistent with the assumption −1 6 β 6 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 13401-15 K. Takechi, K. Yoshida, T. Arimitsu with the relative velocity c. Let b be the impact parameter and Θ(−π/2 6 Θ 6 π/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 6 b 6 a, we have for the probability density P̃ (ϑ) of ϑ (−π/2 6 ϑ 6 π/2) as P̃ (ϑ) = 1 2 1 (1 + tan2 Θ)( √ 1 + tan2 Θ− ς tanΘ) cosϑ , (6.15) with tanΘ = tanϑ+ ς √ tan2 ϑ+ (1− ς2) 1− ς2 . (6.16) In figure 5, the symmetrized probability density function P (ϑ) := P̃ (ϑ) + P̃ (−ϑ)(0 6 ϑ 6 π/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 ν 6 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, ap- proximate ϑ 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 σ (c) + , σ (k) + , and σ−(= σ (c) − ) normalized by mγ̇2 as functions of R/γ̇ obtained in the simulation are given for various ν(= 0.1, 0.7 and 0.8) and µ(= 0.3 and 0.8) in figures 7 and 8. The constitutive equations (3.26), (3.27) and (3.28) from the kinetic theory are also given in the figures. Here, the values obtained in the simulations are used for the “temperatures” Tt and Tr, and the value of the radial distribution function at r = 2a, g0 = g0(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 σ (k) + where the relative deviations are within 5%. The relative deviations for σ (c) + is about 20%. However, since \|σ(k) + \| is much larger than \|σ(c) + \|, the deviations for σ+(= σ (k) + + σ (c) + ) are almost determined by those of σ (c) + . 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 σ (c) + /σ (k) + increases. We have σ (c) + /σ (k) + ∼ 20 for ν = 0.7. In spite of this high ratio, the kinetic theory still gives fairly good estimates of σ (k) + and σ (c) + whose relative deviation from the simulation results are about 20% for σ (k) + and 10% for σ (c) − . As expected from the fact that the peak of P (nint) is shifted from nint = 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., C(β) = ∑ i { [ σth + (Ri;β)− σsim + (Ri) ]2 σsim + (Ri)2 + [ σth − (Ri;β)− σsim − (Ri) ]2 σsim − (Ri)2 } , (6.17) 13401-16 Constitutive equations for granular flow Figure 7. Stresses σ (c) + , σ (k) + [(a),(c)], and −σ− [(b),(d)] normalized by mγ̇2 as functions of −R/γ̇ for the area fractions ν = 0.1 and 0.7. The kinetic friction coefficient µ = 0.3 for all the figures. Shown in (a),(c) are σ (c) + and σ (k) + in the simulation (white circle and square), in the kinetic theory (3.26), (3.27) with β = β̄ (black circle and square, interpolated with line) and in the kinetic theory with β = βeff (vertical and diagonal cross). Shown in (b),(d) are σ− in the simulation (white circle), in the kinetic theory (3.28) with β = β̄ (black circle, interpolated with line) and in the kinetic theory with β = βeff (vertical cross). 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 Ri = 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 +Cmin. We determined the error ∆β of β to satisfy C̃(β ±∆β) = 2Cmin. 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. 6.3. 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 (nint) is not nint = 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. 13401-17 K. Takechi, K. Yoshida, T. Arimitsu Figure 8. Stresses σ (c) + , σ (k) + [(a),(c)], and −σ− [(b),(d)] normalized by mγ̇2 as functions of −R/γ̇ for the area fractions ν = 0.1 and 0.7. The kinetic fraction coefficient µ = 0.8 for all the figures. Shown in (a),(c) are σ (c) + and σ (k) + in the simulation (white circle and square), in the kinetic theory (3.26), (3.27) with β = β̄ (black circle and square, interpolated with line) and in the kinetic theory with β = βeff (vertical and diagonal cross). Shown in (b),(d) are σ− in the simulation (white circle), in the kinetic theory (3.28) with β = β̄ (black circle, interpolated with line) and in the kinetic theory with β = βeff (vertical cross). Figure 9. Dissipation function Φ normalized by mγ̇3 as a function of −R/γ̇ in the simulation for µ = 0.3 (white circle) and µ = 0.8 (white square). Dashed lines are fitted functions of the form (6.20). Also shown are the normalized dissipation function (4.10) in Kanatani’s theory with the value of pc from the simulation for µ = 0.3 (black circle) and µ = 0.8 (black square). As a guide, interpolations of the calculated data points are given in solid lines. 13401-18 Constitutive equations for granular flow Figure 10. Normalized stresses σ±/mγ̇2 as functions of −R/γ̇ in the simulations and those based on the decomposition by Kanatani (6.19) for µ = 0.3 [(a)] and µ = 0.8 [(b)]. ν = 0.8 for both figures. Circle and square symbols denote normalized σ+ and σ− in the simulation, respectively. Thick and thin dashed lines are the fitted lines (6.27), (6.28) to σ+ and σ−, respectively. Thick and thin solid lines indicate the decompositions (6.25), (6.26) for σ+ and σ−, respectively. Under the conditions (6.4) and (6.5), the dissipation function Φ of (2.17) reduces to Φ(γ̇, R) = σ+(γ̇, R)γ̇ + 2σ−(γ̇, R)R. (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 pc. 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 pc. 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: σ+ = 1 3 ∂Φ ∂γ̇ , σ− = 1 6 ∂Φ ∂R . (6.19) The normalized dissipation function Φ̃ := Φ/mγ̇3 in the simulation can be fitted by a polynomial of r := R/γ̇, Φ̃ = A+Br2, (6.20) with A ≃ 45, B ≃ 9.4, (µ = 0.3), (6.21) A ≃ 52, B ≃ 15, (µ = 0.8), (6.22) (see figure 9). Then, the choice (6.19) implies that σ̃+ = A+ +B+r 2, σ̃− = A−r, (6.23) with A+ = A, B+ = B 3 , A− = B 3 . (6.24) 13401-19 K. Takechi, K. Yoshida, T. Arimitsu From (6.21), (6.22) and (6.24), one obtains A+ ≃ 45, B+ ≃ 3.1, A− ≃ 3.1, (µ = 0.3), (6.25) A+ ≃ 52, B+ ≃ 5.0, A− ≃ 5.0, (µ = 0.8). (6.26) On the other hand, we obtain A+ ≃ 45, B+ ≃ 1.0, A− ≃ 4.2, (µ = 0.3), (6.27) A+ ≃ 52, B+ ≃ 2.3, A− ≃ 6.3, (µ = 0.8), (6.28) 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 fig- ure 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. 7. 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 nint, β̄(ϑ) 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. Acknowledgements We are grateful to the anonymous referee for his/her enlightening comments and notification of several important papers, which led to significant expansion and improvement of the paper. 13401-20 Constitutive equations for granular flow References 1. Condiff D.W., Dahler J.S., Phys. Fluids, 1964, 7, 842; doi:10.1063/1.1711295. 2. Eringen A.C., J. Math. Mech., 1966, 16, 1. 3. Savage S.B., Jeffrey D.J., J. Fluid Mech., 1981, 110, 255; doi:10.1017/S0022112081000736. 4. Jenkins J.T., Savage S.B., J. Fluid Mech., 1983, 130, 187; doi:10.1017/S0022112083001044. 5. Jenkins J.T., Richman M. W., Phys. Fluids, 1985, 28, 3485; doi:10.1063/1.865302. 6. Jenkins J.T., Zhang C., Phys. Fluids, 2002, 14, 1228; doi:10.1063/1.1449466. 7. Yoon D.K., Jenkins J.T., Phys. Fluids, 2005, 17, 083301; doi:10.1063/1.2000768. 8. Saitoh K., Hayakawa H., Phys. Rev. E, 2007, 75, 021302; doi:10.1103/PhysRevE.75.021302. 9. Lun C.K.K., Savage S.B., Trans. ASME E: J. Appl. Mech., 1987, 54, 47; doi:10.1115/1.3172993. 10. Lun C.K.K, J. Fluid Mech., 1991, 233, 539; doi:10.1017/S0022112091000599. 11. Goldshtein A., Shapiro M., J. Fluid. Mech., 1995, 282, 75; doi:10.1017/S0022112095000048. 12. Mitarai N., Hayakawa H., Nakanishi H., Phys. Rev. Lett., 2002, 88, 174301; doi:10.1103/PhysRevLett.88.174301. 13. Hatano T. , J. Phys. Soc. Jpn., 2008, 77, 123002; doi:10.1143/JPSJ.77.123002. 14. Otsuki M., Hayakawa H., Prog. Theor. Phys., 2009, 121, 647; doi:10.1143/PTP.121.647. 15. Otsuki M., Hayakawa H., Phys. Rev. E, 2009, 80, 011308; doi:10.1103/PhysRevE.80.011308. 16. Silbert L.E., Ertaş D., Grest G.S., Halsey T.C., Levine, D., Phys. Rev. E, 2002, 65, 031304; doi:10.1103/PhysRevE.65.031304. 17. Zhang H.P., Makse H.A., Phys. Rev. E, 2005, 72, 011301; doi:10.1103/PhysRevE.72.011301. 18. Shundyak K., Hecke M., SaarLoos W., Phys. Rev. E, 2007, 75, 010301(R); doi:10.1103/PhysRevE.75.010301. 19. Kanatani K, Int. J. Eng. Sci., 1979, 17, 419; doi:10.1016/0020-7225(79)90078-8. 20. Kano J., Shimosaka A., Hidaka J. J. Soc. Powder Technol. Jpn., 1996, 33, 95 (in Japanese). 21. Huthmann M., Zippelius A., Phys. Rev. E, 1997, 56, R6275; doi:10.1103/PhysRevE.56.R6275. 22. Goldhirsh I., Noskowicz S.H., Bar-Lev O., Phys. Rev. Lett., 2005, 95, 068002. 23. Kranz W.T., Billiantov N.V., Pöschel T., Zippelius A., Eur. Phys. J., 2009, 179, 91; doi:10.1140/epjst/e2010-01196-0. 24. Walton O.R., Braun R. L., J. Rheol., 1986, 30, 949; doi:10.1122/1.549893. 25. Onsager L., Phys. Rev., 1931, 37, 405; 38, 2265; doi:10.1103/PhysRev.38.2265. 26. Labous L., Rosato A.D., Dave R.N., Phys. Rev. E, 1997, 56, 5717; doi:10.1103/PhysRevE.56.5717. 27. Lees A.W., Edwards S.F., J. Phys. C: Solid State Phys., 1972 5, 1921. 28. Bagnold R.A., Proc. Roy. Soc. London A, 1954, 225, 49; doi:10.1098/rspa.1954.0186. 29. Hirschfelder J.O., Curtis C.R., Bird R.B., The molecular theory of gases and liquid. New York Wiley, 1954. 30. Edwards S.F., Physica A, 1998, 249, 226; doi:10.1016/S0378-4371(97)00469-X. 31. Otsuki M., Hayakawa H. Preprint arXiv:1006.3597, 2010. 13401-21 http://dx.doi.org/10.1063/1.1711295 http://dx.doi.org/10.1017/S0022112081000736 http://dx.doi.org/10.1017/S0022112083001044 http://dx.doi.org/10.1063/1.865302 http://dx.doi.org/10.1063/1.1449466 http://dx.doi.org/10.1063/1.2000768 http://dx.doi.org/10.1103/PhysRevE.75.021302 http://dx.doi.org/10.1115/1.3172993 http://dx.doi.org/10.1017/S0022112091000599 http://dx.doi.org/10.1017/S0022112095000048 http://dx.doi.org/10.1103/PhysRevLett.88.174301 http://dx.doi.org/10.1143/JPSJ.77.123002 http://dx.doi.org/10.1143/PTP.121.647 http://dx.doi.org/10.1103/PhysRevE.80.011308 http://dx.doi.org/10.1103/PhysRevE.65.031304 http://dx.doi.org/10.1103/PhysRevE.72.011301 http://dx.doi.org/10.1103/PhysRevE.75.010301 http://dx.doi.org/10.1016/0020-7225(79)90078-8 http://dx.doi.org/10.1103/PhysRevE.56.R6275 http://dx.doi.org/10.1140/epjst/e2010-01196-0 http://dx.doi.org/10.1122/1.549893 http://dx.doi.org/10.1103/PhysRev.38.2265 http://dx.doi.org/10.1103/PhysRevE.56.5717 http://dx.doi.org/10.1098/rspa.1954.0186 http://dx.doi.org/10.1016/S0378-4371(97)00469-X K. Takechi, K. Yoshida, T. Arimitsu Матер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шнього крутильного моменту. Результат чисельних симуляцiй проаналiзовано на основi моделi мiкрополярного плину. В симуляцiях реалiзується поле однорiдного середнього зсуву, яке не є пiдпорядковане полю вихоровостi. Оцiнки напружеь, зробленi на основi кiнетичної теорiї Люна [Lun, J. Fluid Mech. 233(1991) 539], добре узгоджуються з результатами симуляцiй в областi низьких часток ν = 0.1, але узгодження погiршується, якщо ця величина зростає. Проте, оцiнки, зробленi в кiнетичнiй теорiї можуть бути пiдiгнанi до результатiв симуляцiй аж до ν = 0.7 шляхом ренормалiзацiї коефiцiнта шорсткостi. Для вiдносно густого гранульованого потоку (ν = 0.8), результати симуляцiй також порiвнюються з теорiєю Канатанi [Kanatani, Int. J. Eng. Sci 17(1979) 419]. Знайдено, що дисипативна функцiя i її декомпозицiя в матерiальнi рiвняння в теорiї Канатанi не узгоджується з результатами симуляцiй. Ключовi слова: гранульований потiк, матерiальнi рiвняння, мiкрополярний плин, кiнетичне рiвняння 13401-22 Introduction Equations for micropolar fluid Kinetic theory for collisional granular flow Kanatani's theory Comments on the kinetic theory and Kanatani's theory Numerical simulations Setting Comparison with the kinetic theory Comparison with Kanatani's theory Discussion

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

Репозитарії

Схожі ресурси