See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/1833732

Molecular Dynamics of Comminution in Ball Mills

Article in The European Physical Journal B · June 2000

DOI: 10.1007/PL00011052 · Source: arXiv

CITATIONS READS

20 40

3 authors:

Volkhard Buchholtz Jan Alfred Freund

21 PUBLICATIONS 376 CITATIONS Carl von Ossietzky Universität Oldenburg
86 PUBLICATIONS 1,298 CITATIONS
SEE PROFILE
SEE PROFILE

Thorsten Pöschel
Friedrich-Alexander-University of Erlangen-Nürnberg
313 PUBLICATIONS 5,003 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Granular matter View project

Organized Conferences View project

All content following this page was uploaded by Thorsten Pöschel on 29 October 2014.

The user has requested enhancement of the downloaded file.

Eur. Phys. J. B 16, 169–182 (2000)
THE EUROPEAN
PHYSICAL JOURNAL B
EDP Sciences
c Società Italiana di Fisica
Springer-Verlag 2000

Molecular dynamics of comminution in ball mills

V. Buchholtz1 , J.A. Freund2 , and T. Pöschel2,3,a
1
Logos Verlag Berlin, Michaelkirchstraße 13, 10179 Berlin, Germany
2
Humboldt-Universität, Institut für Physik, Invalidenstraße 110, 10115 Berlin, Germany
3
ICA1, Universität Stuttgart, Pfaffenwaldring 27, 70569 Stuttgart, Germany

Received 19 January 2000

Abstract. We investigate autogenous fragmentation of dry granular material in rotating cylinders using
two-dimensional molecular dynamics. By evaluation of spatial force distributions achieved numerically for
various rotation velocities we argue that comminution occurs mainly due to the existence of force chains. A
statistical analysis of these force chains explains the spatial distribution of comminution efficiency in ball
mills as measured experimentally by Rothkegel [1] and Rolf [2]. For animated sequences of our simulations
see http://summa.physik.hu-berlin.de/∼kies/mill/bm.html

PACS. 81.05.Rm Porous materials; granular materials – 83.70.Fn Granular solids

1 Introduction chains. Our results will hint at how to improve the effi-
ciency of a milling machinery widely applied in industry.
In the past years molecular dynamics simulations of gran-
ular systems – such as sand, fertilizer, grain and others –
have been developed towards a reliable tool for an investi- 2 The simulation model
gation of granular systems. Unlike a few years ago, when
physicists started to investigate granular assemblies con- 2.1 Molecular dynamics
sisting of just a few hundred spherical granular particles
in two dimensions, e.g. [3–9], and less than 100 particles For the simulation we assumed pairwise interaction forces
in three dimensions, e.g. [10], today, owing to the rapid and applied the model by Cundall and Strack [4] and Haff
evolution of computer facilities, we are able to simulate and Werner [3]: A pair of two-dimensional spherical par-
complex systems of many thousands of particles in two ticles i and j interacts through the force
and three dimensions even accounting for non-spherical
Fij = FijN n + FijT t, (1)
grain features.
Due to this progress molecular dynamics of granu- if their distance is smaller than the sum of their radii
lar systems by now can be applied to the simulation of Ri + Rj − |ri − rj | = ξij > 0. The unit vectors in normal
technologically important processes striving for a bet- and tangential direction are defined as
ter understanding of relevant details of the process and,
henceforth, for an optimization of industrial technologies n = rij / |rij | (2)
(e.g. [11–17]). Molecular simulation techniques, therefore,  
0 −1
might develop to a prospective engineering tool partly re- t= · rij / |rij | , (3)
placing costly laboratory experiments by “computer ex- 1 0
periments”. with rij ≡ ri − rj . The forces in normal and tangential
In the present paper we apply the method of molecular direction are respectively given by
dynamics to the simulation of autogenous dry comminu-
tion in a tumbling ball mill [18]. To this end we devised FijN = Y ξij − meff
ij γN ṙij · n (4)
a novel algorithm which accounts for the fragmentation

of particles. First we will demonstrate that experimen- FijT = sign(vij
rel
) min meff rel N
ij γT vij , µ Fij (5)
tally known results, in particular those by Rothkegel and
with
Rolf [1,2], can be reproduced up to a good degree of ac-
curacy. We will show that the efficiency of comminution rel
vij = ṙij · t + Ri Ωi + Rj Ωj (6)
in ball mills is mainly determined by the presence of force mi mj
eff
mij = · (7)
a
e-mail: thorsten@physik.hu-berlin.de mi + mj
170 The European Physical Journal B

Here Ri , mi , ri , ṙi and Ωi are respectively the radius, the the material properties of the particles [24,27]. Besides
mass, the position, the velocity and the rotation velocity some conceptional difficulties which can be overcome by
of the ith particle. Y = 8×106 g s−2 is the elastic constant applying numerical tricks [28] these algorithms are only
and γN = 800 s−1 , γT = 3000 s−1 are the parameters of applicable in case the duration of collisions is negligible in
damping in normal and tangential direction. The Coulomb comparison with the mean free flight time, i.e. the time
friction constant was assumed to be µ = 0.5. These empir- a particle moves without interaction. In this limit colli-
ical values which we used throughout all simulations have sions can be assumed as instantaneous events. Moreover,
proven to render realistic behaviour for a typical granular this implies that three-particle interactions are very rare
material. events. However, in the system under study, namely the
The normal force (4) is composed of an elastic repulsive ball mill, this premise does not hold because the majority
part and a dissipative part which acts against the direc- of particles permanently is in close contact with neigh-
tion of motion. For a collision of two-dimensional spheres, bouring particles.
3/2
i.e. disks, the Hertz contact law FijN ∼ ξij [19,20] reduces
to [21]
! 2.2 Fragmentation probability of grains
2 4(Ri + Rj )Er
ξij ∼ Fij N
+ ln , (8)
3 FijN It can be observed in experiments that different particles
of identical size and material vary with respect to their
where Er is the reduced elastic module, i.e. a material
tensile strengths. That means when compressing particles
constant. Equation (8) provides a relation ξ = ξ F N , of a sample by applying a fixed force only a fraction of the
however, in the molecular dynamics simulations we need sample actually will break. For a sufficiently large sam-
the inverse F N = F N (ξ). To calculate the forces due to ple this fraction can be rephrased in terms of a fracture
equation (8) one would have to invert (8) numerically probability.
for each particle contact in each time step, i.e., to solve The reason for this different behaviour can be ex-
a transcendental equation. Alternatively one could tabu- plained by the existence of flaws which provide sites
late the function F N (ξ) but due to large force gradients, for stress concentration and the initiation of a crack
in particular in the instant of contact, we have to simu- which, subsequently, propagates through the material.
late with double precision accuracy. Therefore, the table These flaws are distributed throughout the whole volume,
would be extremely large. Apart from the small logarith- however, several investigations have revealed that surface
mic term −FijN ln FijN in equation (8) one finds that the flaws activated by high tensile stress play the dominant
force is proportional to the compression ξij . Therefore, in role for the initiation of a crack [29]. Hence, the fracture
our simulations we used the linear law (4) which is a good probability is intimately related to the statistical distribu-
approximation in the force interval of interest. tion of surface flaws and the resulting stress distribution
Equation (6) describes the relative velocity of the par- over the particle’s surface.
ticle surfaces at the point of contact which results from Starting from some plausible assumptions concerning
both the tangential part of the relative particle velocity the flaw distribution and the behaviour of the material the
and the velocity of particle spin. fracture probability P of a single particle can be derived
The Coulomb friction law is taken into account in employing statistical reasoning. A particle of radius R,
equation (5): If the tangential force of two colliding par- subjected to a force which stores the specific elastic energy
ticles exceeds µ times the normal force the particles slide Wm = W/m, will break with probability [30]
upon each other feeling constant friction. In this way the

Coulomb law formulates a maximum transferable shear P (R, Wm ) ∼ 1 − exp −cR2 Wm z
. (9)
force.
As numerical integration scheme we used the Gear Here, c and z are material constants. The original deriva-
predictor-corrector method of sixth order (e.g. [22]). tion (considering rods of a brittle material) goes back to
Let us remark that there exist various models with di- Weibull [31]. This simple law has proven in practice by
verse descriptions of the interaction forces, e.g. [23]. In fitting various experimental data. The constants c and
three-dimensional calculations it is √crucial that the dissi- z can be extracted from measurements when plotting
pative term in (4) is proportional to ξ ξ̇ [24] (in viscoelas- [ln ln(1 − P )−1 ] vs. [ln Wm ]. Typical values of z are found
tic approximation). In three dimensions the used term ∼ ξ˙ to lie in the range 1.5, ..., 2.5.
in equation (4) yields wrong results [23,25]. Notice that in agreement with experimental results the
Another class of molecular dynamics are event-driven term R2 in the exponent (for fixed Wm ) in general predicts
simulations, e.g. [26] which require much less numerical ef- an anticorrelation between particle size and resistance to
fort than the method described above (force-driven molec- breakage. Clearly, this fact is rooted in the diminished
ular dynamics). In event-driven simulations one does not probability to find a sufficiently large flaw on a smaller
integrate the equations of motion explicitly but rather ex- surface.
presses the loss of mechanical energy in normal and tan- For spherical particles (or, more generally, for particles
gential direction due to a collision via coefficients of resti- having convex surface), assuming linear elastic behaviour,
tution N and T both of which can be determined from Hertz-theory [19,20] can be employed to derive a relation
 V. Buchholtz et al.: Molecular dynamics of comminution in ball mills 171

unique thickness L. The fracture probability correspond-

ing to equation (11) can be obtained in analogy to (9)
1 yielding
Fn
Fn
3
P (R, Wm ) ∼ 1 − exp (−cRLWm
z
). (12)
Note that the surface of a sphere ∝ R2 has to be replaced
4 by the surface of a cylinder ∝ RL (excluding top and
2
Fn bottom wall). Correspondingly equation (10) turns into
Fn
ξ = Y F N and, hence, the stored elastic energy for spheres
has to be replaced by
1 N (F N )2
Fig. 1. A typical situation of one particle being compressed W = ξF = ∼ (F N )2 . (13)
2 2Y
by four impacting particles. The forces are considered as inde-
pendent stress sources. Therefore, in simulations we evaluated Taking into account that the mass of a cylinder is propor-
only the maximum of all exerted forces to compute the fracture tional ∝ R2 , one obtains
probability P (R, Wm ).   N 2 z 
(F )
P (R, F ) ∼ 1 − exp −c̃R
N
. (14)
R2
between the elastic energy W stored in a sphere and the In Section 3 we will motivate the choice z = 2. With this
exerted repulsive contact force F specification one finally arrives at the following expression
for the fracture probability
  13 
2 2 D2 F 5
4 
W = ξF = , (10) P (R, F N ) ∼ 1 − exp −c̃R−3 F N . (15)
5 5 R

where ξ quantifies the particle’s deformation and D is an-

other material constant (expressed by the Poisson ratio 3 Fragment size distribution
and constant of elasticity) [19,20]. It should be kept in
Once a flaw has been activated it forms an initial crack
mind that formula (10) remains valid in the linear elastic
which rapidly propagates through the particle. Typical
regime only.
crack velocities in the range of 1500 m/s have been mea-
Applied to our simulation of ball mills we make use
sured [29]. The formation of a variety of differently sized
of these laws in the following way: Consider the situation
fragments occurs by virtue of branching cascades. Branch-
sketched in Figure 1: A particle placed in the centre is
ing is governed by a balance between the energy release
compressed by four impacting particles.
rate G and the so-called crack resistivity B [29] the lat-
The highest compressive stress occurs around the con- ter accounting for the creation of new surfaces. Because
tact point; a fact well-known to engineers [32] and also of the extremely short duration of the rupture process
found in numerical simulations [33,34]. Hence, the fracture one can safely neglect external energetic contributions to
probability of the centre particle is linked to the probabil- this energy balance. The energy release rate G depends
ity of finding a flaw of sufficient size in the very vicinity on the crack velocity and the crack length whereas the
of any of the four contact regions. Thus, the action of dif- crack resistivity B depends on the crack velocity and on
ferent impacting particles can be regarded as independent the number of branches. Due to a maximal crack velocity,
stress sources. That means we consider the elastic energy viz. the speed of sound, the balance between both terms
stored by each impacting particle separately, calculate the requires the formation of new cracks. From these consid-
related fracture probability according to formulae (9, 10) erations it can be shown [29] that the fragment size distri-
as a function of the particle’s size and the maximum of all bution (cumulative mass distribution) Q is a function of
acting normal forces F N the product Rf Wm , where Rf denotes the fragment size,

5z  i.e. Q = Q(Rf Wm ).
P (R, F N ) ∼ 1 − exp −c̃R(2−z 3 ) F N 3 ,
10
(11) Experimental evidence [30] for a scaling law Q ∼ Rfβ is
in agreement with the well-known empirical Schuhmann
and actually decide the question of breakage by subse- law [35]
quently drawing a random number for each of the contact  β
sites. Rf
Q∼ , (16)
The constant c̃ is related to the empirical constant c k
in equation (9); the precise connection is elaborated by
which itself can be regarded as an approximation of the
inserting equation (10) into (9) and expressing m by the
Rosin-Rammler law [36]
particle radius m ∝ R3 . "   #
Up to this point we discussed the fracture probabil- Rf
β
ity for three-dimensional spheres. In our simulations we Q ∼ 1 − exp − . (17)
k
modelled two-dimensional spheres, i.e. circular discs of
172 The European Physical Journal B

The variable k has the dimension of a length and, thus,

can only be identified with the size of the original particle. 11111111
00000000
00000000
11111111 00000
11111
0000000
1111111
111111
000000 00000
11111
A rigorous derivation of equation (17) – compatible 00000000
11111111
00000000
11111111
00000000
11111111
00000
11111
000000
111111
00000
11111
000000
111111
00000
11111
1111111
0000000
00000
11111
0000000
1111111
00000
11111
00000000
11111111 000000
111111
00000
11111 0000000
1111111
00000
11111
0000000
1111111
with an exponent β = 1 – starting from rather mild as- 00000000
11111111 000000
111111
00000
11111
000000
111111
00000
11111
0000000
1111111
00000
11111
00000000
11111111 00000
11111
000000
111111
00000
11111 0000000
1111111
00000
11111
00000000
11111111
00000000
11111111 000000
111111
00000
11111 0000000
1111111
00000
11111
sumptions and applying Poisson statistics was performed 00000000
11111111 000000
111111 0000000
1111111
0000000
1111111
00000000
11111111 000000
111111 0000000
1111111
00000000
11111111
by Gilvary [37] (see also [38]).
It has to be mentioned that the derived laws
(16, 17) describe the distribution of fragments only Fig. 2. Sketch of the particle fragmentation procedure. Im-
below a certain size (endoclastic vs. exoclastic distri- mediately after the fragmentation the particles must overlap
bution [37]). For the largest fragments an equivalently due to mismatch of geometry of the spherical fragments (mid-
well-defined quantitative statement does not exist. The dle), a small repulsive force helps to overcome this unphysical
deviation between theoretical predictions and experimen- transient situation.
tal observations typically occurs for grains which together
collect about 75–90% of the total mass [29]. Moreover,
some diverging theoretical mean values can be understood in general, it will be impossible to place the two frag-
from the fact that the assumption of equidistributed flaws ments avoiding an overlap. However, since a considerable
on all length scales neglects the effect of flaw depletion ap- overlap corresponds to a situation of extreme compres-
plying to tiny fragments. Flaw depletion thus sets a lower sion an unrealistically strong repulsive normal force (4)
bound to the range of applicability. Insofar the size of frag- together with an exceptional deformation energy would
ments which is satisfactorily described by the distributions be the consequence. We bypassed this problem by mod-
(16, 17) is restricted to the intermediate range. ifying the interaction between the two fragments in the
When fitting both the more refined Rosin-Rammler following way: Up to the moment of first complete sepa-
law (17) and the simplified Schuhmann law (16) to diverse ration of two such fragments i and j the normal force (4)
experimental data it was found that in some situations at the (k + 1)st time step is replaced by
the Rosin-Rammler law proved superior [29] (mill prod-
ucts [38]) whereas in other cases (specimen fractured in  ∗
 Y ξ (k) − meff
ij γN (ṙi − ṙj ) · n
gelatin [38]) the Schuhmann law was even more adequate.  ij
Finally, we mention the population balance model [39] FijN = if (ṙi − ṙj ) · n < 0 (19)


which assumes that the fragment size distribution can be  otherwise
normalized to the initial particle size, i.e. Q = Q(Rf /R).
This assumption is based on empirical observations. More- with
over, its central assertion can also be derived when start-
ing from a Weibull exponent z = 2 and applying statistical
∗
reasoning similar to the one sketched above. ξij (k) = |ri (k − 1) − rj (k − 1)| − |ri (k) − rj (k)| . (20)
In our simulations we employed the Rosin-Rammler
law (17) with k = R and an exponent β = 1. Since ex- This means that the two fragments, on the one hand, re-
perimental data have shown the mean dimension of the sist to further compression like standard particles and, on
largest fragment to be roughly 75% of the original parti- the other hand, experience a small repelling force  which
cle size [38], i.e. Rfmax = 34 R, the normalized fragment size drives them apart. The repelling force is chosen sufficiently
distribution simply reads small in order to suppress unphysical energy gain of the
  system; this means that the unavoidable energy, which due
R to this force is pumped into the system anyhow, must be
  1 − exp − f
Rf R negligible in comparison with the mean kinetic energy of
Q =   · (18)
R 3 the grains. After the moment of first complete separation
1 − exp − of the two fragments they also interact in the standard
4
way, as dictated by equation (1). Of course, the inter-
action of two such fragments with all other particles is
3.1 Molecular dynamics modelling of fragmentation never modified. By trimming of  we could not only sup-
press the elastic energy gain but also evidenced that the
In the present study we exclusively utilize spherical par- period requiring a modified interaction force was compara-
ticles (disks), which means that also the fragments of tively short. The procedure of fragmentation is sketched in
a disrupted particle must be spherical. The single frac- Figure 2. We want to mention that multiple fragmenta-
ture statistics, as discussed in Section 3, is put into prac- tion, i.e. further fragmentation of the fragments is in-
tice by dismembering a cracking disk into two fragments. cluded in the model and is frequently observed in sim-
Whereas the size of one fragment is chosen in accordance ulations.
with the Rosin-Rammler law the size of the second one is Let us remark that our algorithm was flexible enough
determined by mass conservation, i.e. for two-dimensional to include multiple fracture events: By this we mean the
particles by area conservation. observation that a fragment was apt to break itself even
In the majority of situations a breaking particle before its first complete separation from its original twin
is closely surrounded by neighbouring ones. Hence, fragment.
 V. Buchholtz et al.: Molecular dynamics of comminution in ball mills 173

An animated sequence of fragmentation can be found

in the internet under
http://summa.physik.hu-berlin.de/∼kies/mill/bm.html
Finally, we mention that the proposed algorithm is not
the only way to model particle fragmentation. An alter-
native algorithm has been put forward by Åström and
Herrmann [40].

3.2 Specification of the model system

As depicted in Figure 3 the cylindrical wall of our system

has been modelled by an ensemble of spheres possessing
the same material constants as the grains inside the con-
tainer. However, the motion of wall spheres is not affected Fig. 3. The cylindrical container is modeled by particles
by impacts but, instead, strictly governed by the continu- placed along the circumference of the container. Roughness is
ous rotation of the cylinder. Moreover, in spite of identical taken into account by stochastic variation of the sizes of these
elastic and dissipative constants the wall particles were particles.
considered to be much harder so they never would break.
A suitable size distribution of wall particles generates the accumulated dust mass exceeded an upper threshold,
the effect of surface roughness. This type of modelling of course, obeying mass conservation. In this way we con-
has been used in several numerical simulations and is also served the average filling level throughout the simulation
used in experiments to guarantee well-defined boundary and, thus, could reach a steady state of the continuous
properties. The wall particles have been placed along the operation mode.
circumference of the cylinder. Just like in realistic ball
mills we evenly distributed 16 toolbars along the inner
periphery of the cylinder which serve to lift particles and 4 Simulation results
thus prevent them from sliding downhill. We composed
the toolbars again from a collection of rigid spheres. As 4.1 Mass distributions
reported in literature [32] the maximum milling efficiency
is achieved for a filling degree, i.e. the ratio of material As a first result we studied the frequency dependent frag-
and container volume, of ca. 40%. ment size distribution. This distribution is not necessarily
The simulation results to be presented in the next sec- equal to the distribution of Section 3, which was the result
tion have been done with a system consisting on average of single fracture experiments.
of about 1000 particles. Due to fragmentation the num- It is known [41–43] that for an arbitrary fracture mech-
ber of particles was not constant, i.e. the particle number anism, which hierarchically is applied on all length scales,
fluctuates. The container was modeled by about 500 wall one will always generate a log-normal distribution of frag-
particles as described above. ment sizes; this assertion is a pure consequence of the cen-
Initially the spheres which model the grist were ran- tral limit theorem. However, in our modelling the ball mill
domly positioned inside the container avoiding any over- is operated as an open system, i.e. there is a particle flow
lap. This artificial initial condition had to be relaxed by into and out of the mill. Therefore, analyzing the asymp-
evolution of the system. In this way we achieved a realistic totic size distribution is also of theoretical interest.
configuration for any chosen rotation velocity. Obviously Furthermore, the explicit knowledge of the distribution
the typical asymptotic configuration depends on the ro- plays also a practical role, because from that knowledge
tation frequency ν, therefore, we started our numerical one can decide whether a multi-level fragmentation pro-
monitoring only after transients had died out. In Figure 4 cess is to be preferred to a single-level one or vice versa.
we present snapshots of typical configurations depending A multi-level process can be either a distribution of the
on four different rotation frequencies. comminution process to several ball mills, where each is
In technologically relevant ball mills there is a con- filled with granular material of a specified narrow size dis-
tinuous transport of raw material into the mill and of tribution, or a spatial separation of material of different
milled material out of the device. In most industrial ap- size in one mill.
plications this exchange of raw material and final product Figure 5 shows the differential fragment size distribu-
is accomplished by axial transport (z-direction) of mate- tion, i.e. the mass fraction of all particles with size within
rial. Of course, with our two-dimensional model we cannot a certain interval, for four different rotational frequencies.
account for axial transport, therefore, we simulated the The black circles in each column correspond to measure-
material exchange in continuous operation mode in the ments at different times and thus reflect the temporal vari-
following way: Whenever fragments were generated whose ance. The time period of each simulation was 10 seconds
size dropped below a minimum we removed them from real time (beyond transients).
the container. After exact bookkeeping of the mass of re- Obviously the lowest curve (2.5 Hz) changes only very
moved “dust” particles we reinserted a big particle when weakly during the simulation, which can be interpreted
174 The European Physical Journal B

Fig. 4. Typical dynamical configurations depend on the rotation speed: 0.5 Hz (top left), 1.0 Hz (top right), 1.25 Hz (bottom
left), 2.5 Hz (bottom right). Grey scale codes the instantaneous maximum compression (normal force) experienced by each
particle through contact with its neighbours (light: small normal force, dark: large force). Obviously the largest forces occur in
the range of intermediate rotation velocities (1.0 Hz, 1.25 Hz). Moreover, we identify the location of largest forces which can
initialize fragmentation to be close to the bottom of the rotating cylinder. This important observation will be discussed in detail
in Section 5.

0.50 Hz as a low rate of fragmentation. This curve – and with

differential fragment size distribution

restrictions also the top curve (0.5 Hz) – is determined

by the initialization. Unfortunately, the middle curves are
1.00 Hz not significantly different, which mainly is a consequence
of the short simulation time (10 seconds). To answer
the interesting questions (theoretical distribution, multi-
1.25 Hz level process) it will be necessary to simulate considerably
longer time periods.

2.50 Hz
4.2 Optimization of the efficiency
msieve mass intervals mraw
In the context of industrial applications the notion of ef-
ficiency mainly focuses on the aspects of maximum com-
Fig. 5. The differential fragment size distribution for different
minution rate and, perhaps, of power consumption. There
rotation velocities. The black circles represent the values for a
are many free parameters which can be varied to maximize
fixed mass interval at different times. The variance of the circles
the efficiency: speed of rotation, filling degree, size distri-
illustrates fluctuations of the distribution. The solid lines are
tracing the averages of each interval. The left border of the
bution, cylinder diameter, grist size distribution, various
scale is fixed by the size of the sieve (msieve ), the right border mill types, etc.
by the size of the refilled big grains (mraw ). In our analysis we investigated the dependence of
the comminution rate on the speed of rotation. In order
 V. Buchholtz et al.: Molecular dynamics of comminution in ball mills 175

to compare our numerical results with experimental data

we normalized quantities in the following way: The 1.0
speed of rotation is normalized to the velocity nc =
p

rel. commin. rate ρ max

g/M(2π)−1 for which the centrifugal force balances the 0.8
gravitational force whereas the comminution rate is nor-
malized to its measured maximum. The result of our simu-
0.6
lation is plotted in Figure 6 using filled circles and dashed
line. Clearly, a maximum occurs around n ≈ nc whereas
the rate rapidly diminishes for larger or smaller velocities. 0.4
The full line in Figure 6 depicts a similar curve follow-
ing from experiments. The curve shows the relative power 0.2
consumption as a function of the normalized speed of ro-
tation. Qualitatively both curves are in good agreement;
systematic deviations can be explained by the following 0.0 0.5 1.0 1.5 2.0
observation: Strictly speaking both curves relate to dif- rel. speed of rotation n/nc
ferent physical quantities. The fraction rate only accounts
for the energy dissipated through fractures while the ex- Fig. 6. The dependence of the fraction rate on the rescaled
perimental curve includes all dissipative processes, for in- rotation velocity, simulation (filled circles plus dashed line) and
stance friction between grains, friction between the grains experiment (full line) taken from [32] (original in [44]). The
and the wall, dissipative impacts, heating, etc.). As a con- error bars for the simulation visualize the fluctuations during
sequence, the curve for the simulation should be found the simulation period of approximately 10 seconds real time.
below the experimental curve. Moreover, from the same
reasoning it is also clear that observed deviations are mi-
nor if the main portion of the overall dissipated energy As can be concluded from Figures 7 and 8, the major-
results from the fractures, i.e. around the maximum. ity of fragmentation events does not occur near the free
surface of the granular material but deep inside the ma-
4.3 Preferred fragmentation locations terial near to the walls. This strange behaviour was ob-
served already in experiments by Rothkegel and Rolf [1,2]
The series of pictures in Figure 7 illustrates the spatial dis- in two-dimensional model mills. They used instrumented
tribution of fragmentation locations for four different rota- balls which emitted a flash of light whenever the strain ex-
tion velocities. We restricted our presentation to the most ceeded a fixed threshold. The positions of recorded flashes
important lower right segment of the ball mill (clockwise were digitized and from the statistics of these data they
rotation). To compute this distribution the coordinates of inferred the stationary strain distribution inside the ball
breaking particles were stored during the simulation. Af- mill. For various threshold values they always found the
ter the simulation the fraction of all fragmentation events peak strain deep inside the material (see Fig. 8). This was
occurring at places inside little cells (coarse graining) was surprising because inside the material relative grain veloc-
computed. In Figure 7 grey values code the frequency: ities are rather low in comparison with relative velocities
dark means high frequency and light low frequency at the surface.
respectively. This experimental finding, together with its validation
In the bottom right figure one finds a few fragmen- by our simulation (see the previous section), calls for an
tation events near the cylinder wall. The restriction to explanation. We will solve this challenging, so far open
medium-grey values results from the fact, that the rare question using the powerful possibilities of molecular dy-
events are equidistributed across the inner perimeter. namics, namely the option to monitor physical quantities
The predominance of rather light grey (or even white) which are hardly accessible through experiments.
pixels in the rest of the figures hints at regions where In the following simulations we ignored fragmentation
fragmentation events are rare as compared to the sparse events – which are not relevant for our explanation – sim-
medium-grey to dark pixels which indicate regions of high ply for the sake of computational efficiency, so we could
fracture probability. The observation that most fragmen- afford longer simulation time. This time the radius of the
tation events occur deep inside the material and not near two-dimensional ball mill was set to M = 4 cm. The
the surface, which is heavily agitated by impacts, is, at mill was filled with N = 800 spherical grains with radii
first sight, somewhat counter-intuitive. An explanation for equally distributed in the interval [0.05, 0.11] cm. We per-
this fact will follow from a closer analysis of the local formed simulations with three different rotation velocities
pressure concentration in the next section. ΩI = 2 Hz, ΩII = 10 Hz and ΩIII = 19 Hz.
Figure 9 shows snapshots of our performed simula-
5 Spatial distribution of pressure in ball mills tions. The pictures in the top row are taken form a sim-
ulation with rotation velocity Ω = ΩI . The velocity is
In this chapter we will dwell on the physical origin sufficient to observe a continuous flow at the surface. The
of the unusual pressure distribution [45] measured in surface shape is similar to a plane. The regime with dis-
experiments. continuous flow (for lower rotation velocities) will not
176 The European Physical Journal B

Fig. 7. Spatial distribution of fragmentation locations for different rotation velocities: 0.5 Hz, 1.0 Hz, 1.25 Hz, 2.5 Hz (from top
left to bottom right). Grey values code the frequency of fragmentation in a certain region (cell), dark meaning high frequency
and light low frequency respectively. Preferred fragmentation regions are found near the bottom of the mill, i.e. deep inside the
granular material. This observation is in good agreement with experimental results and the results from Section 5.

Fig. 8. Spatial distribution of the strain for different rotation velocities and filling grades. Again, grey levels code strain values.
The regions of maximal strain are located inside the material. This figure is redrawn from reference [1].

be investigated in this paper, because it is not relevant in these pictures code the local pressure Pi acting on grain
for technological applications. For a detailed discussion i given by
of the interesting observations at the transition between
X
continuous and stick-slip flow we refer to the work by Pi = FijN . (21)
Rajchenbach [46]. j
The picture in the middle row of Figure 9 depicts snap-
shots of a simulation with Ω = ΩII . The surface of the The index j is running over all particles which are in con-
material no longer shapes a plain but the grains fall down tact with grain i. Dark circles denote high pressure, light
a steep slope hitting a flat surface. spheres low pressure.
The bottom row of Figure 9 exhibits the snapshot of Obviously many of the heavier stressed particles are
a simulation with the highest rotation velocity Ω = ΩIII . deep inside the material, which is in good agreement
Here the grains are carried away with the fast rotation with the experimental observation by Rothkegel and
of the cylinder. Beyond a certain point they are likely Rolf [1,2]. Most importantly, the pictures reveal that
to loose contact with the wall and follow a free parabola grains experiencing peak stress often form linear struc-
before impacting back down on the surface. tures. In the following we focus on the development and
Moreover, the snapshots of Figure 9 provide an intu- the properties of such force chains. By numerical analy-
itive understanding, what mechanism is responsible for the sis we will prove that those force chains are the crucial
pressure distribution mentioned above. Again, grey values physical reason for the observed pressure distribution.
 V. Buchholtz et al.: Molecular dynamics of comminution in ball mills 177

Fig. 9. Snapshots of the simulation for different rotation velocities. Grey values code the local pressure Pi (see text). Linear
segments connect grains belonging to a force chain. ΩI = 2 Hz (top), ΩII = 10 Hz (middle), ΩIII = 19 Hz (bottom).

We defined force chains by a set of three self-evident larger as compared to the average pressure of the neigh-
conditions: Grains i, j and k are considered to be members bouring grains not belonging to a force chain, i.e. the force
of the same force chain if: distribution is strongly inhomogeneous.
1) particles i, j and j, k are next neighbours, The development of such force chains is not restricted
to grains in a ball mill but has been observed in differ-
2) the pressure acting on each of the grains exceeds a
ent granular systems by various experimentalists [47–50]
certain threshold (Pi > 1000 g cm s−2 ),
and numerically [51–53]. It could be shown, that one of
3) the connecting lines between i, j and j, k form an angle the granular phases is characterized by the observation
larger than 150◦ , i.e. the centres of three grains almost of force chains [52]. From this one might conclude that
fall on a line. the occurrence of force chains is an inherent feature of all
These three conditions are evaluated by a computer algo- granular matter, which may substantially contribute to its
rithm. In Figure 9 all grains belonging to a force chain are specific characteristics.
graphically connected by lines. Figure 10 comprises six plots which come in three
We see that most of the harder stressed grains could pairs; the lower plot of each pair always depicts the fre-
be assigned to a force chain. From this we conclude that quency distribution of force chains whereas the related
the main part of the static and dynamic pressure prop- upper one presents the (statistically evaluated) maximal
agates along the force chains. The pressure acting on a pressure both as a function of the length L of a force chain
highly stressed grain inside a force chain is up to 100 times (which is nothing but the number of particles belonging
178 The European Physical Journal B

6000 in a force chain (see Sect. 2.2). For lower rotation velocities
ΩI and ΩII a monotonous increase of the maximal pres-
5000
sure with the length L can be observed. For a discussion
4000
of these curves one has to take into account that the pairs
of plots in Figure 10 do not constitute independent quan-
3000 tities. Choosing a sample of L particles from an ensemble
I 10000 of particles on which forces with a given distribution are
1000 acting, the maximal force in this sample will increase with
increasing L. Therefore, in order to prove significance, the
100
result in Figure 10 has to be weighted against this purely
10 statistical effect.
1 The probability to measure a maximal pressure P max
5 10 15 20 in the interval P max ∈ [x; x + dx] in a sample of size L is
6000
equal to the probability that all particles feel a pressure
P < x + dx minus the probability that all L particles feel
5000
a pressure P < x.
4000
pL (P max ∈ [x; x+dx]) = [p(P < x+dx)]L −[p(P < x)]L .
3000 (22)
II 10000
1000 Assuming equally distributed pressures P ∈ [0 : Pm ],
p(P < x) means
100
10 Zx
1 x
1 p(P < x) = dx0 = · (23)
5 10 15 20 Pm Pm
6000 0

5000 By inserting equation (23) into equation (22), one finds in

first order of dx
4000
1
pL (P max ∈ [x; x + dx]) = L
L xL−1 dx. (24)
III 3000 (Pm )
10000
1000 Therefore, the average maximal pressure as a function of
sample size L is given by
100
10 ZPm
1
1 hP max
(L)i = L
x L xL−1 dx
5 10 15 20 (Pm )
0
Fig. 10. Within each pair of the triplet (I, II, III) the lower 1L L+1
plot presents the frequency distribution whereas the upper one = L
Pm
shows the maximal pressure inside a force chain both as a
(Pm ) L + 1
function of its length L (abscissa for all plots). ΩI = 2 Hz (top), L
= Pm . (25)
ΩII = 10 Hz (center), ΩIII = 19 Hz (bottom) (explanation of L+1
lines see text).
The dashed lines in Figure 10 show this statistical effect
for P max with Pm = 4800 g cm s−2 . Obviously, these
curves do not reach up to our numerical data of the top
to the chain). These pairs are shown for three different ro- and center figure; especially for larger L the dashed line
tation velocities varying from ΩI (top) to ΩIII (bottom). falls far below. The numerical data points in the bottom
The frequency of a force chain drops exponentially figure (high rotation velocity) to some extent coincide with
with increasing length L. For the highest rotation velocity the dashed curve, which only represents the statistical ef-
(ΩIII ) the statistics becomes unreliable for chains longer fect discussed above. Hence, we can conclude that for lower
than L = 13. This behaviour, which is in contrast to obser- velocities the maximum stress significantly concentrates
vations for slower measured rotations, is due to the higher within long force chains whereas for high rotation veloci-
energy input, which causes a decompactification of the ties a significant effect of force chains is not evident.
material. From experimental and numerical results (for instance
Concerning the earlier mentioned spatial strain distri- [49,54,55]) it is known, that the observation of force chains
bution the maximal pressure measured in each force chain is connected to a strongly inhomogeneous pressure dis-
is of paramount interest rather than the average pressure tribution. The distribution is essentially an exponentially
 V. Buchholtz et al.: Molecular dynamics of comminution in ball mills 179

decreasing function with increasing pressure. A simplified 1300

expression for this distribution is given by
1200
p(P ) = a exp{−a P }. (26)

[g cm s ]
−2
In analogy to equation (23) the probability to measure a 1100
pressure P < x is given by
Zx 1000

max
0 0
p(P < x) = a exp{−a x }dx = 1 − exp{−a x}. (27)

F
0 900

Inserting this into equation (22) we obtain in first order 800

of dx 1 3 5 7
nchain
 L
L
pL (P max ∈ [x; x+dx]) = 1 − e−a(x+dx) − 1−e−a x
Fig. 11. The maximal force acting on a particle as a function
= L(1−e−a x )L−1 e−a x a dx. (28) of the relative position in the force chain. The monotonous
behaviour indicates that forces are accumulating downwards
The average of the maximal pressure can be calculated by along the chain.
integrating

Z∞

L−1 Further insight into the mechanism of force chains is
hP max
(L)i = x L 1 − e−a x e−a x a dx = achieved through Figure 11 which shows the average of the
0 maximal force acting on a grain as a function of its posi-
  tion within the force chain. The average was performed for
 ZN 

L N
L all force chains of length L = 8. The top particle of a chain,
lim x 1 − e−a x − 1 − e−a x dx . (29)
N →∞ 0  corresponding to the maximal y-coordinate, was defined
0 as position 1. The monotonous increase of this maximal
force indicates that, on average, strain accumulates down-
An iterative splitting
wards along the force chain. In this respect a force chain
(1 − e−a x )i = (1 − e−a x )i−1 (1 − e−a x ) (30) is similar to a beam in a framework: The masses and mo-
menta of all particles above a force chain are supported
results in a summation of integrals of the type by this chain, with the consequence that lower particles of
the chain are heavier loaded than upper ones. This trend
ZN is clearly visible in the data of Figure 11. Furthermore,
1 a force chain can shield neighbouring particles, i.e. these
(1 − e−a x )i−1 e−a x dx = · (31)
ai grains experience a much weaker force which can be ob-
0
served in the simulation data as well.
Thus, we obtain Figure 12 shows the spatial pressure distribution for
different rotation velocities. Here grey values code the ab-
XL
1
solute local pressure, with light pixels indicating low and
hP max (L)i = · (32) dark pixels high values, in order to allow for a direct com-
i=1
ai parison of all three images.
The presented data are temporal averages. To extract
A restriction of the statistics to values of the pressure the spatial distribution of this data, the two-dimensional
P > 1000 g cm s−2 , performing an analogous calculation, space was coarse-grained. For each time step the positions
results in of all particles were mapped to the grid and the result-
ing instantaneous local pressure was added to its related
XL
1
hP max (L)i = + 1000 g cm s−2 . (33) cell. Notice that the averaged data does not reflect the
i=1
a i occurrence of force chains.
Surprisingly, for lower rotation velocities one finds that
From the simulations we extracted the pressure distribu- direct impacts with high relative velocities at the material
tion yielding – as expected – an exponential decrease. The surface do not result in high local pressure. Therefore,
prefactor is a = 0.0009 for Ω = 10 Hz and a = 0.001 these impacts are not relevant for the fragmentation of
for Ω = 2 Hz. The corresponding curves are plotted in grains in ball mills. Regions of high pressure can be found
Figure 10 as dotted lines and reasonably describe the mea- mainly near to the wall, deep inside the material. This
sured values. Thus, we indirectly can conclude from the observation is in good agreement with the experimental
measurement that we have an exponential distribution, result of Rothkegel and Rolf [1,2]. The extreme right figure
which indicates the occurrence of force chains. in Figure 12 reveals a different behaviour. The maximum
180 The European Physical Journal B

Fig. 12. Spatial distribution of the pressure. Grey values code the absolute pressure, with light pixels indicating low and dark
pixels high pressure values. The grey scale is unique for all three figures to allow for a direct comparison. From left to right:
Ω = 2 Hz, Ω = 10 Hz, Ω = 19 Hz. While for lower rotation velocities the maximum is deep inside the material near to the wall,
for higher rotation velocity it is near to the free surface, which is heavily agitated by impacts.

Fig. 13. Spatial frequency distribution of particles, exerted to a pressure P > 3000 g cm s−2 . Grey values, unique throughout
the images, code the frequency: dark pixels high frequency and light pixels low frequencies. The top row only accounts for
particles which do not belong to a force chain while, in contrast, the bottom row includes only particles being member of a force
chain. From left to right: Ω = 2 Hz, Ω = 10 Hz, Ω = 19 Hz. By visual inspection it becomes obvious that most of the heavier
loaded particles belong to force chains and are found near to the wall of the mill.

is near to the material surface where heavy impacts occur. belonging to a force chain contribute only weakly to the
The absolute values of the pressure are much lower than comminution. For the higher velocity ΩIII direct impacts
in the left figures. This observation relates to the existence of particles at the surface are the dominant part. Never-
of rather weak force chains. theless, absolute values are much lower than for ΩI and
ΩII . Again, this observation is in good agreement with the
Figure 13 shows the spatial distribution of the num-
experiment [1].
ber of particles, which experience a pressure Pi >
3000 g cm s−2 . In these figures the particles which be-
long to a force chain are only considered in the bottom
row whereas particles not belonging to a force chain are
contained in the top row. Figures 13 can be directly com-
6 Conclusion
pared with the figures of Rothkegel [1] (see Fig. 8); both
are found to be in good agreement. The maxima of the
We have simulated the process of autogenous dry com-
distribution for ΩI and ΩII lie near to the wall of the mill,
minution in a ball mill by the method of molecular dy-
inside the material as observed in the experiment. Those
namics. To this end we developed a novel algorithm which
maxima are an exclusive result of force chains. Grains not
 V. Buchholtz et al.: Molecular dynamics of comminution in ball mills 181

accounts for fragmentation of particles. Throughout our References

simulations all particles, including the fragments, were
modelled as spheres. To accomplish this idealization we
introduced rather short-lived unphysical transient situa- 1. B. Rothkegel, Örtliche Verteilung der Stoßenergien ≥
tions in which the fragments resulting from a cracked par- 22 mJ und die zugeordneten Bewegungszustände von Mod-
ellmahlkörpern in einer Modellkugelmühle, Ph.D. thesis,
ticle penetrate each other. The results of our simulations
Technical University Berlin, 1992.
are available in the internet as animated sequences under
2. L. Rolf (personal communication).
http://summa.physik.hu-berlin.de/∼kies/mill/bm.html
3. P.K. Haff, B.T. Werner, Powder Technol. 48, 239 (1986).
Employing our molecular dynamics algorithm we could
4. P.A. Cundall, O.D.L. Strack, Géotechnique 29, 47 (1979).
satisfactorily reproduce experimental phenomena, e.g. the 5. J.A.C. Gallas, H.J. Herrmann, S. Sokolowski, Phys. Rev.
normalized comminution rate as a function of the normal- Lett. 69, 1371 (1992).
ized rotation speed and the spatial pressure distribution. 6. Y.-h. Taguchi, Phys. Rev. Lett. 69, 1367 (1992).
From the achieved numerical data we were able to ex- 7. O.R. Walton, in Numerical Methods in Geomechanics,
plain an experimental result which was poorly understood edited by Z. Eisenstein (Balkema, Rotterdam, 1982).
so far [1]: Comminution processes in ball mills occur deep 8. O.R. Walton, R.G. Braun, in Joint DOE/NSF Workshop
inside the material rather than close to the surface where on flow of particulates and fluids (Ithaka, 1993), pp. 1-17.
the relative collision velocity is maximal. We explained 9. V. Buchholtz, T. Pöschel, in Friction, Arching and Contact
this effect by the formation of force chains. Particles which Dynamics, edited by D.E. Wolf, P. Grassberger (World
are members of one or more force chains experience a sig- Scientific, Singapore, 1997), p. 265.
nificantly larger average force, thus substantially enhanc- 10. J.A.C. Gallas, H.J. Herrmann, T. Pöschel, S. Sokolowski,
ing the fragmentation probability, as compared to parti- J. Stat. Phys. 82, 443 (1996).
cles which do not belong to any force chain. 11. G.A. Kohring, Physica A 195, 1 (1993).
From this we concluded that the experimentally ob- 12. C. Salueña, S.E. Esipov, T. Pöschel, S. Simonian, SPIE
served and numerically achieved spatial pressure and frag- 3327, 19 (1998).
mentation distributions are crucially determined by the 13. C. Salueña, S.E. Esipov, T. Pöschel, Phys. Rev. E 59, 4422
formation of force chains. Direct collisions of particles are (1999).
of minor importance for comminution. Without the ex- 14. G.A. Kohring, S. Melin, H. Puhl, H.-J. Tillemans, W.
istence of force chains ball mills would work much less Vermöhlen, Appl. Mechanics Eng. 124, 273 (1995).
efficient. In the context of a ball mill efficiency optimiza- 15. T. Yokohama, K. Tamura, G. Jimbo, Intern. Chem. Eng.
tion it would be intriguing to investigate the possibility 34, 611 (1994).
of enhancing the formation of force chains, e.g. by modi- 16. T. Yokohama, K. Tamura, G. Jimbo, in Proceedings
fications of the geometrical properties of mills, additional 8. European Symposium on Comminution, edited by E.
application of vibration etc. Forssberg (Elsevier, Amsterdam, 1996), p. 413.
The method developed in this article may be elabo- 17. P.K. Songfack, S. Agrawala, B.K. Mishra, R.K. Rajamani,
rated in many different aspects: in Proc. 2nd Intern. Conf. on Discrete Element Methods,
edited by J.R. Williams, G.G.W. Mustoe (IESL, Cam-
– The numerical investigations were done in two dimen- bridge, 1993), p. 205.
sions. Ball mills are three-dimensional devices. 18. A tumbling ball mill is a tubular vessel which is kept rotat-
– Axial transport of the material was not simulated. ing around its own axis. The grinding medium, here balls,
– All particles, including fragments, were spherical. is elevated by the vessel wall and then rolls, slides, or falls
– Only a single-level process was simulated. Modern ball toward the lowest point of the vessel. The feed material
mills frequently operate in multi-level mode. is trapped between the media particles and, thus, ground
– Adhesive forces were not considered in the force law. by their movement. In case the grinding media is the feed
– Variation of material properties which originates from material itself one speaks of autogeneous grinding. Even
an increase of temperature due to comminution was though the autogeneous comminution in a tumbling mill
not considered. is considered throughout this paper we, nevertheless, will
use the simple term ball mill.
We expect that molecular dynamics will rapidly de-
19. H. Hertz, J. Reine Angewandte Math. 92, 156 (1882).
velop towards an important tool for the powder technol- 20. L.D. Landau, E.M. Lifschitz, Elastizitätstheorie (Akademie
ogy machinery. Due to the rapid progress in hardware and Verlag, Berlin, 1983).
software technology, including parallel algorithms, it will 21. P.A. Engel, Impact Wear of Materials (Elsevier, Amster-
soon be possible to construct and optimize powder tech- dam, 1976).
nology facilities. 22. M.P. Allen, D.J. Tildesley, Computer Simulations of Liq-
uids (Clarendon Press, Oxford, 1987).
23. J. Schäfer, S. Dippel, D.E. Wolf, J. Phys. France I 6, 5
(1996).
The authors thank S. Bernotat, R. Drögemeier, E. Gommeren, 24. N.V. Brilliantov, F. Spahn, J.-M. Hertzsch, T. Pöschel,
G. Gudehus, H.J. Herrmann, L. Rolf, J. Schwedes and R. Phys. Rev. E 53, 5382 (1996).
Weichert for discussion and for providing relevant literature. 25. R. Ramı́rez, T. Pöschel, N.V. Brilliantov, T. Schwager,
The work was supported by German Science Foundation Phys. Rev. E 60, 4465 (1999).
through project Po472/3-2. 26. D.C. Rapaport, J. Comp. Phys. 34, 184 (1980).
182 The European Physical Journal B

27. T. Schwager, T. Pöschel, Phys. Rev. E 57, 650 (1998). bei der Naßmahlung in Kugelmühlen, Ph.D. thesis, Techn.
28. S. Luding, S. McNamara, Granular Matter 1, 113 (1998). Universität München, 1985.
29. R. Weichert, Zement Kalk Gips 45, 1 (1992). 45. The term pressure distribution does not relate to the spa-
30. R. Weichert, Int. J. Mineral Processing 22, 9 (1988). tially averaged force vector field, but to the absolute value
31. W. Weibull, The Phenomenon of Rupture in Solids. Proc. of forces acting between particles. Strictly speaking, this
151 and 153. Ing. Vetensk. Akad. (Stockholm, 1939). quantity has the unit of a force, however, to avoid confu-
32. S. Bernotat, K. Schönert, Size reduction, Ullmann’s Ency- sion with the averaged force vector field we adopt the term
clopedia of Industrial Chemistry B2, 1 (1988). pressure from engineering literature.
33. F. Kun, H.J. Herrmann, Int. J. Mod. Phys. 7, 837 (1996). 46. J. Rajchenbach, Phys. Rev. Lett. 65, 2221 (1990).
34. F. Kun, H.J. Herrmann, Comp. Meth. Appl. Mech. Eng. 47. J. Duran, T. Mazozi, S. Luding, E. Clément, J.
138, 3 (1996). Rajchenbach, Phys. Rev. E 53, 1923 (1996).
35. R. Schuhmann, Am. Inst. Mining, Met., Petrol. Engrs. 48. R. Behringer, Nonlinear Science Today 3, 1 (1993).
49. C.H. Liu, S.R. Nagel, D.A. Schecter, S.N. Coppersmith,
Tech. Publ. No. 1189 (1940).
S. Majumdar, O. Narayan, T.A. Witten, Science 269, 513
36. P. Rosin, E. Rammler, J. Inst. Fuel 7, 29 (1933).
(1995).
37. J.J. Gilvary, J. Appl. Phys. 32, 391 (1961).
50. P. Dantu, Ann. Ponts Chausses 4, 144 (1967).
38. J.J. Gilvary, B.H. Bergstrom, J. Appl. Phys. 32, 400 51. S. Schöllmann, Phys. Rev. E. 59, 889 (1999).
(1961). 52. S.E. Esipov, T. Pöschel, J. Stat. Phys. 86, 1385 (1997).
39. J.A. Herbst, M. Siddique, K. Rajamani, E. Sanchez, Soc. 53. D.E. Wolf, in Computational Physics, edided by K.H.
Min. Eng. AIME Trans. 272, 1945 (1981). Hoffmann, M. Schreiber (Springer, Heidelberg, 1996),
40. J.A. Åström, H.J. Herrmann, Eur. Phys. J. B 5, 551 p. 64.
(1998). 54. D. Howell, B. Miller, C. O’Hern, R.P. Behringer, in Fric-
41. A.N. Kolmogorov, Doklad. Akad. Nauk. S.S.S.R. 31, 99 tion, Arching and Contact Dynamics, edided by D.E. Wolf,
(1941). P. Grassberger (World Scientific, Singapore, 1997), p. 293.
42. P.R. Halmos, Ann. Math. Stat. 15, 182 (1944). 55. F. Radjai, M. Jean, J.J. Moreau, S. Roux, Phys. Rev. Lett.
43. B. Epstein, J. Franklin Inst. 244, 471 (1947). 77, 274 (1996).
44. J. Leluschko, Untersuchungen zum Einfluß der Flüssigkeit

View publication stats