- Email: [email protected]
Local fluctuations in a fluidized granular medium
PHYSICA ELSEVIER
Physica A 240 (1997) 361-373
Local fluctuations in a fluidized granular medium K. H e l a l , T. B i b e n * ,
J.P. H a n s e n
Laboratoire de Physique Ecole Normale SupOrieure de Lyon, 69364 Lyon cedex 07, France 1
Abstract
Molecular dynamics simulations have been carried out on a two-dimensional model granular material made up of inelastically colliding discs subjected to gravity. The granular system is fluidized by a constant energy input from a vibrating base, and a stationary non-equilibrium state is achieved. Density and granular temperature profiles are determined; while the former goes through a maximum at intermediate altitudes, the latter exhibits a minimum at a higher altitude, beyond which the granular temperature is found to increase. Local density fluctuations are characterized by the time-displaced density autocorrelation function. Longitudinal compression (sound) waves propagating in the horizontal direction are found to be strongly damped, while they are overdamped in the vertical direction at the longest accessible wavelength. P A C S : 46.10.--z; 05.60.+w; 05.40.÷j
1. Introduction
A granular material, like a sandpile, may be set in motion under the action of gravity, as in a hopper, or of externally imposed vibrations. Under appropriate conditions, the applied vibrations may induce surface or bulk fluidization of the granular medium, which then exhibits a behavior reminiscent of fluid flow. At the "microscopic" level, the grains undergo inelastic collisions, and a stationary non-equilibrium state is achieved when the dissipated energy is exacty compensated by the continuous energy input from a vibrating base. A number o f recent experiments and numerical simulations carried out on two-dimensional model granular media have examined global and collective properties o f such flows, like the granular temperature, the mean height of the centerof-mass, the propagation of compression waves, or the energy flux through the medium [ 1-6]. Use o f a high-speed camera allows the motion of individual grains confined between vertical glass plates to be monitored as a function o f time [7,8], and very * Corresponding author. E-mail: [email protected]. I Unit6 de Recherche Associre 1325 du CNRS. 0378-4371/97/$17.00 Copyright @ 1997 Elsevier Science B.V. All rights reserved PI[ S0378-4371 (97)00 159-3
362
K. Helal et al./Physica A 240 (1997) 361-373
recently this technique has been used for a quantitative characterization of local density fluctuations in a two-dimensional, vibrated array of steel beads via time-displaced correlation functions familiar from the study of molecular motions in fluids in thermodynamic equilibrium [9]. Some striking analogies between the relaxation of density fluctuations in the fluidized granular medium and in an atomic fluid, but also some significant differences, were discovered by the experiments reported in Ref. [9]. In this paper we present results of molecular dynamics (MD) simulations [13] of a system of inelastically colliding discs subjected to external vibrations and gravity; these simulations were carried out in close connection with the experiments of Ref. [9]. The main objective is to gain a more complete understanding of density and energy fluctuations in a simple model which captures the essential features of the experimental system, rather than to seek a faithful reproduction of the experimental data. In particular, the model under consideration is free of the complications which may arise from the friction experienced by the beads on the confining glass plates. The MD simulations also allow an easier control of the sample-size-dependence of the results as well as better statistics.
2. Model and simulations
We consider a system of N impenetrable discs of diameter tr and mass m confined to the x - y plane. The centers of the discs are confined to a periodically repeated slab of width L in the horizontal (x-direction), and to altitudes y > 0 in the vertical direction (see Fig. 1). The discs are subjected to the constant vertical gravitational force - m y , and suffer instantaneous collisions when the centres of a pair are exactly separated by a distance cr. The pre-collisional and post-collisional (primed) velocities of discs labeled
\_~
\
/
"i
"
s'
= I
..:
/ -
,÷-1
' _.
'
-r--7--T-7"- ~-'-'~7~--'
[
T 1"
.... ,
I
--r--
:+ ,]
k
u
Fig. 1. S c h e m a t i c r e p r e s e n t a t i o n o f the r e c t a n g u l a r o b s e r v a t i o n w i n d o w s L × ( L / 4 ) a n d ( L / 2 ) x ( L / 4 ) u s e d in the simulations. T h e discs are not s h o w n to scale; in practice L = 40tr.
K, Helal et al. I Physica A 240 (1997) 361~73
363
1 and 2 are related by -+-r)[tJl2" (Vl - - I J 2 ) ] R I 2 , • (vl - 1~2)]/112 ,
U'1 =
U 1 -- ½ ( l
U~ =
U2 q- 1 ( 1 -}- r ) [ h l 2
(1)
where hi2 is the unit vector along the direction joining the centres of the colliding discs and r is the coefficient of restitution. Elastic collisions are recovered when r = 1, while in model granular systems, r is typically of the order of 0.95 [7]. The collision rules (1) conserve momentum, but kinetic energy is dissipated when r < 1. This dissipation leads to a rapid "cooling" and "clumping" of grains with given initial velocities [10], but may be compensated by injecting kinetic energy at the y = 0 base through an external vibrator, resulting in a stationary non-equilibrium state. The vertical vibrations imposed in real experiments are sinusoidal in time, but the resulting disc-base collision kinematics are not easily reconciled with the "event-driven" kinematics of binary collisions between discs. For the sake of simplicity, we have modeled the base vibrations and resulting energy input via either of the following two schemes: (a) The base is fixed at y = 0 , but the discs undergo inelastic collisions with the base characterized by a coefficient of restitution R > 1; under these conditions the post-collisional velocity components are vx' = vx; ty" = - R v v , resulting in a relative gain in kinetic energy Ae¢/ec = R 2 - 1 > 0 of a colliding disc. (b) The sinusoidal variation of the base altitude is replaced by an asymmetric sawtooth characterized by its amplitude A and its period to: y ( t ) = A [ t / t o ] - A/2, where Ix] denotes the fractional part of x. Under these conditions a disc gains kinetic energy of at least 2mA2to 2 upon each collision with the base. The simulated sample involved N = 400 discs, corresponding to 10 horizontal layers of 40 discs, comparable to the experimental situation explored by Warr et al. [7-9]. Most simulations were carried out for r = 0.97. A well fluidized stationary state was then obtained with R = 1.5 or ~ = A / ( g t 2 )~_ 5-10. Results obtained with the two-energy input schemes were very similar for comparable rates of base-particle energy transfer; most results reported below were obtained from simulations using energy input scheme (a). Starting from some reasonable initial condition, a typical MD run involved some 10 7 collisions between discs, and statistical averages were taken over the resulting phase space trajectory, beyond an initial "equilibration" period.
3. Density and granular temperature profiles The granular temperature T is defined by the mean-square velocity of the discs according to m
2
r = 5 (v )
(2)
and it is convenient to use the quantity r ' - - kBT - - 1 (1)2)
m
2
(3)
K. Helal et al./Physica A 240 (1997) 361-373
364 0.50
e - - - - o r=O.95;R=l.7 m . - - ~ r=O.97;R=1.5
,o,,o_o. e, 0.40 tO
/,
~,
~f
-'-
-'- r=O.97;vibrator
"\
]
I
0.30 :3 ~
.~_
-'~ 0.20 ¢.~ 0.10
0.00 0.00
10.00
20.00
30.00
40.00
y/o Fig. 2. P a c k i n g fraction profiles r/(y) as a function o f altitude y/a. Circles a n d dash-dotted curve: r = 0.95 and R = 1.7; squares a n d d a s h e d curve: r = 0.97 a n d R = 1.5; triangles a n d full curve: r = 0.97 and saw-tooth vibrator, with amplitude A = 2 / 3 6 a n d period to = 0.02s.
expressing the granular temperature in units of velocity squared. A properly reduced temperature, T* = -TI -,
9a
(4)
characterizes the ratio of the mean kinetic energy of the discs to their gravitational potential energy. The local density, or density profile, p(y) is defined such that p(y)L dy is the mean number of discs in a horizontal strip of area L dy, at the altitude y. It is convenient to define the reduced density profile p*(y)= p ( y ) a 2, and the local packing fraction, q ( y ) = ¼np*(y). Since the granular fluid is "heated" from below, through the energy input from the base, the granular temperature varies with altitude y. We have determined a "temperature profile" T*(y), by averaging the kinetic energy of discs contained in horizontal strips of thickness Ay=2a. Finally, we have computed the rate of energy dissipated per unit mass, w(y), via the loss of kinetic energy at each inelastic collision between discs, as a function of the altitude y. Typical examples of the three profiles r/(y), T*(y) and w(y) are shown in Figs. 2-4, and for two pairs of restitution coefficients (r = 0.97; R = 1.5) and (r = 0.95; R = 1.7). In both cases, the density profile goes through a maximum at an altitude Ym > 10a, in qualitative agreement with the experimental findings of Warr et al. [8]. In the more inelastic case (r =0.95), fluidization is less pronounced, i.e. the discs tend to concentrate at lower altitudes, where the local packing fractions are significantly enhanced relative to the less dissipative system (r = 0.97). More surprisingly, the temperature profile T*(y) is found to go through a minimum at an altitude well above Ym, and to increase at the higher altitudes, where the local density is very low. A faint indication of this unexpected behavior is found in one of the granular temperature profiles measured by Wan" et al.
K. Helal et al. IPhysica A 240 (1997) 361 373
365
25.0 o---~ r=0.95;R=1.7 "~----'~-~rr=00:97;;Rb~a?o r
20.0 \ ~ ~ , , ~
o
15.0
E o
10.0
to'} 5.0 ~"O,,
/O
o---o..e.e..e_..o_o_e_o_.e.-e-e--e-'~
0.0 0.0
L
10.0
20.0 y/~
L
30.0
40.0
Fig. 3. Granular temperature profiles T*(y) versus altitude y/a. The symbols have the same meanings as in Fig. 2.
1.0 ~--. 0.9 e r ~ 0.8 ~ ~\'~,\
.
.
. I" ,~ r=0.97;vibratorI Io---4) r=0.95;a=1.7 I
0.6
0.4
0.1 t 0.0 0.0
, 10.0
"13,, . ~ . , , ~ ~ 20.0 y/o
30.0
~ A 40.0
Fig. 4. Rate of energy dissipation, w(y), normalized by its maximum value Wmax, as a function of altitude. The symbols have the same meanings as in Fig. 2; the results for r = 0.97; R = 1.5 are not shown for visual clarity.
K. Helal et al./Physica A 240 (1997) 361-373
366
19.5 ~
.1¢
1-
14.5
0.2
\
E -~\ e'-
~
ca" 0.0
0.0
9.5
4.5
0.0
~
10.0
20.0
20.0
y/c
4( .0
30.0
40.0
Fig. 5. Comparison between the theoretical temperature and packing fraction profiles, based on the phenomenologicalcalculationof the Appendix (full curves) and the MD data (circles) for r = 0.97; R = 1.5. [8]. A phenomenological calculation based on hydrodynamic equations for dissipative granular media [11], reproduces the observed behavior of the density and temperature profiles very well, as shown in Fig. 5. Details of this calculation are given in the Appendix. The energy dissipation profiles shown in Fig. 4 decay to zero after a slight initial rise involving the first few particle layers above the base. Once more the decay is significantly faster in the case of the smaller coefficient of restitution. The concept of a local granular tempera~re, playing a role similar to the equilibrium temperature in atomic or molecular systems, is meaningful only if the velocity distribution is strictly Gaussian (or Maxwellian), i.e. ~b(v) °c exp ( - 2~7) exp
-2T']
(5)
To test this assumption, we performed statistical averages of the distributions of the horizontal (Vx) and vertical (Vy) components of the particle velocities over strips of thickness A y = 2 a at several altitudes. Typical examples are shown in Fig. 6. They show a significant difference between the distributions of the horizontal and vertical velocity components. While the former are symmetric around vx = 0 and Gaussian within statistical error at all altitudes, the distributions of the vertical velocities are seen to be close to Gaussian only at intermediate altitudes, while they are significantly skewed to the left (negative Vy) near the base, and to the right (positive Vy) at the highest altitudes. The width of the distributions, proportional to x/T-7, first shrinks with altitude, but broadens again beyond y-~ 20a, in agreement with the behavior of the
K. Helal et al./Physica A 240 (1997) 36l 373
50.0
•
•
i
367 !
(A) 50.0
25.0
0.0 -5.(
0.0
5.0
0.0 -5.0
0.0
5.0
i
(c)
(D)
50.0
Vv
5.0
0.0 -5.0
0.0
5.0
0.0 -5.0
0.0
i.O
Fig. 6. Velocity distribution functions (O(Vx) (crosses) and ~b(Vv) (triangles) at altitudes y/a= 1 (A), I1 (B), 23 (C) and 39 (D), from an MD simulation with r = 0 . 9 7 and R = 1.5. Velocities are in units where L = 1 and 9 : 10. The full curves are the best Gaussian fits to 4~(v~).
temperature profile shown in Fig. 3. Inspection of Fig. 6 leads us to conclude that the concept of a local granular temperature is reasonably well founded.
4. Local density fluctuations Microscopic density fluctuations in an atomic fluid in thermal equilibrium are routinely characterized by the time-displaced autocorrelation function of the Fourier components of the microscopic density, conventionally referred to as the "intermediate scattering function" in connection with inelastic neutron scattering experiments [t2]. For a closed or periodic system containing N atoms, the density autocorrelation function is defined as F ( k , t) = 1 ( p k ( t ) p - k ( O ) ) ,
(6)
where the brackets denote a thermal average and the Fourier components of the microscopic density are N
p k ( t ) = ~-~e ik~'(') , j=|
(7)
K. Helalet al./Physica A 240 (1997) 361-373
368
if rj(t) denotes the position of the jth atom at time t. Experimental or numerical measurements of F(k, t) provide a powerful tool to investigate collective particle motion on a scale 2 = 2n/k. In particular, propagating longitudinal sound (or compression) waves show up as damped oscillations in F(k, t) for sufficiently long wavelengths (2>>a). At shorter wavelengths (2 _~ a), density fluctuations are found to decay monotonously to zero. In a very recent experiment [9], F(k,t) has been measured in a 2d vibrated granular medium with the help of a high-speed camera, which allows the precise determination of individual grain positions within a rectangular window for a large number of successive instantaneous configurations. We have computed F(k, t) in our MD simulations, for several wave vectors k, and under physical conditions similar to those in Ref. [9]. In the simulations we used windows of area Lx x Ly. Lx was taken to be the total horizontal width of the system, with periodic boundary conditions at ±L/2, or one-half that width, with open boundaries at ±L/4. In the vertical direction, the window extended from y = 7.5a to y = 17.5a above the base, i.e. Ly 10or. Referring back to the density profile in Fig. 2, for ( r = 0 . 9 7 ; R = 1.5), it may be seen that this choice was made to study a sample as uniform as possible in the vertical direction. The accessible wave vectors k are of the form =
k
z
(nx ny ~ 27r\Lx L y J ' - - ,
- -
(8)
where nx, ny are integers. There are some fundamental differences between the density autocorrelation functions in a thermal fluid and in a fluidized granular medium: (a) The statistical average is to be taken over initial conditions which belong to an equilibrium ensemble for fluids, and to an ensemble of stationary non-equilibrium states in the case of fluidized granular matter; in both cases stationarity implies time translation invariance of the correlation functions. (b) In an isotropic fluid, rotational invariance implies that F(k,t) depends only on k = Ikl while translational invariance implies that only Fourier components of opposite wave vectors have non-vanishing correlations. Both invariances are broken in granular matter due to the preferred (vertical) direction, and the vertical inhomogeneity imposed by the base vibrations and gravity. We have nevertheless restricted our investigations to ( k , - k ) pairs and considered density fluctuations with horizontal and vertical wave vectors separately. (c) For equilibrium fluids, one normally considers closed or periodic systems, such that the number of atoms contributing to the density fluctuations (7) is fixed. In the present simulations, and in the experiments of Ref. [9], the granular sample within the window under observation forms an open system. This means that N, and the identity of the discs or beads within the window, fluctuate in time. Consequently, N is replaced by its mean value (N) in definition (6) (other conventions could be adopted). Moreover, the lack of space inversion (k ~ - k ) symmetry, combined with the open nature of the system, imply that F(k,t) is in general no longer a real function, except for t = 0, or for horizontal k-vectors, since the x ~ - x inversion symmetry is preserved. In the t = 0 limit, F(k,t) reduces to the static
K. Helal et al./ Physica A 240 (1997) 361-373
369
2.0
"5
o
,=
1.0
09 0
[
¢/)
~
0
0.0
Baus & Colot 0 Full window x Half window
I
I
i
I
I
2.0
4.0
6.0
8.0
10.0
12.0
q=kc Fig. 7. Static structure factor S(q) versus reduced wave number q = ka determined from particle coordinates in the observation windows shown in Fig. 1 and described in the text; the MD simulations are for r = 0.97; R = 1.5 and the mean packing fraction in the window is ~7= 0.462. Circles and crosses correspond to the full and half-windows, respectively; data points are for wave vectors of variable direction, stressing the isotropy of the structure factor; the full curve is the structure factor for a homogeneous fluid of hard discs at t / = q [14].
structure factor 1
S(k) = F(k, t = O) = 7-~ (PkP-k)
(9)
while in the long-time limit 1
2
t l i m F ( k , t ) = -~-~[(p~)] .
(10)
For horizontal wave vectors, this limit vanishes, since the mean density is constant in the x-direction, while for vertical wave vectors, (Pk) is the Fourier transform of the vertical density profile: jfl, "¢2
(Pk) =
eikYp(y) dy,
( 11 )
where Yl and Y2 are the lower and upper altitudes of the simulation window. Since p(y) is not uniform, (Pk) is non-zero, and the resulting limit (10) of F(ky, t) is non-zero. A few examples of results from the MD simulations are shown in Figs. 7-9. In Fig. 7 the computed static structure factor S(k) is compared to that of a homogeneous
K. Helal et aL /Physica A 240 (1997) 361-373
370
1.0
~" 0.05 (/) v
09
0.5
0.00
v'
,
i
I
0.0
i
ii
~
-~''-~- -
25.0
ii
50.
(O.tE
0'00.0
20,0
40,0
60.0
80.0
100.0
t/t E
Fig. 8. Normalized density autocorrelation function F(k,t) versus time (in units of the Enskog collision time te), for k = (27r/L, 0), and the full observation window in Fig. 1 (MD run with r = 0.97; R = 1.5). The corresponding spectrum S(k, 02) is plotted versus ~o • tE in the insert.
fluid of elastic hard discs in thermal equilibrium, as parametrized by Baus and Colot [14]. The agreement is seen to be excellent. Significant discrepancies were found in a similar comparison between the elastic hard disc results and experimental data [9]; these discrepancies are probably due to the fact that the granular medium in the experimental window was much more inhomogeneous than the corresponding simulated system. Both simulation and experimental data are surprisingly isotropic, i.e. for given [k], S(k) is practically independent of the orientation of the wave vector. Fig. 8 shows the real part of F(k, t) for the smallest horizontal wave vector, k---(2rt/L, 0). A small amplitude oscillation is seen to be superimposed on the slow decay of the correlation function, characteristic of a strongly damped longitudinal collective mode, which we tentatively identify with a compression (or "sound") mode propagating through the model granular medium. The spectrum of the correlation function obtained by Fourier transformation with respect to time, exhibits a very broad resonance at an angular frequency o2"~ 27.5 rads - l (in unit where L = 1) corresponding to a sound velocity c = og/k _~ 1.9 m s - 1 for particles of diameter tr = 5 mm roughly compatible with that predicted from the adiabatic compressibility of a homogeneous fluid of elastic hard discs with density equal to the mean density of the inelastic system considered here. A MD simulation carried out for elastic discs shows that the oscillations in F(k,t) are much more pronounced in that case, leading to a much sharper resonance in the density fluctuation spectrum. Fig. 9 compares the real parts of F(k,t) computed for horizontal and vertical wave vectors of the same modulus. As foreseen in our earlier discussion, the correlation function for the horizontal k decays to zero, while that for the vertical k goes to a non-zero limit, which agrees within statistical errors with the value computed from
371
K. Helal et al./Physica A 240 (1997) 361-373 1.0
1.0 0.9
v
iJ_
0,9
(A)
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0.0
0.0 -0.1
(B)
w
0.0
10.0 t/t E
20.0
-0.1
i
0.0
10.0
20.0
t/t E
Fig. 9. Normalized density autocorrelation function F ( k , t ) versus t/tE for: (A) k=(O, 87z/L); (B) k = (8~z/L,O), under conditions identical to those in Fig. 8.
Eq. (10). Data obtained for several wave vectors show that the rate of relaxation initially increases with Ikl, as in the case of fluids in thermal equilibrium. For all correlation functions which we have investigated, the imaginary part is very small at all times, remaining below the statistical noise level.
5. Conclusion The simulations described in this paper point to some strong similarities between the collective behavior of fluidized granular matter on a local scale (i.e. on the scale of a few grain diameters) and the corresponding properties of equilibrium fluids on a microscopic scale. A similar conclusion had been drawn on the basis of experiments reported in Ref. [9], and semi-quantitative agreement is found between the experimental and simulation data for comparable physical conditions. Both sets of results are very sensitive to the amplitude and frequency of the external vibrations; it is a matter of delicate balance between system size, rate of energy input and rate of bulk-energy dissipation to achieve a stationary fluidized granular state. In particular, the present work shows that the static and dynamical properties of a system of inelastic discs are very sensitive to the coefficient of restitution r, which has, in particular, a strong incidence on the density profile. Preliminary simulations for r = 0.92 point to the formation of a slowly varying inhomogeneity (long-lived "clumps") in the horizontal direction. The key findings of the present work may be summarized as follows: (a) The granular temperature profile goes through a minimum at intermediate altitudes, and increases in the low-density region at higher altitudes. This unexpected behavior is well reproduced by a hydrodynamic calculation (Appendix).
372
K. Helal et al. / Physica .4 240 (1997) 361-373
(b) The static structure factor measured over a range of altitudes where the system is nearly homogeneous (i.e. around the maximum in the density profile) is in perfect agreement with the known structure factor of a homogeneous fluid of hard discs in thermodynamic equilibrium. This indicates that the local structure in a stationnary non-equilibrium state of fluidized granular matter strongly resembles that of fluids in thermal equilibrium. (c) Longitudinal compression (or "sound") waves are strongly damped even for the longest accessible wavelength, compared to the corresponding behavior in an equivalent system of elastic discs; a small degree of inelasticity (r = 0.97 in the present calculations), coupled with the vertical inhomogeneity, is thus found to provide a strong damping mechanism for propagating compression waves. We conclude that well fluidized granular samples do exhibit local behavior reminiscent of fluids in thermal equilibrium and provide a unique opportunity to study individual and collective particle dynamics in the presence of ultra-strong gradients of density and temperature, which cannot be experimentally achieved in atomic fluids.
Acknowledgements
The authors are grateful to Steven Warr for providing the experimental counterpart of the simulations which was a source of constant inspiration, and to Lyderic Bocquet for suggesting the calculation in the appendix.
Appendix In a stationary state, the conditions of thermal and mechanical equilibria of a slab of width L and thickness dy read LCJQCy) - JQCy + d y ) ) = w C y ) L d y ,
CA.1)
L ( P ( y + dy) - P ( y ) ) =
(A.2)
- p(y)mgL dy,
where JQ(y) denotes the heat flux at altitude y, w ( y ) is the rate of energy dissipation and P ( y ) denotes the pressure at altitude y. These quantities are given by the following constitutive equations [11]: w ( y ) = 7[q(y)T(y)] 3/2 ,
(A.3)
(y) J Q ( y ) = -- ~cdT(Y)dy = - ~[~I(Y)T(y)]I/2 d Tdy
(A.4)
where 7 and t¢ (the thermal conductivity) are phenomenological coefficients; the latter is proportional to ~ and approximately proportional to ~ since at high
373
K. Helal et al./Physica A 240 (1997) 361373
densities the mean free path of the discs is proportional to their mean spacing. For the pressure we use the scaled particle equation-of-state [14]
p(y)ksT(y)
(A.5)
P(Y)= (1 - q(y))2 "
Combination of the relations (A.1)-(A.5) leads to the following set of two coupled differential equations for the packing fraction and temperature profiles t/(y) and T(y):
dq(y)
l
d~-
r(y)
~/(y)[1 - r/(y)] (
l+~-(~S
d2T(y) - ~;~(y)T(y)
dy2
~ +
1
\g(1-~(Y))2+ q(y)
dr(y)~ ay / '
(A.6)
(dT(y)) 2
r(y) 1+ .(y) \ dy J 9 [1 -- r/(y)] 3 dT(y) 2T(y) 1 + q(y) dy
(A.7)
Eqs. (A.6) and (A.7) are solved numerically subject to boundary conditions at y = 0. These boundary conditions are provided by imposing the values of the density, the granular temperature and its gradient (or equivalently the heat flux) within the lowest horizontal slice of width 2a, as determined in the MD simulations. Since the ratio 7/~ of phenomenological coefficients is unknown, it was determined by a least-squares fit of the theoretical temperature profile to the simulation data. The resulting profiles T(y) and r/(y) (the latter did not involve any fitting procedure) shown in Fig. 5, are in excellent agreement with the MD curves. It is worth mentioning that the minimum in the theoretical temperature profile vanishes above a threshold value of the ratio 7/~.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
E. C16ment, J. Rajchenbach, Europhys. Lett. 16 (1991) 133. S. Luding, H.J. Hermann, A. Blumen, Phys. Rev. E 50 (1994) 3100. A. Goldstein, M. Shapiro, L. Moldavsky, M. Fischman, J. Fluid. Mech. 287 (1995) 349. S. Luding, Phys. Rev. E 52 (1995) 4442. Y. Taguchi, H. Takayasu, Europhys. Lett. 30 (1995) 499. S. McNamara, J.L. Barrat, to be published. S. Warr, G.T. Jacques, J.M Huntley, Powder. Tech. 81 (1994) 41. S. Wart, J.M. Huntley, G.T. Jacques, Phys. Rev. E 52 (1995) 5583. S. Warr, J.P. Hansen, to be published. I. Goldhirsch, G. Zanetti, Phys. Rev. Lett. 70 (1993) 1619; S. McNamara, W.R. Young, Phys. Rev. E 50 (1994) R28. [ll] J.T, Jenkins, S.B. Savage, J. Fluid Mech. 130 (1983) 187; P.K. Haft, J. Fluid Mech. 134 (1983) 401. [12] See e.g.J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, 2nd edition, Academic Press, London, 1986.
[13] M.P. Allen, D.J. Tildesley, The computer simulations of liquids, Clarendon Press, Oxford, 1987. [14] M. Baus, J.L. Color, J. Phys. C 19 (1986) L643.
Physica A 240 (1997) 361-373
Local fluctuations in a fluidized granular medium K. H e l a l , T. B i b e n * ,
J.P. H a n s e n
Laboratoire de Physique Ecole Normale SupOrieure de Lyon, 69364 Lyon cedex 07, France 1
Abstract
Molecular dynamics simulations have been carried out on a two-dimensional model granular material made up of inelastically colliding discs subjected to gravity. The granular system is fluidized by a constant energy input from a vibrating base, and a stationary non-equilibrium state is achieved. Density and granular temperature profiles are determined; while the former goes through a maximum at intermediate altitudes, the latter exhibits a minimum at a higher altitude, beyond which the granular temperature is found to increase. Local density fluctuations are characterized by the time-displaced density autocorrelation function. Longitudinal compression (sound) waves propagating in the horizontal direction are found to be strongly damped, while they are overdamped in the vertical direction at the longest accessible wavelength. P A C S : 46.10.--z; 05.60.+w; 05.40.÷j
1. Introduction
A granular material, like a sandpile, may be set in motion under the action of gravity, as in a hopper, or of externally imposed vibrations. Under appropriate conditions, the applied vibrations may induce surface or bulk fluidization of the granular medium, which then exhibits a behavior reminiscent of fluid flow. At the "microscopic" level, the grains undergo inelastic collisions, and a stationary non-equilibrium state is achieved when the dissipated energy is exacty compensated by the continuous energy input from a vibrating base. A number o f recent experiments and numerical simulations carried out on two-dimensional model granular media have examined global and collective properties o f such flows, like the granular temperature, the mean height of the centerof-mass, the propagation of compression waves, or the energy flux through the medium [ 1-6]. Use o f a high-speed camera allows the motion of individual grains confined between vertical glass plates to be monitored as a function o f time [7,8], and very * Corresponding author. E-mail: [email protected]. I Unit6 de Recherche Associre 1325 du CNRS. 0378-4371/97/$17.00 Copyright @ 1997 Elsevier Science B.V. All rights reserved PI[ S0378-4371 (97)00 159-3
362
K. Helal et al./Physica A 240 (1997) 361-373
recently this technique has been used for a quantitative characterization of local density fluctuations in a two-dimensional, vibrated array of steel beads via time-displaced correlation functions familiar from the study of molecular motions in fluids in thermodynamic equilibrium [9]. Some striking analogies between the relaxation of density fluctuations in the fluidized granular medium and in an atomic fluid, but also some significant differences, were discovered by the experiments reported in Ref. [9]. In this paper we present results of molecular dynamics (MD) simulations [13] of a system of inelastically colliding discs subjected to external vibrations and gravity; these simulations were carried out in close connection with the experiments of Ref. [9]. The main objective is to gain a more complete understanding of density and energy fluctuations in a simple model which captures the essential features of the experimental system, rather than to seek a faithful reproduction of the experimental data. In particular, the model under consideration is free of the complications which may arise from the friction experienced by the beads on the confining glass plates. The MD simulations also allow an easier control of the sample-size-dependence of the results as well as better statistics.
2. Model and simulations
We consider a system of N impenetrable discs of diameter tr and mass m confined to the x - y plane. The centers of the discs are confined to a periodically repeated slab of width L in the horizontal (x-direction), and to altitudes y > 0 in the vertical direction (see Fig. 1). The discs are subjected to the constant vertical gravitational force - m y , and suffer instantaneous collisions when the centres of a pair are exactly separated by a distance cr. The pre-collisional and post-collisional (primed) velocities of discs labeled
\_~
\
/
"i
"
s'
= I
..:
/ -
,÷-1
' _.
'
-r--7--T-7"- ~-'-'~7~--'
[
T 1"
.... ,
I
--r--
:+ ,]
k
u
Fig. 1. S c h e m a t i c r e p r e s e n t a t i o n o f the r e c t a n g u l a r o b s e r v a t i o n w i n d o w s L × ( L / 4 ) a n d ( L / 2 ) x ( L / 4 ) u s e d in the simulations. T h e discs are not s h o w n to scale; in practice L = 40tr.
K, Helal et al. I Physica A 240 (1997) 361~73
363
1 and 2 are related by -+-r)[tJl2" (Vl - - I J 2 ) ] R I 2 , • (vl - 1~2)]/112 ,
U'1 =
U 1 -- ½ ( l
U~ =
U2 q- 1 ( 1 -}- r ) [ h l 2
(1)
where hi2 is the unit vector along the direction joining the centres of the colliding discs and r is the coefficient of restitution. Elastic collisions are recovered when r = 1, while in model granular systems, r is typically of the order of 0.95 [7]. The collision rules (1) conserve momentum, but kinetic energy is dissipated when r < 1. This dissipation leads to a rapid "cooling" and "clumping" of grains with given initial velocities [10], but may be compensated by injecting kinetic energy at the y = 0 base through an external vibrator, resulting in a stationary non-equilibrium state. The vertical vibrations imposed in real experiments are sinusoidal in time, but the resulting disc-base collision kinematics are not easily reconciled with the "event-driven" kinematics of binary collisions between discs. For the sake of simplicity, we have modeled the base vibrations and resulting energy input via either of the following two schemes: (a) The base is fixed at y = 0 , but the discs undergo inelastic collisions with the base characterized by a coefficient of restitution R > 1; under these conditions the post-collisional velocity components are vx' = vx; ty" = - R v v , resulting in a relative gain in kinetic energy Ae¢/ec = R 2 - 1 > 0 of a colliding disc. (b) The sinusoidal variation of the base altitude is replaced by an asymmetric sawtooth characterized by its amplitude A and its period to: y ( t ) = A [ t / t o ] - A/2, where Ix] denotes the fractional part of x. Under these conditions a disc gains kinetic energy of at least 2mA2to 2 upon each collision with the base. The simulated sample involved N = 400 discs, corresponding to 10 horizontal layers of 40 discs, comparable to the experimental situation explored by Warr et al. [7-9]. Most simulations were carried out for r = 0.97. A well fluidized stationary state was then obtained with R = 1.5 or ~ = A / ( g t 2 )~_ 5-10. Results obtained with the two-energy input schemes were very similar for comparable rates of base-particle energy transfer; most results reported below were obtained from simulations using energy input scheme (a). Starting from some reasonable initial condition, a typical MD run involved some 10 7 collisions between discs, and statistical averages were taken over the resulting phase space trajectory, beyond an initial "equilibration" period.
3. Density and granular temperature profiles The granular temperature T is defined by the mean-square velocity of the discs according to m
2
r = 5 (v )
(2)
and it is convenient to use the quantity r ' - - kBT - - 1 (1)2)
m
2
(3)
K. Helal et al./Physica A 240 (1997) 361-373
364 0.50
e - - - - o r=O.95;R=l.7 m . - - ~ r=O.97;R=1.5
,o,,o_o. e, 0.40 tO
/,
~,
~f
-'-
-'- r=O.97;vibrator
"\
]
I
0.30 :3 ~
.~_
-'~ 0.20 ¢.~ 0.10
0.00 0.00
10.00
20.00
30.00
40.00
y/o Fig. 2. P a c k i n g fraction profiles r/(y) as a function o f altitude y/a. Circles a n d dash-dotted curve: r = 0.95 and R = 1.7; squares a n d d a s h e d curve: r = 0.97 a n d R = 1.5; triangles a n d full curve: r = 0.97 and saw-tooth vibrator, with amplitude A = 2 / 3 6 a n d period to = 0.02s.
expressing the granular temperature in units of velocity squared. A properly reduced temperature, T* = -TI -,
9a
(4)
characterizes the ratio of the mean kinetic energy of the discs to their gravitational potential energy. The local density, or density profile, p(y) is defined such that p(y)L dy is the mean number of discs in a horizontal strip of area L dy, at the altitude y. It is convenient to define the reduced density profile p*(y)= p ( y ) a 2, and the local packing fraction, q ( y ) = ¼np*(y). Since the granular fluid is "heated" from below, through the energy input from the base, the granular temperature varies with altitude y. We have determined a "temperature profile" T*(y), by averaging the kinetic energy of discs contained in horizontal strips of thickness Ay=2a. Finally, we have computed the rate of energy dissipated per unit mass, w(y), via the loss of kinetic energy at each inelastic collision between discs, as a function of the altitude y. Typical examples of the three profiles r/(y), T*(y) and w(y) are shown in Figs. 2-4, and for two pairs of restitution coefficients (r = 0.97; R = 1.5) and (r = 0.95; R = 1.7). In both cases, the density profile goes through a maximum at an altitude Ym > 10a, in qualitative agreement with the experimental findings of Warr et al. [8]. In the more inelastic case (r =0.95), fluidization is less pronounced, i.e. the discs tend to concentrate at lower altitudes, where the local packing fractions are significantly enhanced relative to the less dissipative system (r = 0.97). More surprisingly, the temperature profile T*(y) is found to go through a minimum at an altitude well above Ym, and to increase at the higher altitudes, where the local density is very low. A faint indication of this unexpected behavior is found in one of the granular temperature profiles measured by Wan" et al.
K. Helal et al. IPhysica A 240 (1997) 361 373
365
25.0 o---~ r=0.95;R=1.7 "~----'~-~rr=00:97;;Rb~a?o r
20.0 \ ~ ~ , , ~
o
15.0
E o
10.0
to'} 5.0 ~"O,,
/O
o---o..e.e..e_..o_o_e_o_.e.-e-e--e-'~
0.0 0.0
L
10.0
20.0 y/~
L
30.0
40.0
Fig. 3. Granular temperature profiles T*(y) versus altitude y/a. The symbols have the same meanings as in Fig. 2.
1.0 ~--. 0.9 e r ~ 0.8 ~ ~\'~,\
.
.
. I" ,~ r=0.97;vibratorI Io---4) r=0.95;a=1.7 I
0.6
0.4
0.1 t 0.0 0.0
, 10.0
"13,, . ~ . , , ~ ~ 20.0 y/o
30.0
~ A 40.0
Fig. 4. Rate of energy dissipation, w(y), normalized by its maximum value Wmax, as a function of altitude. The symbols have the same meanings as in Fig. 2; the results for r = 0.97; R = 1.5 are not shown for visual clarity.
K. Helal et al./Physica A 240 (1997) 361-373
366
19.5 ~
.1¢
1-
14.5
0.2
\
E -~\ e'-
~
ca" 0.0
0.0
9.5
4.5
0.0
~
10.0
20.0
20.0
y/c
4( .0
30.0
40.0
Fig. 5. Comparison between the theoretical temperature and packing fraction profiles, based on the phenomenologicalcalculationof the Appendix (full curves) and the MD data (circles) for r = 0.97; R = 1.5. [8]. A phenomenological calculation based on hydrodynamic equations for dissipative granular media [11], reproduces the observed behavior of the density and temperature profiles very well, as shown in Fig. 5. Details of this calculation are given in the Appendix. The energy dissipation profiles shown in Fig. 4 decay to zero after a slight initial rise involving the first few particle layers above the base. Once more the decay is significantly faster in the case of the smaller coefficient of restitution. The concept of a local granular tempera~re, playing a role similar to the equilibrium temperature in atomic or molecular systems, is meaningful only if the velocity distribution is strictly Gaussian (or Maxwellian), i.e. ~b(v) °c exp ( - 2~7) exp
-2T']
(5)
To test this assumption, we performed statistical averages of the distributions of the horizontal (Vx) and vertical (Vy) components of the particle velocities over strips of thickness A y = 2 a at several altitudes. Typical examples are shown in Fig. 6. They show a significant difference between the distributions of the horizontal and vertical velocity components. While the former are symmetric around vx = 0 and Gaussian within statistical error at all altitudes, the distributions of the vertical velocities are seen to be close to Gaussian only at intermediate altitudes, while they are significantly skewed to the left (negative Vy) near the base, and to the right (positive Vy) at the highest altitudes. The width of the distributions, proportional to x/T-7, first shrinks with altitude, but broadens again beyond y-~ 20a, in agreement with the behavior of the
K. Helal et al./Physica A 240 (1997) 36l 373
50.0
•
•
i
367 !
(A) 50.0
25.0
0.0 -5.(
0.0
5.0
0.0 -5.0
0.0
5.0
i
(c)
(D)
50.0
Vv
5.0
0.0 -5.0
0.0
5.0
0.0 -5.0
0.0
i.O
Fig. 6. Velocity distribution functions (O(Vx) (crosses) and ~b(Vv) (triangles) at altitudes y/a= 1 (A), I1 (B), 23 (C) and 39 (D), from an MD simulation with r = 0 . 9 7 and R = 1.5. Velocities are in units where L = 1 and 9 : 10. The full curves are the best Gaussian fits to 4~(v~).
temperature profile shown in Fig. 3. Inspection of Fig. 6 leads us to conclude that the concept of a local granular temperature is reasonably well founded.
4. Local density fluctuations Microscopic density fluctuations in an atomic fluid in thermal equilibrium are routinely characterized by the time-displaced autocorrelation function of the Fourier components of the microscopic density, conventionally referred to as the "intermediate scattering function" in connection with inelastic neutron scattering experiments [t2]. For a closed or periodic system containing N atoms, the density autocorrelation function is defined as F ( k , t) = 1 ( p k ( t ) p - k ( O ) ) ,
(6)
where the brackets denote a thermal average and the Fourier components of the microscopic density are N
p k ( t ) = ~-~e ik~'(') , j=|
(7)
K. Helalet al./Physica A 240 (1997) 361-373
368
if rj(t) denotes the position of the jth atom at time t. Experimental or numerical measurements of F(k, t) provide a powerful tool to investigate collective particle motion on a scale 2 = 2n/k. In particular, propagating longitudinal sound (or compression) waves show up as damped oscillations in F(k, t) for sufficiently long wavelengths (2>>a). At shorter wavelengths (2 _~ a), density fluctuations are found to decay monotonously to zero. In a very recent experiment [9], F(k,t) has been measured in a 2d vibrated granular medium with the help of a high-speed camera, which allows the precise determination of individual grain positions within a rectangular window for a large number of successive instantaneous configurations. We have computed F(k, t) in our MD simulations, for several wave vectors k, and under physical conditions similar to those in Ref. [9]. In the simulations we used windows of area Lx x Ly. Lx was taken to be the total horizontal width of the system, with periodic boundary conditions at ±L/2, or one-half that width, with open boundaries at ±L/4. In the vertical direction, the window extended from y = 7.5a to y = 17.5a above the base, i.e. Ly 10or. Referring back to the density profile in Fig. 2, for ( r = 0 . 9 7 ; R = 1.5), it may be seen that this choice was made to study a sample as uniform as possible in the vertical direction. The accessible wave vectors k are of the form =
k
z
(nx ny ~ 27r\Lx L y J ' - - ,
- -
(8)
where nx, ny are integers. There are some fundamental differences between the density autocorrelation functions in a thermal fluid and in a fluidized granular medium: (a) The statistical average is to be taken over initial conditions which belong to an equilibrium ensemble for fluids, and to an ensemble of stationary non-equilibrium states in the case of fluidized granular matter; in both cases stationarity implies time translation invariance of the correlation functions. (b) In an isotropic fluid, rotational invariance implies that F(k,t) depends only on k = Ikl while translational invariance implies that only Fourier components of opposite wave vectors have non-vanishing correlations. Both invariances are broken in granular matter due to the preferred (vertical) direction, and the vertical inhomogeneity imposed by the base vibrations and gravity. We have nevertheless restricted our investigations to ( k , - k ) pairs and considered density fluctuations with horizontal and vertical wave vectors separately. (c) For equilibrium fluids, one normally considers closed or periodic systems, such that the number of atoms contributing to the density fluctuations (7) is fixed. In the present simulations, and in the experiments of Ref. [9], the granular sample within the window under observation forms an open system. This means that N, and the identity of the discs or beads within the window, fluctuate in time. Consequently, N is replaced by its mean value (N) in definition (6) (other conventions could be adopted). Moreover, the lack of space inversion (k ~ - k ) symmetry, combined with the open nature of the system, imply that F(k,t) is in general no longer a real function, except for t = 0, or for horizontal k-vectors, since the x ~ - x inversion symmetry is preserved. In the t = 0 limit, F(k,t) reduces to the static
K. Helal et al./ Physica A 240 (1997) 361-373
369
2.0
"5
o
,=
1.0
09 0
[
¢/)
~
0
0.0
Baus & Colot 0 Full window x Half window
I
I
i
I
I
2.0
4.0
6.0
8.0
10.0
12.0
q=kc Fig. 7. Static structure factor S(q) versus reduced wave number q = ka determined from particle coordinates in the observation windows shown in Fig. 1 and described in the text; the MD simulations are for r = 0.97; R = 1.5 and the mean packing fraction in the window is ~7= 0.462. Circles and crosses correspond to the full and half-windows, respectively; data points are for wave vectors of variable direction, stressing the isotropy of the structure factor; the full curve is the structure factor for a homogeneous fluid of hard discs at t / = q [14].
structure factor 1
S(k) = F(k, t = O) = 7-~ (PkP-k)
(9)
while in the long-time limit 1
2
t l i m F ( k , t ) = -~-~[(p~)] .
(10)
For horizontal wave vectors, this limit vanishes, since the mean density is constant in the x-direction, while for vertical wave vectors, (Pk) is the Fourier transform of the vertical density profile: jfl, "¢2
(Pk) =
eikYp(y) dy,
( 11 )
where Yl and Y2 are the lower and upper altitudes of the simulation window. Since p(y) is not uniform, (Pk) is non-zero, and the resulting limit (10) of F(ky, t) is non-zero. A few examples of results from the MD simulations are shown in Figs. 7-9. In Fig. 7 the computed static structure factor S(k) is compared to that of a homogeneous
K. Helal et aL /Physica A 240 (1997) 361-373
370
1.0
~" 0.05 (/) v
09
0.5
0.00
v'
,
i
I
0.0
i
ii
~
-~''-~- -
25.0
ii
50.
(O.tE
0'00.0
20,0
40,0
60.0
80.0
100.0
t/t E
Fig. 8. Normalized density autocorrelation function F(k,t) versus time (in units of the Enskog collision time te), for k = (27r/L, 0), and the full observation window in Fig. 1 (MD run with r = 0.97; R = 1.5). The corresponding spectrum S(k, 02) is plotted versus ~o • tE in the insert.
fluid of elastic hard discs in thermal equilibrium, as parametrized by Baus and Colot [14]. The agreement is seen to be excellent. Significant discrepancies were found in a similar comparison between the elastic hard disc results and experimental data [9]; these discrepancies are probably due to the fact that the granular medium in the experimental window was much more inhomogeneous than the corresponding simulated system. Both simulation and experimental data are surprisingly isotropic, i.e. for given [k], S(k) is practically independent of the orientation of the wave vector. Fig. 8 shows the real part of F(k, t) for the smallest horizontal wave vector, k---(2rt/L, 0). A small amplitude oscillation is seen to be superimposed on the slow decay of the correlation function, characteristic of a strongly damped longitudinal collective mode, which we tentatively identify with a compression (or "sound") mode propagating through the model granular medium. The spectrum of the correlation function obtained by Fourier transformation with respect to time, exhibits a very broad resonance at an angular frequency o2"~ 27.5 rads - l (in unit where L = 1) corresponding to a sound velocity c = og/k _~ 1.9 m s - 1 for particles of diameter tr = 5 mm roughly compatible with that predicted from the adiabatic compressibility of a homogeneous fluid of elastic hard discs with density equal to the mean density of the inelastic system considered here. A MD simulation carried out for elastic discs shows that the oscillations in F(k,t) are much more pronounced in that case, leading to a much sharper resonance in the density fluctuation spectrum. Fig. 9 compares the real parts of F(k,t) computed for horizontal and vertical wave vectors of the same modulus. As foreseen in our earlier discussion, the correlation function for the horizontal k decays to zero, while that for the vertical k goes to a non-zero limit, which agrees within statistical errors with the value computed from
371
K. Helal et al./Physica A 240 (1997) 361-373 1.0
1.0 0.9
v
iJ_
0,9
(A)
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0.0
0.0 -0.1
(B)
w
0.0
10.0 t/t E
20.0
-0.1
i
0.0
10.0
20.0
t/t E
Fig. 9. Normalized density autocorrelation function F ( k , t ) versus t/tE for: (A) k=(O, 87z/L); (B) k = (8~z/L,O), under conditions identical to those in Fig. 8.
Eq. (10). Data obtained for several wave vectors show that the rate of relaxation initially increases with Ikl, as in the case of fluids in thermal equilibrium. For all correlation functions which we have investigated, the imaginary part is very small at all times, remaining below the statistical noise level.
5. Conclusion The simulations described in this paper point to some strong similarities between the collective behavior of fluidized granular matter on a local scale (i.e. on the scale of a few grain diameters) and the corresponding properties of equilibrium fluids on a microscopic scale. A similar conclusion had been drawn on the basis of experiments reported in Ref. [9], and semi-quantitative agreement is found between the experimental and simulation data for comparable physical conditions. Both sets of results are very sensitive to the amplitude and frequency of the external vibrations; it is a matter of delicate balance between system size, rate of energy input and rate of bulk-energy dissipation to achieve a stationary fluidized granular state. In particular, the present work shows that the static and dynamical properties of a system of inelastic discs are very sensitive to the coefficient of restitution r, which has, in particular, a strong incidence on the density profile. Preliminary simulations for r = 0.92 point to the formation of a slowly varying inhomogeneity (long-lived "clumps") in the horizontal direction. The key findings of the present work may be summarized as follows: (a) The granular temperature profile goes through a minimum at intermediate altitudes, and increases in the low-density region at higher altitudes. This unexpected behavior is well reproduced by a hydrodynamic calculation (Appendix).
372
K. Helal et al. / Physica .4 240 (1997) 361-373
(b) The static structure factor measured over a range of altitudes where the system is nearly homogeneous (i.e. around the maximum in the density profile) is in perfect agreement with the known structure factor of a homogeneous fluid of hard discs in thermodynamic equilibrium. This indicates that the local structure in a stationnary non-equilibrium state of fluidized granular matter strongly resembles that of fluids in thermal equilibrium. (c) Longitudinal compression (or "sound") waves are strongly damped even for the longest accessible wavelength, compared to the corresponding behavior in an equivalent system of elastic discs; a small degree of inelasticity (r = 0.97 in the present calculations), coupled with the vertical inhomogeneity, is thus found to provide a strong damping mechanism for propagating compression waves. We conclude that well fluidized granular samples do exhibit local behavior reminiscent of fluids in thermal equilibrium and provide a unique opportunity to study individual and collective particle dynamics in the presence of ultra-strong gradients of density and temperature, which cannot be experimentally achieved in atomic fluids.
Acknowledgements
The authors are grateful to Steven Warr for providing the experimental counterpart of the simulations which was a source of constant inspiration, and to Lyderic Bocquet for suggesting the calculation in the appendix.
Appendix In a stationary state, the conditions of thermal and mechanical equilibria of a slab of width L and thickness dy read LCJQCy) - JQCy + d y ) ) = w C y ) L d y ,
CA.1)
L ( P ( y + dy) - P ( y ) ) =
(A.2)
- p(y)mgL dy,
where JQ(y) denotes the heat flux at altitude y, w ( y ) is the rate of energy dissipation and P ( y ) denotes the pressure at altitude y. These quantities are given by the following constitutive equations [11]: w ( y ) = 7[q(y)T(y)] 3/2 ,
(A.3)
(y) J Q ( y ) = -- ~cdT(Y)dy = - ~[~I(Y)T(y)]I/2 d Tdy
(A.4)
where 7 and t¢ (the thermal conductivity) are phenomenological coefficients; the latter is proportional to ~ and approximately proportional to ~ since at high
373
K. Helal et al./Physica A 240 (1997) 361373
densities the mean free path of the discs is proportional to their mean spacing. For the pressure we use the scaled particle equation-of-state [14]
p(y)ksT(y)
(A.5)
P(Y)= (1 - q(y))2 "
Combination of the relations (A.1)-(A.5) leads to the following set of two coupled differential equations for the packing fraction and temperature profiles t/(y) and T(y):
dq(y)
l
d~-
r(y)
~/(y)[1 - r/(y)] (
l+~-(~S
d2T(y) - ~;~(y)T(y)
dy2
~ +
1
\g(1-~(Y))2+ q(y)
dr(y)~ ay / '
(A.6)
(dT(y)) 2
r(y) 1+ .(y) \ dy J 9 [1 -- r/(y)] 3 dT(y) 2T(y) 1 + q(y) dy
(A.7)
Eqs. (A.6) and (A.7) are solved numerically subject to boundary conditions at y = 0. These boundary conditions are provided by imposing the values of the density, the granular temperature and its gradient (or equivalently the heat flux) within the lowest horizontal slice of width 2a, as determined in the MD simulations. Since the ratio 7/~ of phenomenological coefficients is unknown, it was determined by a least-squares fit of the theoretical temperature profile to the simulation data. The resulting profiles T(y) and r/(y) (the latter did not involve any fitting procedure) shown in Fig. 5, are in excellent agreement with the MD curves. It is worth mentioning that the minimum in the theoretical temperature profile vanishes above a threshold value of the ratio 7/~.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
E. C16ment, J. Rajchenbach, Europhys. Lett. 16 (1991) 133. S. Luding, H.J. Hermann, A. Blumen, Phys. Rev. E 50 (1994) 3100. A. Goldstein, M. Shapiro, L. Moldavsky, M. Fischman, J. Fluid. Mech. 287 (1995) 349. S. Luding, Phys. Rev. E 52 (1995) 4442. Y. Taguchi, H. Takayasu, Europhys. Lett. 30 (1995) 499. S. McNamara, J.L. Barrat, to be published. S. Warr, G.T. Jacques, J.M Huntley, Powder. Tech. 81 (1994) 41. S. Wart, J.M. Huntley, G.T. Jacques, Phys. Rev. E 52 (1995) 5583. S. Warr, J.P. Hansen, to be published. I. Goldhirsch, G. Zanetti, Phys. Rev. Lett. 70 (1993) 1619; S. McNamara, W.R. Young, Phys. Rev. E 50 (1994) R28. [ll] J.T, Jenkins, S.B. Savage, J. Fluid Mech. 130 (1983) 187; P.K. Haft, J. Fluid Mech. 134 (1983) 401. [12] See e.g.J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, 2nd edition, Academic Press, London, 1986.
[13] M.P. Allen, D.J. Tildesley, The computer simulations of liquids, Clarendon Press, Oxford, 1987. [14] M. Baus, J.L. Color, J. Phys. C 19 (1986) L643.