Starting from a monodisperse 
pile with bumpy bottom and
= 97, see Fig. 2(b),
we change the particle size of each particle slowly to the
diameter 
, where 
is a random number
homogeneously distributed in
the interval [-r/2,r/2]. We present the vertical stress in
Fig. 8, for simulations with r = 2/3000 (a), 2/300 (b),
and 1/30 (c).
We plot the result of one run (solid line) and compare it with
the monodisperse case (dashed line) and the average over 40 runs (a)
or 100 runs (b) and (c) (symbols). The fluctuations in stress
increase with increasing r. In fact we observe fluctuations
much larger
than the total stress for the monodisperse pile. With increasing
r the shape of the averaged vertical stress changes in the
center from a hump [see r=2/3000],
to a dip [see r = 1/30]. The averaged stress in
Fig. 8(c) is similar to the stress obtained (after many
averages) from a cellular automaton model for the stress propagation
in the presence of randomly opened contacts [17].
Figure 8:
Vertical stress V(1), in row M=1, vs. X for a
pile with bumpy bottom and 
. The particle diameter is
homogeneously distributed in the interval 
.
The values of r are r=2/3000 (a), r=2/300 (b), and
r=1/30 (c). The dashed line
gives the result of Fig.2(b) with r=0 and .
The solid line gives the result of one run and the symbols
correspond to an average over 40 runs for (a), or 100 runs and three
particles for (b) and (c).
In Fig. 9(a) we give the contact network of one run as
presented in Fig. 8(c). The line thickness indicates the
magnitude of forces active at a contact. In Figs. 8(b)
we present a part of this contact network, and in Figs. 8(c)
we plot the principal axis of the stress tensor for the same part.
In Figs. 9(a) and (b)
each line represents the normal direction of one contact and
each particle center is thus situated at the meeting point of several
lines. Note that in Fig 9(a) some particles inseide the
pile have no contacts to their above neighbors, i.e. they lie below
an arch. Comparing the contact network (b) with the stresses in (c)
one may again relate the structure of the contcat network to the
angle 
and the ratio of the principal axis, as discussed
above in Sec. ivA2.
Finally, we calculate the probability distribution P for
vertical stresses V.
We average only over the lowermost row M=1 and also neglect
the outermost particles of this row. In detail we average
all particles 
[counting from the lower
left end to the right] over 100 runs.
Since the stress in row M=1 is a function of X in the
case of r=0, we scale the stresses for r>0 by the
stresses found for r=0, i.e. we use the scaled stress
with the minimum of all T, obtained from particles
which are shielded and thus feel only their own weight, see Fig.\
9.
Note that this occurs frequently, even inside the pile.
We checked that the probability distribution of T does not depend
on the specific choice of the interval, i.e. we also
averaged over a smaller interval 
, or
over particles in row M=2 with 
, and
we found no difference besides fluctuations.
Figure:
(a) Contact network of one pile from Fig. 8(c).
The line thickness indicates the magnitude of the contact force.
(b) Part of the contact network from (a).
(c) Principal axis of stress from the
same simulation as in (b).
We plot the distribution function P(T) in Fig. 10. The dashed line in Fig. 10(a) shows a power law for small stresses T, while the dotted line in Fig. 10(b) shows an exponential decay for large T.
Figure:
(a) Double logarithmic
plot of the probability distribution of small vertical stresses
T in row M=1. We skip 10 particles at the right and the left
and average over 100 runs. see Eq. 12.
(b) Log-linear plot of the probability distribution for large stresses
from the same data as in (a).
Thus our results are in agreement with the theoretical predictions of Ref.\ [5], and the numerical findings of Ref. [6], at least for large T, see Fig. 10(b). Note that the probabilty to find large T is greater for r=1/30 than for r=2/300, corresponding to stronger fluctuations. For small T we find a power law with exponent -1/2 for both values of r, see Fig. 10(a).
