A simulational study of granular boundary flows in two dimension

COMPUTATIONAL MATERIALS SCIENCE ELSEVIER

Computational Materials

Science4 (1995)125-132

A simulational study of granular boundary flows in two dimension Sang Rak Kim Physics Department, Kyonggi University, Suwon, 440-760, South Korea

Received28 November1994;accepted23 December1994

Abstract

Granular shear flows between flat rigid boundary walls in two dimension have been studied using particulate dynamics computer simulation. The walls, modelled as flat rigid planes with an inelasticity and a surface roughness, are sheared as an energy source. Density, translational/rotational velocity and temperature profiles across the walls are calculated. The granular particles are aligned parallel to the walls if the wall widths are within a few times larger than the particle diameter, irrespective of the granular properties of constituent particles and walls. Significant velocity slips at the walls are observed for the smooth particles. At large separations between the walls, a nonuniform distribution of granular temperatures across the walls and thus a phase separation is observed.

1. Introduction

A granular particle is a fine discrete solid particle. A granular material is a collection of a large number of granular particles. It shows many interesting features such as particle size segregation [l], self-organized criticality [2]. When the particles in the granular material are in contact with another, they behave like an elastic solid forming an angle of repose in sand piles since the frictional bonds between granular particles can support the applied loads. As the stress in the granular system is increased, the granular blocks move along their slip planes, and the system shows plastic behaviors as in soil mechanics [3]. This is called a quasi-static region. If the stress is further higher, each granular particle moves freely and independently of each other, and the particles collide instantly. This region of motion is called a rapid-flow or grain-inertia flow [4]. We 0927-0256/95/$09.50

0 1995

SSDlO927-0256(95)00021-6

will consider here this rapid-flow region. Industrial examples are flows of grains, ores and pharmaceutical powders through pipes and channels; geophysical examples are snow avalanches [51 and slides of rock debris [6]. A granular temperature may be defined as T= 1/2(V*> where V=V- (v), the randomvelocity of a granular particle relative to the mean velocity (Y) and (...> is an appropriate ensemble average. The shear work on the granular system as an energy source is done through the external driving force. The collisions between granular particles which have an inelasticity and a surface roughness dissipate away the granular temperature into real thermodynamic heat. The magnitude of the granular temperature is thus determined by a trade-off between the two mechanisms of heat generation and dissipation [4bl. It is well-known that the fluids in confined regions show remarkable behaviors [7-91 which

ElsevierScience B.V. All rights reserved

are quite different from that of homogenous bulk fluids. They show, for example, a nonuniform density distribution and adsorption effect. There are many common examples of fluids confined in the spaces of the size of a constituent particle in natural and technological products and processes, i.e. the molecular-level lubrication and friction. Experimental characterization of fluids in this regime is, however, quite difficult. Due to complex interactions and structures, large deviations from equilibrium and the lack of symmetry, computer simulation on model systems becomes an important tool to test ideas and to supplement real experiments in trying to understand the properties of fluids confined on the molecular scale. There are some theoretical efforts [10] to predict the fields of density, velocity, granular temperature and stress components in Couette flow. These theories do not include the rotational motion of granular particle, but they report good agreements with simulational and experimental results. Campbell and Brennen [ll] calculated, using computer simulation, the density velocity and granular temperature for the case of without rotation of granular particles. And Zhang and Campbell [12] reported the yield condition of the primary transition from solid to fluid behavior in two-dimensional granular flows. Recently, Campbell [13] studied in detail the ways that torques are applied to the particles at the flat boundary, assuming two boundary conditions, Type A and Type B. Two-dimensional flows of granular materials are convenient to study the behaviors of rapid granular flows under shear. In this work, we consider just a two-dimensional shear cell of granular particles confined by flat parallel walls. To make a shear, walls at top and at bottom are moved with constant speed in opposite directions. The walls are rigid and flat, but at collisions with a particle, there is an inelasticity and a surface roughness. In section 2, the computer simulational method is explained. The novel simulational results and the corresponding disscussions are presented in section 3 and finally the conclusions are followed in section 4.

2. Computer simulational

method

‘The granular system is composed of N identical granular particles, which are modelled as a circular disk with a diameter u, a mass of m, and a moment of inertia I. Each particle i has a translational velocity, V, and an angular velocity about its center of mass, 0,. Let us consider any two particles 1 and 2. They collide when the condition 1r, -rz( = CT meets. At collision, the total relative velocity at contact is defined g,, = “12 -$X0,

(I)

where Y,~ =v, - v2 and R=w, +w,, and i is a unit vector connecting from the center of particle 1 at r, to that of particle 2 at rz. The collisions between two granular particles are characterized by a coefficient of restitution E and a coefficient of surface roughness B which are defined as follows. Just after the collision, we have the relations &g;,=

-Ei.g,,,

(2)

f xg;, = -Bt% xg,,,

(3) where the primed ones denote the post-collisional values. Thus EC0 5 E I 1) and B( - 1 2 B 2 I) characterize the incomplete restitution of normal and tangential component of g,,, respectively. In other words, E models the linear velocity changes and B characterizes the rotational behavior during the collision. Note that the case of perfectly smooth particles corresponds to B = - 1. In that case the tangential component of g,, is unchanged by the collision. From Eqs. (l)-(3), we can write v; -v,

= -(vi-v*) =

-772v12

-

(7,

++Cfl,

4

-

Qq~.V,,)

(4)

-w,=w;-w2

mu2

-

?172@

(5)

S.R. Kim/Computational

Materials Science 4 (1995) 125-132

where 7, = (1 + E)/2, nZ = ((1 + B)/2XK/(l+ K)), and K = 41/mo2,a radius of gyration. Particle simulational methods [141 is used to simulate the granular systems. The intergranular force between granular particles is a hard sphere type and the governing equation of motion of a granular particles is a Newtonian. For the sake of simplicity and efficiency, we have considered a two-dimensional system in y-z plane. Further extension to three-dimensional system is quite straightforward. The walls are located at z = 0 and at z = h. The planar shear flow is driven by moving walls. In other words, to generate a shear flow, constant wall velocities Va are applied in the opposite directions. The upper wall is moving in +y direction and the lower in -y direction. Standard periodic boundary condition is applied along y-direction. The collisions with other particles are characterized by E, and B, and the collisions with walls are characterized by E, and B,, respectively. The energy of the granular system is dissipative through the inelastic collisions and thus the system is in a nonequilibrim steady state. Hence, there is no need of an artificial thermostat to maintain a constant temperature, contrary to the usual nonequilibrium molecular dynamics methods [14]. The numbers of particles we have taken are N = 200-2000. The relevant parameters of the

127

system and the simulational results are summarized in Table 1. All the simulational values are expressed in the usual reduced units. Thus nondimensional forms can be achieved using m, u, g as a basic unit of mass, length, and acceleration, respectively, where g is a gravitational acceleration constant. Hence, time t is, for example, expressed as in the unit of (a/g)1/2. g is introduced because it is a important constant in granular system. Campbell [13] studied two fictitious flat boundary types of a two-dimensional granular flow with the aim of trying to understand the possible effects of the boundary conditions on the bounded granular flow. The two boundary conditions, Type A and Type B, differ largely in the way that they apply torques to the particles during collisions with a wall. During a particle-wall collision, the Type A boundary applies the force at the particle surface, and thus the largest possible torque to the granular particle, while Type B boundary applies the force directly to the particle center, resulting in the application of zero torque. Type A boundary condition is rather artificial in that it corresponds to an infinite surface friction coefficient. He reported that through a small change of position of applying force, i.e. only a particle radius apart, there is, however, a large difference to the macroscopic behavior of the granular system on the whole.

Table 1 The simulational parameters and corresponding results. Here y = Vn / k is the applied shear rate, Tr is the translational granular temperature defined by Tt = 1/2(V,? + V,‘), T, is the rotational part of granular temperature defined by T, = 1/21(V,) where VW= w - (w > is a fluctuation component of angular velocity of a granular particle.

#l #2 #3 #4 #5 #6 #7 #8 #9 #lO #ll #12 #13

N

E”

B,

&

BW

K

Y

k

Tt

TW

200 200 200 200 200 200 200 200 200 200 200 200 2000

0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.8 0.9 0.9 0.9

0.9 0.9 0.9 0.0 - 0.9 0.9 0.9 0.9 0.0 -0.9 0.9 0.9 0.9

0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.8 0.9

0.9 0.5 0.0 0.9 0.9 0.9 0.9 0.5 0.9 0.9 0.9 0.9 0.9

0.5 0.5 0.5 0.5 0.5 0.5 0.5 1.0 1.0 1.0 0.5 0.5 0.5

5.0 5.0 5.0 5.0 5.0 1.0 50. 5.0 5.0 5.0 5.0 5.0 5.0

3.96 3.96 3.96 3.96 3.96 3.96 3.96 3.96 3.96 3.96 3.96 3.96 35.89

21.9 16.2 8.83 8.05 1.68 1.06 2827 21.9 11.0 3.07 16.4 22.9 670

33.0 17.9 9.10 9.90 2.19 1.29 3281 24.4 14.4 5.04 21.3 28.3 1064

12x

S.R. Kim /Computational

Matertals Science 4 (1995) 125-132

The system we studied is rather similar to the Type A boundary. We can, however, easily control an arbitrary surface friction through B,,, parameter and thus there is no artificial restriction on the wall boundary condition during a collision. Furthermore, we impose two heat source at bottom and at top and so we do not need, for simplicity, to apply normal stress to control the separation of walls. There is thus a symmetry about z = h/2 plane, which is evidently shown in the following Figs. 1-9. To get the informations of the microscopic behaviors of granular paricles, we compared the case of small wall separation with that of large one. And to see the microscopic profiles of local properties, we divide the system into very thin imaginary slabs parallel to the wall plane. The thickness of slab is taken to be 0.05~ to see the microscopic variation of properties. The particles with its center of mass located within a slab contribute to the properties. Thus we need a large time of sampling, at least, of the order of 1000 collisions per particle. The initial positions of the particles is located at a triangular lattice. The initial velocities are given at random to a specified temperature. But as time evolves, the local temperatures become quite different from the initial set-up value. We averaged the properties after a stationary condition has met.

3. Results and discussions 3.1. Small wall separations

When the separation between two walls is set at h = 3.96, the density profiles along the distance z from the bottom wall are shown in Fig. 1. Roughly speaking, the peak-to-peak distance is around 1.0. There are large peaks at the walls and two small peaks above the average density p = 0.84 at z = 1.6 and z = 2.4, respectively. The exitence of wall provides a layered density distribution corresponding to the successive neighbors from the wall. The gap between two small peaks in central region is narrower than that of between large end peak and nearby small peak due to

0

g

2' 1.5 1' 0.5 0, 0.5

1

1:5

2

2.5

3

3.5

z / sigma

Fig. I. Density profiles at h = 3.96. E, = 0.9, E, = 0.9.: (13) E,,=O.O. B,=O.9; (x) B,= -0.9, q&,=0.9; (a) B,=0.9. B, =o.s; (ml B,=O.9, B,=0.9;(*) B,=O.9, B,=O.O.

structural reconstruction of the colliding granular particles. This is consistent with the results of Davis et al. 1153for a simple fluid. However, the results are largely independent of the granular parameters E,, E,, B, and B,. We can see a broken triangular structure of the colliding particles in Fig. 2. The flat wall acts as a refuge of the flow field and thus the particles are gathered around the walls. This stable layers thus form the large peak at wall. Once the layer at wall is made, the small peaks can be formed about one diameter apart. This distribution of density is mainly determined by the separation h of walls, whether the system is in equilibrium or not. The small peak at z = 1.6 and z = 2.4 in Fig. 2 are due to a weak formation of second layer. The stream (longitudinal) flow properties [16171 of the granular system are shown in Fig. 3. There are three distinct regions in velocity profiles. The first region near center has a linear velocity profile which menifest a characteristic of shear flow, the second region has approximately zero slope and the third one very near at walls

Fig. 2. A snapshot h = 3.96.

of the configuration

of granular

system at

S.R. Kim/Computational

-51 0.5

1

1.5

2

2.5

3

Materials Science 4 (1995) 125-132

3.5

z / sigma

-20 0.5

1

129

1.5

2 z

2.5

3

3.5

/ sigma

Fig. 3. Translational velocity V, profiles at h = 3.96, E, = 0.9, E, = 0.9. The symbols are the same as in Fig. 1.

Fig. 4. Angular velocity VW profiles at h = 3.96, E, = 0.9 Also the same symbols as in Fig. 1.

E, = 0.9,

has again a linear velocity profile. The stream velocity profile is quite different from the linear velocity profile in Couette flow, and from the quadratic in Poiseuille flow of viscous fluids. We have also found that the granular particles seem to have some velocity slips at walls. Surprisingly, for the case #5 in Table 1, there is a large slip velocity with an amount of N 2.5. This is due to the smoothness of colliding particles which gives consequently rise to insufficient momentum transfer between particles. For some cases in Table 1, we have changed the radius of gyration K of the granular particles.

We have found that the larger K, the larger momentum transfer, so that the larger slope of velocity profile. Also the velocity profile is independent of applied shear rate y when y is not so high. The angular velocity profiles are shown in Fig. 4. They are also largely independent of E, and E,. We have, on average, zero angular velocity near center, two positive peaks at z = 1.4 and z = 2.6 and a large negative value at walls. The case #5 in Table 1 shows also quite a different profile due to insufficient transfer of angular velocity. Thus there is a very small chance

Fig. 5. Snapshot of the configuration of granular system at h = 35.89, t = 4.0 after the initial triangular lattice structure. Initially, the darker ones positioned in the lower half part and the brighter ones are in the upper positions. Note a rough surface at central region.

130

S.R. Kim / Compuiadonul

Materials Scrence 4 (1995) 125-132

to lose an angular velocity through a collision with another granular particle and the distribution of angular velocities is quite different. If we assume a local equilibrium in velocity components near wall, the rotational velocity at wall is given by wu = -(l + &)/Cl + K) (2V,/a). This is shown in Fig. 4 except #5 case. Each component of the granular temperature has no structure at all and is nearly uniform across the walls because the effects of driving walls transfer completely to the fluctuations of the velocitiy components. But the case #5 is still not in uniform state. For example, the distribution of T, is quite different and has some structure. This is also explained by the fact that there is a small chance of transfer of angular momentum due to the smoothness of colliding particles.

3.2. Large wall separation From now on, we consider the case with large separation h = 35.89. The initial positions of the granular particles at t = 0.0 are a triangular lattice. A snapshot of the positions of the granular particles at t = 4.0, corresponding to 159.4 collisions per particle, on average, after initialization is shown in Fig. 5. As a shearing on the system is executed, the granular particles are scattered by the walls and moved into the central region and crowded there together. The density profile across the wall distance z is shown in Fig. 6. The density distribution is divided into three separate rigions. The peak at near-wall region is the stable constructive layer with a diameter width. In the region next to the wall with 8-9 times of a diameter, there is no stable layer structure whatsoever due to the destructive interference and thus a large fluctuation of velocities can be developed. In the central region there is an obvious structure of microscopic layering. The particles in the central region have an average coordination number about 5.5 corresponding to the glassy or amorphous structure [12]. These are further confirmed by nonuniform distribution of granular temperatures T, as shown in Fig. 9. It can be argued [181 that there exists an equation of the state of granular system like in a simple fluid. The steady-state stream velocity profile is

(I

5

10

20

15

25

30

35

z/sigma

Fig. 6. Density profile at h = 35.89. The average density is 0.65. At central region, they are something like a glassy or amorphous solid structure.

shown in Fig. 7. It is evident that there are two different slopes of stream velocities. The central solid region is stagnant and has a nearly zero slope. The outer fluid region has roughly a linear velocity profile. The velocity slope is changed about z = 10 where a significant change of density profile also occurs. The collision of a particle with a wall seems to be the same situation in the rarefied gas. There is also a large slip velocity

z/sigma Fig. 7. Stream

velocity

V,, profile at h = 35.89.

S.R. Kim /Computational

- 25. This corresponds, in fact, to the rarefied gas case [19] with a large mean-free path. In the solid region the velocity fluctuations or granular temperatures persist, but the density is too high to flow easily. Thus it forms a quasi triangular lattice structure. Therefore it is possible to see the structural phase change of granular flow from the stagnant amorphous solid to the fluid region. The phase boundary can be determined by a slope change of flow velocity as shown in Fig. 7. There are very small fluctuation of flow velocity, and thus shows a smooth behavior. Thus, roughly speaking, the thickness of stagnant solid is about 16. Also, as shown in Fig. 6, this stagnant region has a triangular lattice structure and thus a significant microscopic layered structure in every distance of a diameter of a particle. An extremely important point to note is some correlation between velocity and density of the bound flows. Thus it may be necessary to analyze the data with local density approximation approach of simple fluid 1201. The steady-state angular velocity profile is shown in Fig. 8. At central region there is effectively no rotational motion, but near the walls the granular particles have an opposite direction of rotation in a distance of particle diameter, which is the same behavior as in the small wall separation. The magnitude of angular velocity at wall is 20

10

0

-10 2 -20

-30

-40

-50 0

I

I

I

I

5

10

15

20

z I

I 25

I 30

131

Materials Science 4 (1995) 125-132

35

sigma

Fig. 8. Angular velocity profile VW at h = 35.89. At center there are nearly no rotation at all.

0

5

10

15

20

25

30

35

z/sigma Fig. 9. Translational granular temperature T,,, profile at h = 35.89. There are granular heat conduction from the end walls to the central region.

only about 45, less than the predicted. This is due to the fact that the particles are not at local equilibrium in their velocities near walls. A colliding particle gets a linear momentum and a rotational momentum rather coherently from the moving walls. The particle experience next collision with another particle with a long mean-free path larger than the diameter of granular particle. Thus the distribution of rotational velocities, shown in Fig. 8, shows a large negative values near the walls, and has a positive peak a few times of a diameter apart from the wall boundary. The next negative peak, but weak, occurs about a diameter apart. The central stagnant region shows nearly zero rotational velocities and small fluctuations. The steady state granular temperature profile is shown in Fig. 9. High granular temperature is observed near the walls and low temperature at central region. Due to the lack of symmetry and the dissipation of the granular system, there is an incomplete equipartition between each degree of freedom, i.e. not the same contributions to the granular temperatures [21]. Thus we have observed, in general, T,, > T,, where T,,,, = 1/2(V,*>, T,, = 1/2(VZ2). It is evident that there should be a granular heat conduc-

tion from both the moving walls to the central region. The moving wall acts as a granular heat source and the heat is dissipated mainly through the collisions between granular particles. The graph of 7’,,) shows the change of slopes. The central stagnant region have nearly zero slope and small fluctuations of granular temperature. The outer fluid region has large slope and strong fluctuations of temperature. The near-wall region have apparently zero slope but much larger fluctuations. These behaviors are the same for T.._. and T,.

4. Conclusions We have considered the granular shear flow dynamics between parallel flat rigid walls in two dimension. There are significant inhomogeneous structures, such as layerings at small wall separation h which is independent of granular properties of constituent particles and walls. Also the nonlinear velocity profiles are observed and thus there are slip velocities of granular particles at boundary walls dependent on the B, parameter. When B, = - 0.9, there is quite a small chance of losing an angular velocity through a collision and thus a large slip velocity. At large wall separation h, there is a large non-uniform temperture distribution and thus a granular heat conduction from the moving walls to the central region. At central region with lower temperatures there seems to be glassy or amorphous solid structure. By comparing between the results with narrow and wide separations of wall boundary, it is observed that there is some phase interface at large separations, a boundary between a fluid state and an amorphous solid state. There are still so many remaining problems; a simple extension to 3D system, an application of gravitation, different size/mass ratio of granular particles, more elaborate wall models and so on.

Acknowledgements

This work has been supported in part by NON-DIRECTED RESEARCH FUND, Korea

Research Foundation and the Basic Science Research Institute Program, Ministry of Education, 1993, Project No. BSRI-94-2430 and by the KOSEF through the SRC program of SNU-CTP.

References J.C. Williams, Powder Tech. 15 (lY76) 245; P. Meakin and R. Jullien, Physica A 180 (1992) I. P. Bak. C. Tang and K. Wiesenfeld, Phys. Rev. Lett. 59 (1987) 381; S.R. Nagal, Rev. Mod. Phys. 64 (1992) 321. A. Mroz, Intl. J. Numer. and Anal. Meth. in Geomech. 4

I I’)XO)4s. (a) S.B. Savage, Adv. Appl. Mech. 24 (1984) 289; (h) <‘.S. Campbell, Ann. Rev. Fluid Mech. 22 (1990) 57. E.J. Hopfinger. Ann. Rev. Fluid Mech. 15 (1983) 47. T.H. Erismann, Rock Mech. 12 (1979) 15. S. Murad, P. Ravi and J.G. Powles, J. Chem. Phys. 98 (190.3) 9771. B.N.J. Persson, Phys. Rev. Lett. 71 (1993) 1212. J.M. Georges, S. Millet, J.L. Loubet and A. Tonck, J. (‘hem. Phys. 89 (1993) 7345. C.K.K. Lun, S.B. Savage, D.J. Jeffrey and N. Chepurniy, J. Fluid. Mech. 140 (1984) 223; I(. Hui, P.K. Haff, J.E. Ungar and R. Jackson, J. Fluid. Mech. 145 (1984) 223; J.T. Jenkins and E. Askari, J. Fluid. Mech. 223 (1991) 497. C.S. Campbell and C.E. Brennen. J. Fluid. Mech. 151 (1985) 167. 1121Y. Zhang and C.S. Campbell, J. Fluid. Mech. 237 (1992) Ul. 1131 C‘.S. Campbell, J. Fluid. Mech. 247 (1993) 111. 1141 M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids (Clarendon. Oxford, 1987); J.M. Haile, Molecular Dynamics, Simulation: elementary methods, (Wiley, New York, 1992). 1151 H.T. Davis, I. Bitsanis, T.K. Vanderlick and M.V. Tirrell. Suercomputer research in chemistry and chemical engineering. ed K.F. Jensen and D.G. Truhlar (Washington D.C‘., Am Chem Sot. 1987). Ih] S. Hess and W. Loose, Physica A 162 (1989) 138. 171 M. Sun and C. Ebner, Phys. Rev. A 46 (1992) 4813. 181 C.S. Campbell and C.E. Brennen, Tran. ASME: J. Appl. Mech. 52 (1985) 172. 191 J. Koplik, J.R. Banavar and J.F. Willemsen, Microscopic Simulation of Complex Flows, ed M. Mareschal (Plenum, New York, 1990). [ZO] J.P. Hansen and I.R. McDonald, Theory of Simple Fluids. 2nd ed. (Academic Press, London, 1986). 1711 S.R. Kim and L.V. Woodcock, J. Stat. Phys. 71 (1993) 143: S.R. Kim, J. Korean Phys. Sac. 26 (1993) 268.