# Stress in frictionless granular material: Adaptive Network Simulations

## I Introduction

The fragility of granular matter is a longstanding preoccupation of engineers  and a recent preoccupation of physicists. By granular matter we mean a static assembly of hard, spheroidal grains whose contact forces may be compressive but not tensile. Thus granular matter is noncohesive. Mohr and Coulomb recognized a fundamental continuum consequence of the noncohesive state. There can be no co-ordinate system in which the shear stress exceeds some fixed multiple $`\mu `$ of the normal stress . This “Mohr-Coulomb” condition limits the stresses that a granular material can support, and thus amounts to a form of fragility. When this condition is violated, building foundations settle and embankments slip.

Modern civil engineering practice views the stress field in a granular medium as divided into elastic and plastic zones. The stresses in these plastic zones are at the Mohr-Coulomb limit, and thus these zones are at the margin of stability. The stress in the elastic zones is within the bounds of stability and thus the stress here is transmitted as in an elastic body.

Recently attention has turned to the microscopic origin of the macroscopic fragility of granular media. The microscopic pattern of contact forces and bead motions shows strong local heterogeneity and history dependence. The history of prior motion in a region clearly influences the way it transmits forces. The prior motion may affect the $`\mu `$ coefficient in the Mohr-Coulomb law, the elasticity tensor, or further constitutive properties. The question is, for a given history of relaxation to a static state, how are these forces transmitted and what range of forces can be supported. The transmission of forces can be expressed as a linear-response property of a granular pack. An infinitessimal force is added to the bead at position $`𝐱_0`$ and the corresponding incremental force on a contact at $`r`$ is determined. For sufficiently small perturbations of a finite pack, this linear response function G is well defined. It depends on the shapes and sizes of the beads, their frictional properties, and how the pack was constructed.

If the perturbing forces become too large, motion occurs. Beads shift their positions and form new contacts. This motion may be reversible, so that the beads return to their original positions when the perturbation is removed. This motion may also be irreversible, with the positions altered after removal of the perturbation. The thresholds for reversible and for irreversible motion are fundamental ways to characterize the nonlinear response of the pack. For any given pack there is a weakest perturbing force distribution that causes motion. This threshold force may go to zero as the size of the pack grows. Many simulations have sought to characterize the above features of a granular pack. These studies model the system in a realistic way that requires detailed specifications and many parameters. This detail makes it difficult to discern which observed features are inescapable consequences of the granular state, and which are properties of the particular realization. In this study we take the opposite approach, sacrificing realism for the sake of simplicity. We seek the simplest system that shows the instabilities of noncohesive material. Thus our system consists of frictionless, spherical beads, which have been deposited into a container one at a time and not moved thereafter. Such a system develops tensile contacts.

To avoid these contacts, we must define some motion that evolves the pack to a more stable state. Again we choose a procedure favoring simplicity rather than realism. We seek the stable state attainable with minimal disturbance from the initial state. Accordingly our procedure does not move the beads, but rather removes and adds contacts one at a time in order to attain a stable contact network. In this sense our simulation is an adaptive network.

This network is isostatic: the contact forces are determined from the applied forces solely through the force equilibrium of each bead, without reference to bead displacements or material deformation.

Our adaptive method demonstrates that frictionless granular materials can be mechanically robust. For a given load, the simulation converges to a state of no tensile contacts. A change in the applied load of order unity can be applied with only minor shifts in the contacts. Further, the three components of the stress in two dimensions obey a constitutive law of the “null stress” type: a weighted sum of the three components vanishes, the weights depending on the packing but not on the loading.

## II Response Function

Our system is a set of spherical beads, whose radii are chosen randomly within a moderate range. These beads are supported on one side, called the bottom, with a layer of fixed spheres. The width $`w`$ of the system is much larger than a bead. The beads are arranged densely in this space up to a height $`h`$. A fixed downward force $`F_0`$ is applied to each bead lying at the upper surface. We choose a configuration of beads that is mechanically stable: the normal forces acting at the bead contacts oppose the applied forces $`F_0`$ and prevent motion. Initially, we allow for tensile (negative) contact forces.

We label the $`N`$ beads by an index $`\alpha `$. Then we may denote the contact force from bead $`\alpha `$ to bead $`\beta `$ by the scalar $`f_{\alpha \beta }`$. Newton’s Third Law dictates that for all $`\alpha `$, $`\beta `$, $`f_{\beta \alpha }=f_{\alpha \beta }`$. The $`N_c`$ contact forces are constrained by the requirement that the total vector force on each bead vanish. In $`d`$ dimensions, there are evidently $`dN`$ such constraints. The bead positions $`𝐱_\alpha `$ are likewise constrained by the geometrical condition that the distance between two contacting beads $`\alpha `$ and $`\beta `$ must be the sum of their radii $`r_\alpha +r_\beta `$. There are $`N_c`$ such constraints for the $`dN`$ quantities $`𝐱_\alpha `$. If all these constraints are independent, the number $`N_c`$ of contacts must be exactly $`dN`$. Then the system is isostatic: the $`dN`$ force balance equations are just sufficient to determine the $`N_c`$ contact forces. The equations of force balance may be written

$$\underset{\beta (\alpha )}{}f_{\alpha \beta }\widehat{n}_{\alpha \beta }=𝐅_\alpha $$
(1)   
Here $`\beta (\alpha )`$ denotes the set of contacting neighbors of bead $`\alpha `$, $`𝐅_\alpha `$ is an external force applied to this bead, $`f_{\alpha \beta }`$ and the unit vector $`\widehat{n}_{\alpha \beta }`$ represent the magnitude and (fixed) direction of the contact force between the two beads $`\alpha `$ and $`\beta `$. Since all the above equations are linear, the response of the system to a given external forcing is determined by the response function $`𝐆`$:

$$f_{\alpha \beta }=𝐆(\alpha \beta |\gamma )𝐅_\gamma $$
(2)   
The response function $`𝐆`$ determines not only the response to an external force but also the global displacement field associated with local geometrical perturbation of the network. In order to see this, let us assume that the packing is subjected to external forcing $`𝐅_\gamma `$. Then we relax exactly one of the $`N`$ geometric constraints, and change infinitesimally the distance between two contacting beads, $`r_{\alpha \beta }|𝐱_\beta 𝐱_\alpha |`$. As long as the connectivity of the network does not change, its motion is non-dissipative. This means that the work done to distort the packing locally, $`\delta r_{\alpha \beta }f_{\alpha \beta }`$ is the work against external forces, i. e.

$$\delta 𝐱_\gamma 𝐅_\gamma =\delta r_{\alpha \beta }f_{\alpha \beta }=\delta r_{\alpha \beta }𝐆(\alpha \beta |\gamma )𝐅_\gamma $$
(3)   
The above equation should be valid for any set of external forces; hence,

$$\delta 𝐱_\gamma =𝐆(\alpha \beta |\gamma )\delta r_{\alpha \beta }$$
(4)   
We conclude that $`𝐆`$ is the response function both for contact force and the displacement field. Note that the displacement discussed here is not due to deformation of the beads. It corresponds to a “soft mode” that preserves the distances between all contacting beads other than the perturbed contact.

In a general case, finding the response function for a given configuration is a non-local problem, which requires solving a set of linear Eq. (1). The task becomes much easier for the case of sequential packing. This is created by adding one bead at a time. The requirement of mechanical stability implies that any newly–added bead has exactly $`d`$ “supporting” contacts (in $`d`$-dimensional space). If all the contacts were permanent and this $`d`$-branch tree structure were not perturbed by the future manipulations, the response function might be found by a simple unidirectional projection procedure. Indeed, since there are exactly $`d`$ supporting contacts for any bead in a sequential packing, the total force $`\stackrel{~}{𝐅}_\alpha `$ including external force $`𝐅_\alpha `$ and that applied from the supported beads can be uniquely decomposed onto the corresponding $`d`$ components, directed along the supporting unit vectors $`𝐧_{\alpha \gamma }`$. This gives the values of the supporting forces. The $`f`$’s may be compactly expressed in terms of a generalized scalar product $`\mathrm{}|\mathrm{}_\alpha `$:

$$f_{\alpha \beta }=\stackrel{~}{𝐅}_\alpha |\widehat{n}_{\alpha \beta }_\alpha $$
(5)   
The scalar product $`\mathrm{}|\mathrm{}_\alpha `$ is defined such that $`\widehat{n}_{\alpha \beta }|\widehat{n}_{\alpha \beta ^{}}_\alpha =\delta _{\beta \beta ^{}}`$. for the supporting contacts $`\beta `$, $`\beta ^{}`$ of bead $`\alpha `$. In general, it does not coincide with the conventional scalar product. The resulting response function, $`𝐆(\alpha \beta |\gamma )`$, can be calculated as the superposition of all the projection sequences (i.e. trajectories), which lead from bead $`\gamma `$ to the bond $`\alpha \beta `$:

$$𝐆(\alpha \beta |\gamma )=\underset{(\gamma \alpha _1\mathrm{}\alpha \beta )}{}|\widehat{n}_{\gamma \alpha _1}_\gamma \widehat{n}_{\gamma \alpha _1}|\widehat{n}_{\alpha _1\alpha _2}_{\alpha _1}\mathrm{}\widehat{n}_{\alpha _k\alpha }|\widehat{n}_{\alpha \beta }_\alpha $$
(6)   
Here the summation is done over all the trajectories $`(\gamma \alpha _1\mathrm{}\alpha _k\alpha \beta )`$ such that any bead in the sequence is a supporting neighbor of the previous one.

## III Adaptive network simulation.

For a large enough system, sequential packing is not compatible with the requirement of non-tensile contacts. Anytime when this requirement is violated, a rearrangement occurs and system finds a “better” configuration. One might expect that this would make the problem of force propagation a dynamic one. However, it is possible to limit oneself to a purely geometrical consideration, following the ideas of the previous section.

Suppose $`\alpha \beta `$ is a “bad bond”, whose contact force is negative (tensile). This means that the network would move in such a way that the two beads, $`\alpha `$ and $`\beta `$ are taken apart. In other words, the soft mode associated with the perturbation of $`\alpha \beta `$ bond is activated, and for small enough displacements all the beads move in accordance with Eq. (4). The motion stops when a replacement contact is created, i.e. when a gap between any two neighboring beads closes.

In this work, we limit ourselves to this linear approximation. It should be understood that Eq. (4) is correct only for infinitesimal displacements, and in a general case one should account for the evolution of the response function in the course of the rearrangement. We avoid the problem of changing G by permitting only infinitessimal motion in the model. We imagine that the “bad bond” gets deactivated, and it is replaced with a rigid “strut” between two neighboring beads that were not in contact in the previous configuration. There is a natural choice for where the strut should be placed to cause minimal disturbance. Each pair of non-contacting neighbors $`\gamma \delta `$ has a gap $`r_{\gamma \delta }r_\gamma r_\delta `$. When the contact $`\alpha \beta `$ is removed, the distance $`r_{\alpha \beta }`$ is allowed to change; this change alters the gaps of other neighbors $`\gamma \delta `$ as specified by Eq. 4. Extrapolating this linear-response equation, motion $`\delta r_{\alpha \beta }`$ required to close gap $`\gamma \delta `$ is given by

$$\delta r_{\alpha \beta }=\frac{r_{\gamma \delta }r_\gamma r_\delta }{\widehat{n}_{\gamma \delta }[𝐆(\alpha \beta |\gamma )𝐆(\alpha \beta |\delta )]}$$
(7)   
(For many choices of $`\gamma \delta `$ the required $`\delta r_{\alpha \beta }`$ is infinite since the $`\alpha \beta `$ contact has no effect on the $`\gamma \delta `$ gap.) Using this formula, we identify the gap $`\alpha ^{}\beta ^{}`$ which would require the smallest change of $`r_{\alpha \beta }`$ in order to close, and we link this pair by a strut.

After we have found the replacement bond, the modified response function can be found without solving the whole set of force balance equations (1)! We denote the response function for the initial packing as $`𝐆_0`$. The new response function G must be such that there is no longer a contact force $`f_{\alpha \beta }`$. In general there is such a force in the initial packing. However, we may alter this unwanted force by adding an external force to some other bead. We choose to add external forces to beads $`\alpha ^{}`$ and $`\beta ^{}`$ that mimic a contact force: the two forces are equal, opposite and directed along the unit vector $`\widehat{n}_{\alpha ^{}\beta ^{}}`$ joining them, with a strength denoted $`f_{\alpha ^{}\beta ^{}}`$. Our choice of the replacement pair guarantees that the effect of $`\alpha ^{}`$ or $`\beta ^{}`$ on the $`\alpha \beta `$ contact is non-vanishing. Then the force on this contact is given by

$`f_{\alpha \beta }={\displaystyle \underset{\gamma }{}}𝐅_\gamma 𝐆_0(\alpha \beta |\gamma )`$ (8)   
$`+f_{\alpha ^{}\beta ^{}}[𝐆_0(\alpha \beta |\alpha ^{})𝐆_0(\alpha \beta |\beta ^{})]\widehat{n}_{\alpha ^{}\beta ^{}}.`$ (9)   
We may make this $`f_{\alpha \beta }`$ vanish by a proper choice of the external force, $`f_{\alpha ^{}\beta ^{}}=_\gamma 𝐅_\gamma 𝐆(\alpha ^{}\beta ^{}|\gamma )`$, where $`𝐆(\alpha ^{}\beta ^{}|\gamma )`$ is the new response function, found by requiring that $`f_{\alpha \beta }`$ vanish in Eq. (8):

$$𝐆(\alpha ^{}\beta ^{}|\gamma )=\frac{𝐆_0(\alpha \beta |\gamma )}{\widehat{n}_{\alpha ^{}\beta ^{}}[𝐆_0(\alpha \beta |\alpha ^{})𝐆_0(\alpha \beta |\beta ^{})]}.$$
(10)   
A contact force on an arbitrary contact is now determined by a combination of external forces $`𝐅_\gamma `$ and the above–determined $`f_{\alpha ^{}\beta ^{}}`$. This results in the following expression for $`𝐆(\lambda \mu |\gamma )`$ ($`\lambda \mu `$ other than $`\alpha ^{}\beta ^{}`$);

$`𝐆(\lambda \mu |\gamma )=𝐆_0(\lambda \mu |\gamma )`$ (11)   
$`𝐆_0(\alpha \beta |\gamma ){\displaystyle \frac{\widehat{n}_{\alpha ^{}\beta ^{}}[𝐆_0(\lambda \mu |\alpha ^{})𝐆_0(\lambda \mu |\beta ^{})]}{\widehat{n}_{\alpha ^{}\beta ^{}}[𝐆_0(\alpha \beta |\alpha ^{})𝐆_0(\alpha \beta |\beta ^{})]}}.`$ (12)   
This prescription gives the response of the pack to a wide class of contact replacements. The prescription does not require either the initial or the final state to be stable: it allows tensile contacts. Using the contact replacement procedure we may investigate the stability of a pack systematically. Our algorithm has two major stages: network preparation and its “mutation” via the contact replacement scheme. By repeating this adaptive procedure sufficiently many times, one may hope to get the stable configuration without tensile forces (for a given loading), just like the real system would do. There is a possibility that the present geometry–preserving algorithm could not stabilize the network. For instance if the tangential component of the surface force is strong enough, it is expected to initiate a macroscopic avalanche, as in a sandpile with slope exceeding the critical angle. This class of rearrangements is beyond the capabilities of our connectivity mutation scheme. This circumstance has even certain advantages: we can determine the critical slope from our simulations as the direction of the surface force at which the algorithm stops working. It should be emphasized that our algorithm can easily be modified to incorporate the change in the network geometry. The major reason why we use the above geometry-preserving (“strut”) approximation is its much higher computational effectiveness.

## IV Simulation details and results

### A Method

We begin by creating a two-dimensional sequential pack of variable-sized discs by adding them one by one. The studied system has the following parameters: polydispersity 10% (bead radii from 1 to 1.1), number of beads, $`N`$ from 250 to 500 (the major limitation is the computation time). Although there is no gravitational force acting on the beads in these simulations, the statistics of the packing can be varied by changing the “pseudo-gravity” direction, $`\widehat{g}`$. Namely, while adding a bead to the packing we require $`\widehat{g}`$ to be directed between the two supporting contacts. Simultaneously, we calculate the response function, by using the sum-over-trajectories formula, Eq. (6).

Then we apply certain load to beads on the surface. In the studied cases, the forces applied to the surface beads were all the same, with the only principal variable being the ratio of the two components $`f_x/f_y`$ ($`y`$ is the vertical direction, the surface in average is parallel to $`x`$ due to the periodic boundary conditions in the horizontal direction). As long as the response function is given, we may find all the contact forces for a given load. As we have found, tensile contacts appear within a few beads from the surface. We analyze the sign of the contact forces one by one, from top to bottom (in the order opposite to the one in which the beads were originally deposited). When a tensile contact is encountered, we follow the contact replacement procedure described in the previous section: find the new bond and modify the response function to account for the connectivity change. Now we repeat the procedure again starting from very top until there is no tensile contact left in the system Figure 1 shows a sequence of four typical steps in this procedure. Evidently the removed and the added contact may be far apart.

It sometimes happens that this prescription does not remove the tensile contacts: the removal of a tensile contact continues to generate others. In this case we may modify our procedure for selecting the next tensile contact to remove. For example we may select the strongest tensile contact instead of that tensile contact having the largest sequence number. Such alternative prescriptions seem to have little effect on the force network, as discussed in the next section.

### B  Variability and reproducibility

While our bond replacement procedure mimics the way the real system should rearrange, our choice of the “bad bond” to be replaced is far more arbitrary. For instance, instead of checking the sign of the contact forces one by one from top to bottom, we could go the other way, or try to replace the contact bearing the largest negative force first. Neither of these prescriptions is very realistic; however, we find that the results are insensitive to the procedure. In order to probe this sensitivity, we compared the results from two different “annealing” procedures. The first was the procedure described in the previous subsection. The second procedure is as follows. We perform exactly the same one-by-one check, as in the previous case, but do not remove a negative bond unless the magnitude of the force exceeds certain tolerance threshold. When no more bonds remain in the system for a given threshold level, we reduce the tolerance and repeat the procedure. The threshold plays a role of temperature: keeping it finite allows us to deviate from the target (non-tensile) state of the system and explore its vicinity at the configuration space. The second annealing algorithm converges considerably faster than the zero-tolerance one. For instance, it took us from 1500 to 3000 iterations to complete the original algorithm with 500 beads, while the annealing procedure reduced the needed time to approximately 500-1000 steps. Interestingly, the variation of the convergence time is of order of the time itself.

We have found that the contact configuration resulting from the annealing procedure does differ from the one generated by the zero-tolerance algorithm (see Figure 2). However, we did not detect any statistically–significant variation of the ensemble–averaged properties of the final state obtained with the two methods. These properties included the average stress and the contact force probability distribution function (PDF), presented below.

### C Macroscopic constitutive equation

One of the crucial results of the simulation is that our geometry-preserving adaptive network algorithm does converge for a considerable range of force direction. It stops working when $`|f_x/f_y|`$ approaches $`0.6`$ (for packing prepared at vertical pseudogravity $`\widehat{g}`$. This suggests that the critical slope for the frictionless packing is about 30 degrees, consistent with simple theoretical arguments and some experiments  Note that this slope may considerably exceed the angle of repose in dynamic experiments and simulations because of hysteresis associated with the lack of damping in the frictionless system. Presumably, the critical slope can be observed by quasi-static tilting of a zero-slope packing.

Another interesting observation is that the eventual connectivity of the packing is not too different from the original sequential packing. For example, the 500-bead system needs up to 3000 iterative steps (rearrangements) to find the stable state, and yet only 150 out of 1000 contacts ($`15\%`$) in the final configuration are different from the original network. This provides us with a solid background for using the sequential packing as the zero-order approximation of the real network. This was one of the major hypothesis used in our earlier work  to derive the constitutive equation of frictionless granular packing.

One more hypothesis, used for derivation of the macroscopic equation for stress is the mean field decoupling Ansatz. The average stress in a region of a sequential packing can be written

$`\sigma ^{ij}(𝐱)=`$ (13)   
$`{\displaystyle \underset{\alpha }{}}{\displaystyle \underset{\beta (\alpha )}{}}\delta (𝐱_\alpha 𝐱)\stackrel{~}{𝐅}_\alpha |𝐧_{\alpha \beta }_\alpha n_{\alpha \beta }^in_{\alpha \beta }^jr_{\alpha \beta },`$ (14)   
The sum $`\beta `$ is over the beads that support the bead $`\alpha `$. Our mean-field hypothesis consists in assuming that the force-related part of this average is independent of the geometrical part. We define

$$𝐟(𝐱)\underset{\alpha }{}\delta (𝐱_\alpha 𝐱)\stackrel{~}{𝐅}_\alpha |,$$
(15)   
and

$$\widehat{\tau }\underset{\beta (\alpha )}{}|𝐧_{\alpha \beta }_\alpha n_{\alpha \beta }^in_{\alpha \beta }^jr_{\alpha \beta }.$$
(16)   
Then our mean-field assumption amounts to the statement

$$\sigma ^{ij}=𝐟\widehat{\tau }$$
(17)   
The mean field hypothesis can be directly checked for unperturbed sequential packing, were both fields $`𝐟`$ and $`\widehat{\tau }`$ are well–defined. The results of such a check are represented on Figure 3. The exact values of various stress components are shown to be in an excellent agreement with their evaluation based on the mean-field Ansatz. We conclude that the mean field is a very good approximation at least for non-adaptive sequential packing.

As long as the rearrangements are switched on, there is no obvious way to define the concept of supporting neighbor, and therefore the “force from subsequent beads”, $`𝐟`$ is ill-defined as well. However, a more general meaning of constitutive Eq. (13), is that the stress is parameterizable with some vector $`𝐟`$ , and the third-rank material tensor $`\widehat{\tau }`$ establishes this parameterization. We now take $`\widehat{\tau }`$ corresponding to the original pre-rearrangement sequential packing and probe the Ansatz by checking whether the total stress (after the adaptive procedure) can be expressed as $`\widehat{\tau }𝐟`$. In other words, we compare the only unknown component of the stress $`\sigma _{xx}`$ (two other are given by the boundary conditions) with its theoretical value obtained from $`\widehat{\tau }`$. We performed this check for two different classes of packing, corresponding to different directions of “pseudogravity”, $`g_x/g_y=0`$ and $`0.2`$. The agreement between the two curves is surprisingly good as long as the direction of the applied force does not deviate too much from the preparation conditions (i. e. from the pseudogravity vector), see Figure 4.

### D  Response function

As noted above, our system transmits forces in accord with a null-stress constitutive property. Given the null-stress law, one may infer the corresponding response function G. The force transmission is transmitted from a point source according to a wave-like equation. In a medium where all non-vertical directions are equivalent, the force should propagate downward along slanting characteristic lines, whose slope is dictated by the only parameter in the null-stress law. The responding region lies within the “light cone” bounded by these lines. In two dimensions, the response consists of two delta-functions traveling along the light cone. Disorder is expected to scatter the wave solutions of the pure system, thus resulting in a widening of the delta-peaks. This scattering could be sufficient to create qualitative new mesoscopic behavior from localization effects. Our simulated system showed strong influences from disorder, as illustrated in Figure 5. Because there can be no vertical-force response at the top of the system, we observe a global anisotropy, with stronger responses below the source than above it. The response is also strongly heterogeneous.

Our simulations allow us to perform ensemble averaging of the response stress field. Figure 6 shows the results of such averaging over 600 realizations of the network. As the perturbation propagates deep into the sample, the response function gets a two–peak shape, in a good agreement with the null-stress law. As expected, the peaks are broadened by the disorder, and one cannot resolve them immediately below the source. Another important observation, which also supports the null-stress approach, is that the average response is virtually zero above the source. Note that we have studied only linear response of the system, so that the perturbation did not change the contact network. This need not to be the case in the experiments involving strong local perturbations 

### E Contact force distribution function

Our simulation allows us to address yet another interesting and widely-discussed problem: the statistics of contact force. Recent experiments indicate that this distribution can be well approximated as exponential, that is, it is considerably wider than a naively-expected Gaussian. This is related to the strong heterogeneities of the mesoscopic stress in granular matter: it appears to be localized to string-like structures known as force chains.

In the initial sequential packing there is no constraint on the sign of the contact force, and its amplitude appears to grow indefinitely with the packing depth. After the rearrangements, there are no negative forces in the system, and therefore their amplitude cannot grow forever(the total transmitted force is fixed). Figure 7 shows the spatial distribution of the contact forces in the systems of two different degrees of polydispersity after the adaptive stage is completed. One can clearly see that our simulations are at least in qualitative agreement with experiment: it is easy to identify the force chains in both cases.

We were also able to make a quantitative comparison between the simulations and experiments. Figure 8 shows the probability distribution function of the contact force taken from our simulations of almost monodisperse system. It apparently agrees with the exponential histogram observed experimentally. An insight into the origin of this exponential behavior is given by the “q–model” due to S. Coppersmith et al. The further discussion of this intriguing result will be published elsewhere 

## V Conclusion

In the study of granular materials, clearcut confirmation of theories has been elusive. One predicted feature of great interest is the null-stress constitutive law postulated by Wittmer et al. We have verified that null-stress behavior occurs in a simplified granular system embodying disorder, perfect rigidity and cohesionless contacts. We have measured the free parameter in the null-stress law for several situations. We have confirmed the validity of our major assumptions used for microscopic foundation of this constitutive law. Our simulation also allowed us to compute directly the ensemble–averaged response function, thus providing an additional check for the adequacy of the null-stress approach.

The simulation method has further interesting features. It demonstrates that stable configurations of isostatic force networks can be found without changing the positions of the nodes. It also reveals order-unity variability in the microscopic force distribution resulting from the relaxation process. Finally, it shows strongly heterogeneous response to point forces—much stronger than that of the geometric contact network. This suggests strong multiple-scattering features in the force propagation.

The observed exponential probability distribution function for the contact force is in a good agreement with the experiments. Since there is also an indirect experimental support for the null-stress law, our choice of the system (hard frictionless spheres) appears to be an adequate simplification to capture the basic physics of granular rigidity. The further simplifications, such as the fixed–geometry adaptive algorithm provide an effective tool for the future studies of this problem. This may include a study of non-linear response of the system to large localized perturbation, effects of polydispersity, and history dependence of the response.

## Acknowledgement

The authors thank R. Ball, S. Coppersmith, D. Mueth, H. Jaeger, S. Nagel, and J.Socolar for valuable discussions. Likewise, we thank the participants in the Jamming and Rheology program of the Institute of Theoretical Physics in Autumn 1997. This work was supported in part by the National Science Foundation under Award numbers PHY-94 07194, DMR-9528957, DMR-9975533 and DMR 94 00379.