## Abstract

If sand is continually and slowly added onto a table, a pile will grow and reach some steady state, whereupon any further sand added to the pile will fall off the table, leaving the surface of the pile unchanged. A construction that gives the surface of a unique steady pile for a table of any given shape was proposed by Arpad L. Nadai. Here, we show that Nadai's construction can be justified mathematically on the basis of a few plausible assumptions on the mechanical behaviour of sand. We also show that the surface of the steady pile given by the construction is the maximal topography supported on the domain of the table; that the surface may contain smooth patches, valleys, peaks and ridges; that any smooth patch is developable; that any valley is smooth, straight and without tributaries, and has a slope equal to the angle of repose of the sand; that any peak is sharp; and that any ridge is sharp and has a slope less than the angle of repose of the sand. Lastly, we perform simple experiments and verify that actual steady piles look similar to those given by the construction, with a precision that depends on the type of sand and how carefully the sand is poured to form a pile.

## 1. Introduction

Piles of granular materials such as sand appear to have slopes bounded by the *angle of repose* (figure 1*a*), which is a property of granular materials. The angle of repose arises from a competition between gravity and the resisting forces due to friction, cohesion and particle roughness; it represents the slope for which the resisting forces are just strong enough to check the gravitational force (Jaeger *et al.* 1996). Thus, if a small lump of sand is added to the surface of a pile at the angle of repose, the lump causes the slope to locally exceed the angle of repose. As a result, sand starts to flow down the slope. The flow is confined near the surface (figure 1*b*) and is oriented downhill in the direction of steepest descent.

Consider adding lumps of sand continually, uniformly and quasistatically to a pile on a table of arbitrary shape. The pile will grow and reach some steady state (Nadai 1950). We wish to characterize such a steady state on the basis of the plausible assumptions on the mechanical behaviour of sand discussed previously.

## 2. Theory

Let ⊂_{2} and define *w* : ^{C}↦ to be the height of the pile as a function of position on ^{C}, where ^{C}=∪∂ is the closure of and ∂ is the boundary of (figure 2). Here, the domain ^{C} represents the top of a horizontal table upon which a pile is to be formed. We now enforce the rules, which are as follows.

*is bounded*. *For convenience, let* *be open also*.

*w*=0 *on* ∂.

*w is C*^{0} *continuous on* ^{C}.

*w is piecewise C*^{1} *continuous on* .

Before we consider the angle of repose, let us define the *directional derivative* .

*where* *the unit circle*.

We may now introduce the angle of repose as a bound on .

*the unit circle*.

At this point, there may be many possible functions that satisfy the requirements for *w*. We wish to choose one that represents a steady state as sand is continually added to the pile. (This *w* will be shown by construction to be unique.) Thus, if sand is added to any point on , we wish *w* to be steep enough that the added material does not accumulate, but flows to another point. Since this flow goes downhill, we must establish a restriction on the minimum directional derivative. We shall see that this restriction is enough to ensure that the added sand flows to the boundary.

.

We may now make some observations about regions where *w* is *C*^{1}. In these regions, the gradient ∇*w* is defined and . Then, rules 2.5 and 2.6 give(2.1)which is known as the *eikonal equation*.

Let *a* be a point in where *w* is *C*^{1} continuous. We define _{a} to be the largest connected open *C*^{1} subset of that contains *a*. We shall refer to _{a} as the *smooth patch* containing *a*. Then, we construct curve *Γ* in _{a}, such that(2.2)We shall refer to any curve that obeys property (ii) as a *streamline*. Now consider a second point *b* on *Γ*. Without loss of generality, suppose that *w*(*b*)>*w*(*a*) (figure 3). Then, we use the properties of *Γ* to note(2.3)where is the length of the streamline between *a* and *b*. Let *Γ*^{*} be any other smooth curve in connecting *a* and *b*. Then, we have(2.4)Comparison of (2.3) and (2.4) shows that for any curve *Γ*^{*} connecting points *a* and *b*. Thus, we conclude that streamlines are straight lines in _{a}.

A straight line in _{a} (an open subset of _{2}) must intersect ∂, the boundary of _{a}, in at least two places. Nevertheless, once the line intersects the boundary of _{a}, the above arguments no longer hold and we cannot continue the construction of the streamline. Thus, the streamlines are straight line segments that intersect ∂_{a} in exactly two places. One of these will have the gradient ∇*w* pointing towards ∂_{a}. We shall call this point the *uphill intersection* and denote it by *x*_{u}. Similarly, the *downhill intersection, x*_{d}, has ∇*w* pointing away from the boundary. By construction, any point in ∂_{a} must be a point in either ∂ or where *w* is not *C*^{1}.

Let us now turn our attention to the points in where *w* is not *C*^{1}. Since *w* is piecewise *C*^{1} on , these points may be either (1) isolated, (2) part of a curve in the *xy* plane of such points, or (3) the intersection of two or more curves. Let where *w* is not *C*^{1} and consider each case separately.

Case (i): ** x** is isolated. Then, for any the unit circle, is defined for sufficiently small

*h*, and(2.5)Now, from rules 2.5 and 2.6 , so that for some unit vector ,(2.6)where

*θ*is the angle between and the nearby gradient . Thus, there is at least one direction where the nearby gradient points towards

**. Then,**

*x***is an uphill intersection of some streamline(s). From here it is straightforward to show that**

*x**w*is locally conical, with a sharp peak at

**(figure 4).**

*x*Case 2: ** x** is part of a curve

*S*∈, such that

*w*is

*C*

^{1}on either side of

*S*, but not on

*S*. Let be a unit vector that is ‘tangent’ to

*S*at

**. In the case that**

*x**S*is not smooth or that

**is at the end of**

*x**S*, we simply construct a vector to be in the direction of a one-sided derivative of

*S*. Then, we may use the continuity of

*w*to obtain(2.7)where ∇

*w*

_{1}and ∇

*w*

_{2}are the nearby gradients on either side of

*S*at

**. Thus, if**

*x**θ*

_{1}and

*θ*

_{2}are the angles between each nearby gradient and , then(2.8)Now if

*θ*

_{1}=

*θ*

_{2}, then there is no discontinuity in slope and the curve is

*C*

^{1}after all, so we conclude that

*θ*

_{1}=−

*θ*

_{2}and

*θ*

_{1},

*θ*

_{2}≠0. The two nearby gradients therefore meet the curve at equal and opposite non-zero angles (figure 5). Both gradients could either (i) point away from the curve, forming a sharp valley at

*S*or (ii) point towards the curve, forming a sharp ridge at

*S*. The sharp valley is prohibited by rule 2.6, since the minimum slope is along the valley, which is greater than

*k*by equation (2.7). Therefore,

*S*must be a sharp ridge. Note that since the nearby gradients on either side of the ridge point towards

*S*,

**represents the uphill intersections of at least two streamlines. There may be more than two streamlines whose uphill intersections are at**

*x***if it is the end of curve**

*x**S*.

Case 3: ** x** is the intersection of curves where

*w*is not

*C*

^{1}. Then

**is simply the intersection of sharp ridges. Each of the intersecting ridges must have nearby gradients pointing towards the ridge, so**

*x***is the uphill intersection of some streamlines, but cannot be a downhill intersection.**

*x*Let us summarize what we have concluded so far. We have shown that streamlines are straight lines within the patches where *w* is *C*^{1}. These streamlines may have uphill intersections at sharp peaks or ridges inside , but the downhill intersection cannot be a non-smooth place within . Thus, the downhill intersection of any streamline must be in ∂. Since the slope of *w* along streamlines is a constant *k*>0, we cannot have both uphill and downhill intersections fall on ∂, so any uphill intersection *must* be at a sharp peak or ridge.

Now we may develop a convenient way of determining the value of *w* anywhere in . Let *b*∈. Let *Γ* be any streamline that contains *b* and let *a* be the downhill intersection of *Γ*. If *w* is *C*^{1} at *b*, there will be exactly one streamline; otherwise, there will be more than one. Then *a* must be in ∂, and by rule 2.2, *w*(*a*)=0. Since *w* is *C*^{0}, equation (2.4) still holds, so that(2.9)Thus, for any ** x**,

*w*(

**) is simply the length of the streamline between**

*x***and the boundary. Note that equation (2.9) implies that for points on multiple streamlines (ridges and peaks), each streamline must have the same length. Now we just need to find a streamline that connects each point in to the boundary. We know that each streamline is a straight line segment from the boundary to an interior point, but we need to be able to choose the streamline(s) from many possible line segments.**

*x*Let *Γ*^{*} be any line segment in ^{C} from the boundary of to *b*. Let *a*^{*} be the intersection of *Γ*^{*} with the boundary. Then an argument similar to the argument that led to (2.4) may be used to show that(2.10)Comparison of (2.9) and (2.10) shows that . Therefore, the streamline through any point is simply the shortest line segment connecting ** x** to the boundary. Furthermore, equation (2.9) now gives(2.11)where is the shortest distance between

**and the boundary. Equation (2.11) provides a useful way of constructing the (unique) function**

*x**w*for any given domain (Alouges & DeSimone 1999).

## 3. Example

For simplicity, we set *k*=1; it turns out (§4) that this assumption is not restrictive. We proceed as follows. (i) We discretize the boundary of as a list of points *a*_{k}. (ii) We cover with a two-dimensional uniform grid *x*_{ij} of spacing comparable with the spacing between boundary points. (iii) For each point in the grid, *x*_{ij}, we find the closest point on the boundary, . (iv) For *x*_{ij}, we compute , where . (v) We add surface shading on the basis of the direction of *r*_{ij}, since is approximately anti-parallel to ∇*w* and parallel to the horizontal component of the surface normal . Figure 6 shows the surface of a pile computed in this form.

## 4. Surface properties

The surface produced by equation (2.11) has a number of interesting properties. We present here a partial list of these properties, along with some proofs or explanations of how these properties arise.

*Of all C*^{0} *functions that vanish on* ∂ *and have slope bounded by k* (*rules* *2.2, 2.3* *and* *2.5*), *w is maximal*.

Let *z* : ^{C}↦ be a *C*^{0} function, such that *z*(∂)=0. Suppose that for some . Let *Γ* be the streamline through ** x**. Then, the average slope of

*z*along

*Γ*is greater than

*k*. Since

*z*is

*C*

^{0}, the local slope must be greater than

*k*somewhere along

*Γ*. Thus, the surface

*w*represents the maximum height that any pile can take on the given domain.

*Smooth patches are ‘developable’ surfaces; that is, any point in a smooth patch has a direction with zero curvature*.

For all practical purposes, a developable surface is one that may be constructed out of paper. A sphere, for example, is not developable because a point on the surface has non-zero curvature in every direction. The pile surface is developable because there is zero curvature along the gradient direction.

Area(*w*())/Area .

This result follows easily from equation (2.1).

*All non-trivial domains give at least one point where w is not C*^{1}.

The gradient ∇*w* always points away from the boundary, so the maximum value of *w* must be in the interior. Also, since there are no *C*^{1} points with zero slope, the maximum must occur at a sharp ridge or peak.

*The construction of w is independent of scale*.

The construction was developed in such a way that expanding or contracting the domain by some factor just causes the height *w* to change by the same factor. In addition, *w* may be scaled independently of the domain by simply changing the critical slope *k* (equation (2.11)). This independent scaling justifies the assumption *k*=1 used for calculations. In practice, both real granular materials and numerical implementation of the construction have a significant length-scale—grain size and point spacing, respectively. For a theoretical pile to match a real pile, both length-scales must be smaller than any features of interest in the domain or pile surface.

⊂, *a*∈⇒(*a*)≤(*a*), *where* *and* *are the pile heights on the domain* *and the subdomain* , *respectively*.

This result comes easily from the definitions of subset and *d*_{min}: the boundary of *C* cannot be farther from *a* than the boundary of ⊃.

*w*=max(*cones with slope k and base in* ^{C}).

This alternative construction is a result of the previous observation. For any , a circle of radius fits in ^{C}. Then, a cone with slope *k* on this circular base will be equal to *w* at ** x** and less than or equal to

*w*everywhere else.

*w*=min(*inverted cones with slope k and vertex on* ∂).

This third construction is just one way to visualize (equation (2.11)). The height on the surface of these cones is just *k*.*d*, where *d* is the distance between and the cone vertex (on ∂).

*Streamlines are perpendicular to* ∂ *if* ∂ *is smooth where they intersect*.

If the boundary is smooth, so is *w* (locally), and level curves such as the boundary must be perpendicular to ∇*w*.

*All valleys are smooth, straight and without tributaries, and have slope equal to k*.

Sharp valleys were already shown to be prohibited. Any ‘valley’ in a theoretical pile is simply the collection of straight streamlines that intersect a concave portion of ∂. Since the downhill intersection of any streamline must be in ∂, all valleys lead directly to ∂, and there are no tributary valleys.

*All sharp ridges have slope less than k*.

This property was already shown, but is repeated here for emphasis.

## 5. Analogy to other physical systems

The sandpile construction is relevant to other physical systems, e.g. elastic and ferromagnetic thin films. Details of the results summarized here may be found in the original papers cited in the references.

### (a) Buckling of compressed thin films

Suppose that a thin film is bonded to the flat surface of a thick substrate everywhere except over a domain . If the substrate is subjected to a sufficiently large compressive isotropic strain *ϵ*^{*} parallel to its surface, the film will buckle into a deflected shape *w*(), with *w*≥0. Ortiz & Gioia (1994) have shown that minimization of energy from membrane stretching leads to the eikonal equation |∇*w*|=*k*, with (where *ν* is the Poisson's ratio for the film). The eikonal equation alone, along with the boundary condition *w*(∂)=0, gives infinitely many minimizers of the membrane energy, each the upper envelope of some set of cones whose circular bases cover ^{C}. (Compare with property 4.7, where all possible cones are considered.) Of these possible minimizers of the membrane energy, the physically preferred solution proposed by Ortiz & Gioia (1994) is found by minimizing the bending energy associated with sharp folds (which correspond to ridges in the sandpile model). This strategy gives the same *w* as found by the sandpile construction. More rigorous results concerning this proposed solution may be found in Jin & Kohn (2000) and Jin & Sternberg (2001). Further developments can be found in Gioia *et al.* (2002).

### (b) Magnetization of ferromagnetic thin films

Ferromagnetic materials have a magnetization vector ** m** (magnetic moment per unit volume) of constant magnitude (which is a property of ferromagnetic materials) but variable direction. In the limit of very thin ferromagnetic films, Stoner & Wohlfarth (1948) showed that

**should be in the plane of the film for thin elliptical films. It has been conjectured by Van den Berg (1986) and justified mathematically by Gioia & James (1997) that**

*m***should be in the plane of any thin ferromagnetic film, regardless of shape. This result reduces the problem of minimizing the potential energy to the eikonal equation in the film and**

*m**w*=const. on the film boundary, where

*w*is the magnetic potential. Van den Berg (1986) constructed solutions by the method of characteristics (which correspond to streamlines in the sandpile model) and found that any given film could have many possible minimizers

*w*of the potential energy, each with a different set of

*domain walls*(which correspond to ridges in the sandpile model). These domain walls represent boundaries between magnetic domains or regions of locally aligned magnetization vectors. Domain walls store energy (the

*exchange energy*of micromagnetics) that should also be minimized (DeSimone

*et al.*2000, 2006). If the wall energy is taken to be a constant per unit length, the physically preferred potential is the same

*w*(or, alternatively, −

*w*) as found by the sandpile construction.

## 6. Experiments

In these experiments, we pour fine sand onto a plywood table until the table is ‘full’. We photograph the resulting pile and compare it with the theoretical pile.

A granular material has to be chosen carefully to give good results. Grains that are too fine or too damp will have significant cohesive forces that will affect pile topography. Grains that are too coarse will require a very large-sized experiment to lessen the influence of individual grains. Also, non-uniform grains will tend to cause segregation. After trying a few materials, we select fine mason's sand owing to its availability, relative uniformity and convenient grain size.

Figures 7–9 show the theoretical and experimental piles for domains of three different shapes. The piles show similar ridge patterns, but there is a notable smoothness to the ridges of the experimental piles. This effect may be explained by the impact of falling grains. We take care to add the sand gently, but the grains still reach the surface with a small momentum. Ridges have less material supporting the grains at the surface, so the momentum of incoming grains is more likely to cause flow.

## 7. Conclusions

On the basis of a few plausible assumptions on the mechanical behaviour of sand, we have justified mathematically Nadai's method for constructing the surface of the steady sandpile on a domain of arbitrary shape. We have verified experimentally that the method gives a good approximation of the surface of actual steady sandpiles. Depending on the domain, the surface may be as simple as a cone or it may be a more complex surface with an elaborate network of sharp ridges. The surface is relevant to physical systems other than sandpiles, and we have discussed briefly two examples of such systems. We have shown that the surface has a number of interesting properties, which may be of some practical value for dealing with granular materials. For example, the volume of a stockpile may be estimated, given the shape of the domain. The method may be extended to describe the case where the domain is not planar, such as the steady sandpile in a container with an open top.

## Acknowledgements

We thank Horacio Porta for helping us simplify the derivation of equation (2.4) and James W. Phillips for suggesting improvements to the text. This work was financed in part by the Critical Research Initiatives Program, University of Illinois at Urbana-Champaign.

## Footnotes

- Received September 27, 2006.
- Accepted January 4, 2007.

- © 2007 The Royal Society