数理解析研究所講究録 1413 巻 2005 年 248-254

248
Langevin equation with Coulomb friction -Simulation and theory of granular particles
under the external vibration
Hisao Hayakawa and Atsushi Kawarada
Department of Physics, Yoshida-South Campus, Kyoto University, Kyoto 606-8501,JapanThe simulation of granular particles in aquasi two-dimensional container under the vertical vibration as an experimental accessible model for granular gases is performed. The velocity distributionfunction obeys an exponential-like function during the vibration and deviates from the exponentialfunction in free-cooling states. It is confirmed that this exponential distribution function is produced
by Coulomb ffiction force. ALangevin equation with Coulomb’s friction is proposed to describe themotion of such the system. Through the analysis of the corresponding Fokker-Planck equation, we
have obtained the steady velocity distribution function under the influence of the external field.
I.

INTRODUCTION

Granular physics has been achallenging subject of statistical physics since the rediscovery of significant nature
An assembly of grains has so strong fluctuations in theor the earlyin granular materials in the lateconfiguration, the force and the motion that mean-field theories cannot be used in most of situations. Atypicalexample of the strong fluctuation appears in the force distribution for astatic granular assembly piled by the gravity.[2]The force propagates along force-chains and the distribution function of the magnitude of the force on the bottom ofacontainer does not obey Gaussian but an exponential function.
Such the strong fluctuation coming from non-Gaussian properties should be relevant even in the dynamics of granularassemblies. However, there are not so many systematic researches to focus on the statistical distribution functions inthe steady state, because real and numerical experiments report no unified results, $i.e.$ , velocity distribution functions(VDF) obey Gaussian-like with the exponential tai1[5-12], functions from Gaussian to the exponential depending on
the density $[13, 14]$ , the stretched exponentia1[15-17] and even power-law functions $[18, 19]$ depending on situations.
As in the standard statistical mechanics, an assembly of grains in agas phase is an idealistic situation to study
what the proper statistical weight is. Approximate granular gases can be obtained by rapid granular flow on an
inclined surface $[22, 23]$ , gas-particles mixtures[25] and the extereal vibration.[ll, 14-16, 18] However, these systems
have defects, because (i) the boundary and the gravity effects are so strong in the rapid granular flows, (ii) thehydrodynamic interaction between particles are so complicated, and (iii) adense cluster appears in the vibrating
experiment. Therefore, it is difficult to achieve free-cooling gases without the effect of gravity in experiments.[15]In this report, to remove such the difficulties, first, we propose an experimental accessible situation to produce
granular gases. Second, we demonstrate that VDF in both dense and dilute granular gases under the vertical vibration
can obey an exponential-like function when Coulomb friction is important through our simulation based on the distinct
element method (DEM). This exponential VDF disappears immediately after vibration is stopped, $i.e$ . the grains
are in afree-cooling process. To explain these results, third, we introduce aphenomenological Langevin equation
to describe the motion of particles and explain the mechanism to appear the exponential VDF. We also investigate
mathematical properties of Langevin equation with Coulomb friction in details.
The organization of this paper is as follows. In the next section, we introduce our setup in the simulation. In
section 3, we surnrnarize the results of our simulation. In section 4, we obtain the steady VDF under the influenceof the steady external field based on Langevin equation with Coulomb friction. In section 5we discuss our result
through the comparison of our result with the simulation of granular particles in aquasi-tw0-dimensional box. In the
final section, we conclude our results. We note that this report is basically joint one from two papers $[26, 27]$
$80’ \mathrm{s}$II.

$90’ \mathrm{s}.[1]$THE SETUP OF OUR SIMULATION

We use athree-dimensional DEM for monodispersed spherical particles.[28] Instead of using the Hertzian contact
force, we adopt the linear spring model to represent repulsion of contacted spheres. We also include the rotation of
spheres and the Coulomb slip for tangential contact. We simulate amonodisperse system with particle diameter .
Because of the limitation of the space we skip the detailed description of model including the choice of parameters.
We can see the parameters and the details of model in ref.[26].
We focus on the following situation: Particles are confined in aquasi tw0-dimensional container in which the heightis 1.Sd and the horizontal plane is asquare (Fig. 1). In this system, particles cannot compose amultilayer configurationin the vertical direction, but the particles have vertical velocities in their motion. The mobile particles are randoml$d$$\mathrm{y}$

243
with

$\mathrm{d}\mathrm{e}\mathrm{P}’|^{\mathrm{b}}\prime \mathrm{o}_{\mathrm{O}^{\mathrm{O}}}\mathrm{o}\mathrm{o}_{\mathrm{o}_{\mathrm{o}^{\mathrm{O}}\mathrm{o}}}^{\mathrm{o}\mathrm{o}_{\mathrm{O}}^{\mathrm{O}}}\mathrm{o}\mathrm{o}\mathrm{o}^{\mathrm{O}}\mathrm{o}^{\mathrm{o}^{\mathrm{O}}}\mathrm{o}_{\mathrm{O}\mathrm{O}}^{\mathrm{O}\mathrm{O}\mathrm{O}}$Top$\mathrm{V}\dot{1}8\mathrm{W}$

FIG. 1: The schematic side view (left) and the top view (right) of our setup: The random scatters are fixed spherical particles
on the top board. The diameter of fixed particles on the top plate (white ones) is uniformly random between $0.6d$ and O.Sd, andtop board. To avoid stack of particles at the corners we introduce four fixed particles
their centers are located $0.15d$ above(black ones in the figure) at the corners.
$\mathrm{t}\mathrm{l}\iota \mathrm{e}$$\mathrm{c}$in steady states under the vertical vibration. Here represents the area fraction. For $=7.21\%$
FIG. 2: The scaled VDF
we use 10,000 particles, and use 1,000 particles for other cases. It should be noted that Fig.2 in ref.[26] use the data of 3,000$\Psi$$\overline{f}(c)$$\Psi$

particles.

scattered by fixed particles on the top of the container if the container is vertically vibrated. The number of fixed
scatters is 3.5 times of that of the mobile particles. The vibration is driven by a sinusoidal force whose amplitude
and the frequency are given by $A=1.2d$ and $f=0.5\sqrt{g}/d$ , respectively. Thus, the acceleration amplitude becomes$\Gamma=A(2\pi f)^{2}/g=11.5$ . If the vibration is stopped, the particles are moving with the rotation on the bottom plate,then the system can simulate a tw0-dimensional free-cooling granular gas.
The simulation starts ffom a stationary state of particles on the bottom plate. The particles gain the kinetic energy
from the external oscillation and the system reaches a steady state in the balance between collisional dissipations andthe gain of the energy from the external force. Typically, the system reaches the statistical steady state after 25 cyclesof the oscillation. We should note that the kinetic energy of particles is a little oscillated depending on the phase of
the extemal oscillation after the system reaches the statistical steady state.
III.

THE RESULTS OF OUR SIMULATIONIn this section, we summarize the reslllt of our simulation.
In the steady state, we confirm that the density is uniform and there is no long-range correlation in contrast to
granular gases in free-cooling states. Thus, the system does not have any systematic flows and any definite clusters.
We also check that the effects of side boundaries can be neglected. This is because particles always hit the top wallor the bottom wall during the vibration. After we stop the vertical oscillation, the correlation grows with time as thefree-cooling process proceeds.
The most characteristic quantity for this gas system is VDF. In the steady state, the vertical component of VDF
has double peaks where each peak corresponds to a lifting or a falling process, while the horizontal VDF has a single
peak. For later discussion, we only use the horizontal VDFs for the analysis. We plot the scaled horizontal VDF$\tilde{f}(c)$250$\mathrm{c}$FIG. 3: The scaled VDF for several situations. In the inset, ’no tangential’ and ’no friction’ correspond to the results ofsimulations without } tangential force
friction force, respectively, where
systems are in the steady state andinclude 1,000 particles with the area fraction 0.0649. Here, ’cooling state’ and ’undulation’ indicate the results of simulation of
the cooling process of 10,000 particles and the undulation, respective . Note that Fig.3 incontains the data for 3,000
particles.
$\mathrm{a}\mathrm{r}\iota \mathrm{d}\mathrm{t}\mathrm{t}\iota \mathrm{e}$$\mathrm{t}$$\mathrm{t}\mathrm{l}\downarrow \mathrm{e}$$\mathrm{l}\mathrm{e}$$\mathrm{r}\mathrm{e}\mathrm{f}.[?

$1\mathrm{y}$]$in Fig.2(1)$f(\mathrm{v},t)=r\iota\tau_{0}’(t)^{-2}\overline{f}(\mathrm{v}/\mathrm{c}_{0}’)$with the density , the average speedwith the granular temperature ,,where
is averaged over the cycles. As in Fig.2,
scaled VDF call be approximated by an exponential function.
fact, the flatness defined by $<c^{4}>/<c^{\sim}.,$is llot far from 6. It should
be noted that tlle flatness with tlle Gaussianis 3 and with the exponential VDF is 6. 111 our simulation.
tlue flatnesses are 8.00, 7.63, 6.89, 6.26, 5.67 alld 5.35 for 1000 particles’ simulation for corresponding area fractions
projected into the horizontal plane 0.0649, 0.100, 0.200, 0.301, 0.400 and 0.501. More de tails results for 3,000 particles) , the dependence of
can be seen in Table I in ref.[26]. In(lellsity is relatively weaker in our
flatnesssituation. For a large system with 10,000 particles with the area fraction 0.0721 has smaller flatness as 6.85. If theexternal oscillation stops, tlle flatness decreases quickly a1lll is saturated at 4.20 for 10,000 particles. As carl be seenin Fig.3,
in tlle cooling process is nearly Gaussian for low energy particles but has all exponential tail for highenergy particles as in the usual gas systems.[5] Note thatarea fraction 0.0721 corresponds to the volume fraction
0.0271.
Let us consider the origin of tlle exponential-like VDF. It is easy to verify that the exponential VDF cannot be
reproduced without tlle existence of Coulomb friction in DEM. In fact, if all component of tangential contact force
is omitted, the flatness becomes 3.43 for 3,000 particles with the area fraction 0.0649. While if we only neglect
effect of Coulomb fficti in, the flatness becomes 3.43 in the same situation. In our system, particles feelthe stroll shear force when particles llit the fixed scatters on the top wall, because
colliding particles contact theside of fixed scatters, general,the directions of their rnotiorl are changed drastically. Thus, the tangential slips
between particles and the fixed scatters are the dominant dissipative processes ill } steady state. Since we specify
the origin of exponential-like VDF, we can understand the weak dependence of VDF on the density. Namely, particles
directly hit the fixed scatters for dilute case, while dense particles collide with each other and rotate without slips
besides the collisions withwalls.
Therefore, the ssence to produce tlle large flatness ill VDF is apparently Coulomb’s friction. Thus, we propose tlle
following Langevin equation to describe horizontal motion of particles:$?\prime_{0}=\sqrt{2T/\pi\iota}$$n$$\tilde{f}(\mathrm{r})$$T$$\mathrm{a}\mathrm{n}\mathrm{d}/|d\mathrm{c}\overline{f}(c)=\int^{1}d\mathrm{c}c^{\underline{)}}.\overline{f}(c)=1$$\mathrm{t}\mathrm{I}_{1}\mathrm{e}$$\ln$$>^{2}=<c^{4}> \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}<\mathrm{r}^{\prime 1}.>=\int d\mathrm{c}c^{n}\tilde{f}(c)$$\backslash \prime \mathrm{D}\prime \mathrm{F}$$\mathrm{i}\lambda \mathrm{d}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{e}$$\mathrm{t}1_{1}\mathrm{e}$$11$$011\mathrm{t}\mathrm{l}\mathrm{l}\mathrm{e}$$\backslash ’\mathrm{D}\mathrm{F}$$\mathrm{t}1_{1}\mathrm{e}$$\mathrm{t}\mathrm{l}\mathrm{l}\mathrm{e}$$()\Pi$$\mathrm{e}(1\cdot(^{\uparrow?})$$\mathrm{t}\mathrm{l}\mathrm{l}\mathrm{e}$$\mathrm{a}|\mathrm{l}\mathrm{d}$$\mathrm{i}_{\grave{1}}\mathrm{l}$$\mathrm{t}$$\downarrow \mathrm{e}$$\mathrm{t}1_{1}\mathrm{e}$$\mathrm{e}$,(2)$\frac{d\mathrm{v}}{dt}=-\zeta\frac{\mathrm{v}}{l}$where$\mathrm{v}$,$v$

are respectively the velocity, and the speed. The friction may be proportional towhich is the random force satisfies
fluctuation-dissipation relation:$=|\mathrm{v}|$colnpollellt of$+\eta_{\backslash }$$\zeta$$l^{\mathit{1}}g$. The$\alpha$$\mathrm{t}1_{1}\mathrm{e}$$\eta$$<\eta_{\alpha}(t)>=0$,$<\eta_{\alpha}(t)\eta_{l\mathrm{J}^{1}}(t’)>=2D\delta_{a\mathrm{d}}\delta(t-t’)$

,

(3)

251

for 1,000 particles. The solid line represents
andFIG. 4: The relation between
can be seen in Fig.4 in ref.[26].
3,000
which
contains until
particles
including the data of$\overline{D}=0.0945\sqrt{\tilde{T}}$$\overline{T}$$\overline{D}$-0.0099. The figure$\sqrt{\overline{T}}\simeq 0.5$$(2)$ with (3) can be converted into
where $D=\gamma\sqrt{T}/3m$ is tlle diffusion coefficient. The Langevin:probability distribution function
Fokker-Planck or Kram ), equation for$\mathrm{t}1_{1}\mathrm{e}$$\mathrm{t}\mathrm{i}\mathrm{q}\iota \mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$$P(\mathrm{x}, \mathrm{u}, t)$$\mathrm{t}\mathrm{l}\mathrm{l}\mathrm{e}$$‘$$\mathrm{r}\mathrm{s}$$\frac{\partial P(\mathrm{u},t)}{\partial t}=\{\zeta\frac{\partial}{\partial \mathrm{v}}\cdot\frac{\mathrm{v}}{\tau}+D\frac{\partial\sim)}{\partial_{11}\sim}.., \}P$

(4)$(\mathrm{v}, t)$’

For steady states, the equation forcall be obtained easily$P$–$P_{st}$is reduced to $(\mathrm{u}/u)P_{st}+\sqrt{T/3m}(d/d\mathrm{v})P_{sl}=0$ for $11\neq 0$ . Its solution.

(5)$P_{st}(\mathrm{v})=2\sqrt{\frac{m}{3T}}\exp[-\sqrt{\frac{3m}{T}}l’]$)VDF.
Thus, the Langevin equation with Colo mb ffiction law obeys theof particles, we evaluate both the diffusion
motionequationthevalidityfor
LaIlgevill
of
new
tlle
To check the. At first,
have confirmed that tlle motion of particles in the horizontal
coefficient $D$ and the friction( -to), whereplane is diffusive. Then,diffusion coefficient is evaluated from tlle relation
we choose as 25 cycles of the oscillation and simulate tlle motion of particles until 75 cycles. For each para meter
ment of VDF. Thus, we obtain the relation between $D$ andofoscillation we cleterlrline from the second
and–0.011, where
speaking, the best fitted relation isas in Fig.4.
$\overline{T}=T/(rngd)$ , but this call be replaced bv-0.0099 for 1,000 particles. If we include the data of
3,000 particles the fitted values are a little different from those presented here, but the tendency is the same. the also
note that the collective vertical motion suppresses for larger , and llo diffusion takes place because of the insufficient
as $=0.164g$ which is independent
call evaluate
kinetic energy for maller . Using The relation betw een andparticles
and the horizontal walls caused
betweenprocess
collisions
isdominant dissipativeof . This is because
believe that the Langevin equation (2) can bevertical direction. Thus,by the collective motion of particles ill
used to describe the motion of particles.
$0\tau \mathrm{l}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}$$\mathrm{e}\mathrm{x}\mathrm{I}$$\mathrm{c}\mathrm{o}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{f}\mathrm{f}\mathrm{i}(_{\mathrm{J}}^{\backslash }\mathrm{i}\mathrm{e}\mathrm{I}1\mathrm{t}$$\mathfrak{n}\cdot \mathrm{e}$$\gamma$$<(\mathrm{r}(t)-\mathrm{r}(t_{0}))\underline{.)}>=4D$$\mathrm{t}1_{1}\mathrm{e}$$t$$t_{0}$$T$$\mathrm{t}\mathrm{l}\mathrm{l}\mathrm{e}$$T$$\mathrm{m}\mathrm{o}$$\tilde{D}=0.0925\overline{T}^{0470}$$\tilde{D}=D/(d^{1/\underline{\cdot)}}g^{3/\cdot)}\sim)$$\mathrm{I}\mathrm{I}\mathrm{o}\mathrm{l}\iota \mathrm{e}\mathrm{s}\mathrm{t}\mathrm{l}\mathrm{y}$$\overline{D}=0.0945\sqrt{\tilde{T}}$

$T$

$T$

$\tilde{D}$$T$$\mathrm{s}$$\tilde{T}\mathrm{v}\prime \mathrm{e}$$\gamma$$\gamma$$\mathrm{t}\mathrm{l}\mathrm{l}\mathrm{e}$$\mathrm{t}1_{1}\mathrm{e}$IV.

$\backslash \iota’\mathrm{e}$MATHEMATICAL ANALYSIS OF LANGEVIN EQUATION WITH COULOMB FRICTION

111 the previous section we suggest that Langevin equation associated with Coulomb friction can be an effectiveof granular particles through their simu lation of granular particles. In this section,
equation a vibratingthe.
investigate mathematical properties of Langevin equation with Coulomb friction following
to our system.
ill eq.(2) to discuss the linear response of
external force
introduce
It is convenient to use dimensionless distribution function. For this purpose, we$\mathrm{v}\cdot \mathrm{e}$$\mathrm{i}\iota 1$$\mathrm{s}\backslash .\cdot \mathrm{s}\mathrm{t}\mathrm{e}\mathrm{l}\mathrm{I}\mathrm{l}$’$\backslash \backslash ’ \mathrm{e}$$\mathrm{r}\mathrm{e}\mathrm{f}.[?

$\mathrm{F}_{ex}$]$$\mathrm{a}\{\mathrm{i}\mathrm{d}$$\mathrm{m}\mathrm{a}3$$\mathrm{F}_{ex}$$1\mathrm{n}\mathrm{a}]’$$P(\mathrm{v}, t)=\mathrm{z}_{0}’(t)^{-d}\overline{P}(\mathrm{c}, t)$,$\mathrm{c}=\mathrm{v}/\mathrm{c}_{0}’(t)$(6)

with the normalization$\int d\mathrm{c}c^{\vee}.)\tilde{P}_{0}(\mathrm{c}, t)=1$

(7)

252
dimension, andis the scaled VDF without. We adopt instead of in eq.(7), because thedetermination of the scaled factor is difficult if we use . In fact, isfunction ofarld cannot be determined
without the complete form of VDF.
Before we proceed to the next section, we give somements
the model. First, Coulomb friction is singularat
as inand (4). Indeed
start and the stop of frictional motions of macroscopic materials are thesingular processes such as the break up and the formation of force networks. In addition, the static friction coefficientis almost always larger than the dynamical friction coefficient. Therefore, we can guess that particles are condensed
in simple setups. Indeed, when we simulate the motion of frictional disks on an inclined surface with denseatrandom scatters, most of particles are trapped by arrays of scatters but VDF for mobile particles see ms to obey an
exponential function in a steady state. Even when there are no condensed particles at, the differentiation of
VDF with respect to in
is discontinuous at. Thus, VDF can be singular at. Here we are onlyinterested in the case ofin this paper.
Let us consider the situation that a particle is moving under the influence of $F_{ex}=rngz$ with the unit vector
parallel to the direction of the external force and Coulomb friction force associated with the random force. We assume
that the system is a tw0-dimensional one, because the particle is typically located on a substrate when Coulomb friction
is important. As mentioned in the previous section is not direction of gravity but the direction of the slope 011 }horizontal plane in the setup of Kawarada and Hayakawa[26]. We are interested in the statistical steady state in thebalance between
ffiction and the external force.
Let be the angle between the direction we consider and , the steady equation forbecomes

where

$d$

is

$\overline{P}_{0}$$\mathrm{t}\mathrm{l}\mathrm{l}\mathrm{e}$$F_{ex}$$\tilde{P}$$\mathrm{c}\mathrm{o}\mathrm{m}$$\mathrm{v}=0$$\tilde{P}$$\tilde{P}_{0}$$T$$F_{ex}$$\mathrm{t}1_{1}\mathrm{e}$$\mathrm{o}\mathrm{r}1$$\mathrm{t}\}_{1}\mathrm{e}$$\mathrm{e}\mathrm{q}\mathrm{s}.(2)$$\mathrm{v}=0$$\mathrm{v}=0$$(^{??})$$\mathrm{v}$$\mathrm{v}=0$$\mathrm{v}=0$$\mathrm{v}\neq 0$$\hat{z}$$\hat{z}$$\mathrm{t}$$\mathrm{l}\mathrm{e}$$\mathrm{t}\mathrm{t}_{1}\mathrm{e}$$\theta$$\hat{z}$$\mathrm{v}\neq 0$$g[ \cos\theta\frac{\partial P}{\partial\iota}, +\frac{\mathrm{s}\mathrm{i}\mathrm{r}1^{2}\theta}{l},\frac{\partial P}{\partial\cos\theta}]=\zeta\{\frac{P}{\mathrm{t}}+\frac{\partial P}{\partial\iota}\}+" D[\frac{\partial^{2}}{\partial_{\mathrm{t}^{2}}}+\frac{1}{\mathrm{t}}\frac{\partial}{\partial}+\frac{1}{2}"\frac{\partial^{2}}{"\partial\theta^{2}}]P|l’v$

where$\mathrm{t}’=|\mathrm{v}|$(8).

. It should be noted that the singularity at$\mathrm{t}’=0$(8)still appears in the irst term of the right hand side of

Now we assume the expansion.

$P(10)$ $\theta)=\sum_{n=0}^{\infty}P_{n}(v)\cos n\theta=P_{0}(l’)+P_{1}(\mathrm{z}’)\cos\theta$(9)

Here therms with $n\geq 2$ are irrelevant, because the contribution of $n\geq 2$ is orthogonal to those of $n=0,1$ , andthe effect of gravity appears in the term of $n=1$ . Thus, the normalizationis reduced to$\mathrm{t}\mathrm{e}$$(^{??})$$2 \pi\int_{0}^{\propto}dl’vP_{0}(\mathrm{c}’)=1$

It is easy to verify$P_{0}$

and$P_{1}$,

and$\int_{0}^{2\pi}d\theta$

can be represented by a function$P_{0}(v)=f(\iota’)e^{-\eta v}$

,

cos2

$\theta\int_{0}^{\infty}d\tau$,$\iota\prime P_{1}(v)=0$

.

(10)

$\mathrm{f}(\mathrm{v})$$P_{1}( \iota’)=\frac{2D}{g}f’(v)e^{-\eta 2\prime}$

.

(11)

With the aid of eq.(ll) eq.(8) can be rewritten as
$f’+$

(, )
$\frac{1}{1}$$-\eta$$f’-( \frac{1}{l^{2}},+\epsilon\eta^{2})f’+\epsilon\eta^{3}f=0$

,

(12)

where $\epsilon=(g/\zeta)^{2}/2.[27]$Since we do not know the general procedure to obtain the solution of the third order ordinary differential equation(12) and it may not be easy to get the numerical solution for eq.(12) around the singular point $v=0$ , we rewrite it
as a set of the second order differential equations under the assumption of small . For this purpose, we expand$\epsilon$$f( \iota’)=\frac{\eta^{2}}{2\pi}[1+\epsilon f^{(1)}+\epsilon^{2}f^{(2)}+\cdots]$.

Skipping the detailed derivation we write the result of analytic calculation of(13)$f^{(1)}$

as[27]$f^{(1)}(v)= \mathrm{c}_{0}\{G_{23}^{31}(\eta v|_{0}^{0} 01 0) -\eta_{l\prime}(\gamma-1)+(\gamma+\eta\iota’)\ln\eta\iota’-\frac{1}{2}(\ln\eta v)^{2}\}$$+ \frac{\eta^{2}}{2}\iota\prime^{2}+c_{1}$

(14)

253
1

0
-1

-2

-3

$\mathrm{P}_{\mathrm{n}}(\mathrm{c})$

-4

-5
-6

-7
-8

0.5

0

1

2

1.5

2.5

3

4

3.5

$\mathrm{C}$FIG. 5: The scaled(solid line) and, where(dashed line) as functions of
the figure. For
we plotto remove the effect of the expansion parameter .
$P_{0}(c)$$P_{1}(c)$where$c_{0}$$c_{1}$

and$c_{1}$

is a constant and$P_{1}(c)$$c=\mathrm{v}/\mathrm{v}\mathrm{o}$$P_{n}(c)$denotes Po(c) or Pn (c) in$f^{(1)’}(c)e^{-\sqrt{6}}\mathrm{C}$$\epsilon$$G_{\sim}^{3}.$)$31(\eta\tau’|00$

$01$

$0)$

is Meijer’s$\mathrm{G}$function [29]. From the normalization condition constantsin eq.(14) become$c_{0}=1$,

and$c_{1}=4- \delta_{1}+\frac{\pi^{\sim})}{12}.-\frac{\gamma^{-}}{2’}.\simeq$-2.98961,(15)

respectively.
V.

DISCUSSION

Let us discuss our result. First, we can obtain higher order terms such as
and $f^{(3)}(v)$ , because the homogeneous solutions inandcan be used in ally order. However, these solutions may have complicated forms
and we had better use numerical integrations to represent inhomogeneous terms. Because of the limitation of thelength of this paper, we have omitted to give the explicit forms of higher order terms here.
Second, Kawarada and Hayakawa[26] assume that the diffusion coefficient in the real space is proportional to the
diffusion constant in the velocity space, but this assumption may not be true for our case.
We also note that the behavior of VDF under the influence of in granular particles in an inclined container
with the vibration does not coincide with the theoretical prediction presented here quantitatively in our preliminary
simulation. As shown in Fig.5,
andhave peaks around $c=1$ and become negative near $c=0$ . The behaviorof the simulation for granular particles is qualitatively similar, but the peak position is located near $c=0.5$ and thenegativity in the vicinity of $c=0$ is not large in the simulation. It is not surprised in the difference between ourtheoretical result and our preliminary simulation mentioned here, because the theoretical model is oversimplified withignoring vertical motion of particles and the collisions among mobile particles. In addition, our search for adequate
parameters may not be enough, and higher order terms $f^{(n)}(v)$ with $n\geq 2$ may be important in the actual situation.
In the simulation, we also find the existence of some peaks around the central peak atof VDF. These extra peaks
may me ffom the steady oscillation of particles among scatters. From the existence of these small peaks the peak
positions of VDF in Po (v) alld
may be entrained into small region. In any case, we will have to check whether
the connection between Langevin equation with Coulomb friction and the vibrated granular particles confined in aquasi-tw0-dimensional box is superficial. The detailed quantitative comparisons will be reported elsewhere.
Here, we comment on other possible situations to reproduce Langevin equation with Coulomb friction. One ofsimplest one is that the motion of frictional disks on an inclined slope with dense random scatters as in sectio$f^{\mathrm{t}- \mathrm{I}}.,(\mathrm{t}’)$$\mathrm{e}\mathrm{q}\mathrm{s}.(??)$$(^{??})$$g$$P_{0}$$P_{1}$$\mathrm{t}’=0$$\mathrm{c}\mathrm{o}$$P_{1}(\iota’)$

$c$

254
2. This syste can reproduce an exponential VDF with macroscopic condensed particles in tlle immobile (zerotemperature) state. However, VDF is highly anisotropic in this situation alld most of particles’ motions occur in the
direction parallel to the external ield, $i.e$ . along the slope. Thus, we cannot use the expansion in eq.(9). We also needto check whether we call use our theory ill the simpler situation of dense particles in a vibrating box. It is interestingthat we look for physical situations that we can use Langevin equation with Coulomb friction.
$\ln$$\mathrm{V}\mathrm{I}$.

CONCLUSION

In conclusion, we propose all experimental accessible system for granular gases. The VDF in a ’steady state’obeys an exponential-like function but changes Gaussian-like distribution function when free-cooling starts. This
exponential VDF is caused by Coulomb’s friction force. Thus, we propose the Langevin equation with Coulomb’s
friction to reproduce the results of our simulation, we have also developed the theory of Langevin equation with
Coulomb friction and obtained the steady solution of VDF under the influence of a steady external field.
This work is partially supported by the Grant-in-Aid for Scientific Research (Grand N0.15540393) of the Ministryof Education, Culture, Sports, Science and Technology, Japan.

[1] H. M. Jaeger, S. R. Nagel and R. P. Behringer, Rev. Mod. Phys. 68, 1259 (1996).
[2] P. Dantu, in Proc. of the 4th. Intemational Conf. on Soil Mechanics and Foundation Eng., $Vol.l$ (London, Butterwoths,

1957) p.144.
[3] F. Radjai, M. Jean, J. J. Moreau and S. Roux, Phys. Rev. Lett. 77, 274 (1996):[4] D. M.Mueth, H. M. Jaeger and S. R. Nagel, Phys. Rev. $57$ , 3164 (1998).
[5] T. P\"oshel and S. Luding edit, Granular Gases (Springer, Berlin, 2001)[6] S. E. Esipov and T. P\"oschel, J. Stat. Phys. 86, 1385 (1997)
Noije and M. H. Ernst, Granular Matter 1, 57 (1998).
[7] T. P. J.
, 2019 (1998)[8] J. J. Brey and D. Cubero, Phys. Rev.
, 010301 (R) (2003)[9] H. Nalmlislli, Phys. Rev.
[10] R. Kawahara arld H. Nakanishi, J. Phys. Soc. Jpn. 74, 68 (2004).
, 5583 (1995).
[11] S. Warr, J. M. Huntley and George T. H. Jacques, Phys. Rev.
, 7066 (2000). Hansen, Phys. Rev.
[12] S. Warr, J. M.
and[13] Y. Murayama a1ld M. Sano, J. Phys. Soc. Jpn. 67, 1826 (1998).
, 2468 (R) (1999).
[14] J. S. Olafsen and J. S. Urbach, Phys. Rev.
[15] W. Losert, D. G. W. Cooper, J. Delour, A. Kudrolli and J. P. Gollub, Chaos 9, 682 (1999).
[16] F. Rouyer and N. Menon, Phys. Rev. Lett. 85, 3676 (2000).
, 011301 (2004).
[17] S. J. Moon, J. B. Swift and H. L. Swinney, Phys. Rev.
, 050301 (R) (2001).
[18] D. L. Blair and A. Kudrolli, Phys. Rev.
[1 ] Y-h. Taguchi and H. Takayasu, Europhys. Lett. 74, 1268 (1995), 658 (1995)[20] K. Ichiki and H. Hayakawa, Phys. Rev.
, 031304 (2001).
[21] M. Isobe, Phys. ffiv.
[22] C. S. Campbell, Ann. Rev. Fluid Mech. 22, 57 (1990)[23] N. Mitarai, H. Hayakawa and H. Nakanishi, Phys. Rev. Lett. 88, 174301 (2002).
[25] T. Kawaguchi, T. Tanaka and Y. Tsuji, Powder Ibchol. 96, 129 (1998)[25] T. Kawaguchi, M. Sakamotom, T. Tanaka and Y. Tsuji, ibid. 109, 3 (2000).
[26] A. Kawarada and H. Hayakawa, J. Phys. Soc. Jpn. 73, 2037 (2004).
[27] H. Hayakawa, c0nd-mat/0407789 (to be published in Physica ).
[28] P. A. Cundall and O. D. Strack, Geotechnique 29, 47 (1979).
[29] C. S. Meijer, Proc. Nederl. Akad, Wetensch. 44, 1062 (1941)$\mathrm{E}$$\backslash ’\mathrm{a}\mathrm{n}$$\mathrm{E}57$$\mathrm{E}67$$\mathrm{E}52$$\mathrm{H}\backslash 1\mathrm{n}\mathrm{t}1\mathrm{e}_{\mathrm{u}}\mathrm{v}$$\mathrm{E}60$$\mathrm{J}.\mathrm{P}$$\mathrm{E}60$$\mathrm{E}69$$\mathrm{E}64$$\mathrm{E}52$$\mathrm{E}64$$\mathrm{D}$