## Abstract

Plastic grains are found to form fractal patterns in elastic-hardening plastic materials in two dimensions, made of locally isotropic grains with random fluctuations in plastic limits or elastic/plastic moduli. The spatial assignment of randomness follows a strict-white-noise random field on a square lattice aggregate of square-shaped grains, whereby the flow rule of each grain follows associated plasticity. Square-shaped domains (comprising 256×256 grains) are loaded through either one of three macroscopically uniform boundary conditions admitted by the Hill–Mandel condition. Following an evolution of a set of grains that have become plastic, we find that it is monotonically plane filling with an increasing macroscopic load. The set’s fractal dimension increases from 0 to 2, with the response under kinematic loading being stiffer than that under mixed-orthogonal loading, which, in turn, is stiffer than the traction controlled one. All these responses display smooth transitions but, as the randomness decreases to zero, they turn into the sharp response of an idealized homogeneous material. The randomness in yield limits has a stronger effect than that in elastic/plastic moduli. On the practical side, the curves of fractal dimension versus applied stress—which indeed display a universal character for a range of different materials—offer a simple method of assessing the inelastic state of the material. A qualitative explanation of the morphogenesis of fractal patterns is given from the standpoint of a correlated percolation on a Markov field on a graph network of grains.

## 1. Introduction

Over the past few decades, various materials have been observed to display fractal patterns (Mandelbrot 1982; Feder 2007). Mechanisms of morphogenesis of fractals for various physical systems were proposed: liquid streams meandering over land landscapes, diffusion-limited aggregation, percolation of microcracks in elastic–brittle solids, formation of shear bands in rocks or formation of plastic ridges in Arctic ice fields (e.g. Sahimi & Goddard 1986; Ostoja-Starzewski 1990; Poliakov & Herrmann 1994; Poliakov *et al*. 1994; Saouma & Barton 1994; Zaiser *et al*. 1999; Sahimi 2003; Shaniavski & Artamonov 2004; Sornette 2004).

In a recent study (Li & Ostoja-Starzewski in press), we found that fractal patterns of plasticized grains form at elastic–plastic transitions in two models of linear elastic/perfectly plastic random heterogeneous materials: (i) a composite made of locally isotropic grains with weak random fluctuations in elastic moduli and/or yield limits and (ii) a polycrystal made of randomly oriented anisotropic grains. By ‘plasticized’ we mean grains that have gone from elastic into plastic response. In each case, the spatial assignment of material randomness has been a non-fractal, strict-white-noise field on a 256×256 square lattice of homogeneous, square-shaped grains, with the flow rule in each grain following associated plasticity. These lattices were subjected to pure shear loading increasing through either one of three macroscopically uniform boundary conditions (BCs) (kinematic, mixed-orthogonal or traction), admitted by the Hill–Mandel condition. Upon following the evolution of a set of grains that become plastic, we found it to have a fractal dimension increasing from 0 towards 2 as the material transitions from elastic to perfectly plastic. While the grains possessed sharp elastic–plastic stress–strain curves, the overall stress–strain responses are smooth and asymptote towards perfectly plastic flows; these responses and the fractal-dimension–strain curves were almost identical for three different loadings. The randomness in elastic moduli alone was sufficient to generate fractal patterns at the transition, but had a weaker effect than the randomness in yield limits. In the model with isotropic grains, as the random fluctuations tended to zero (i.e. the composite becomes a homogeneous body), a sharp elastic–plastic transition was recovered.

In the present case, we turn to more realistic model materials with elastic-hardening plastic, isotropic grains having random yield limits or elastic/plastic moduli. In particular, we pose the following questions: (i) does the elastic–plastic transition occur as a fractal, plane-filling process of plastic zones under increasing, macroscopically uniform applied loading? (ii) Is the fractal character of plastic zones robust under significant changes of material properties? (iii) Can the fractal character of the set of plastic grains be explained through a stochastic model?

In §2, a white-noise random microstructure is set up, the key property of the model being that the material given *a priori* is non-fractal. Then, in §3, we determine fractal dimensions of sets of evolving plastic grains when the body is subjected to monotonic loading(s) consistent with the Hill–Mandel condition; several models and kinds of randomness are examined. Section 4 presents a Markov-field model of stochastic and fractal evolution of plastic grains, where the Markov property is dictated by the nearest-neighbour interactions between the contiguous grains.

## 2. Model formulation

A random heterogeneous material is defined as a set **B**={*B*(*ω*);*ω*∈*Ω*} of deterministic media *B*(*ω*), where *ω* indicates a specific realization and *Ω* is an underlying sample space (e.g. Ostoja-Starzewski 2008). The material parameters of any microstructure, such as the elasticity tensor or the yield tensor, jointly form a random field ** Θ** that is required to be mean ergodic on (very) large scales, that is
2.1
Here, the overbar indicates the volume average and 〈 〉 means the ensemble average.

*P*(

*ω*) is the probability measure assigned to the ensemble {

*G*(

*ω*,

*x*);

*ω*∈

*Ω*,

*x*∈

*V*} and its

**-algebra.**

*σ*We define a homogenized response as that in which there is equivalence between energetically () and mechanically () defined effective responses, . This is the well-known Hill–Mandel condition in linear elastic materials, leading to three types of uniform BCs:

(i) kinematic (displacement) BC (with applied constant strain

*ε*^{0}): 2.2(ii) traction (static) BC (with applied constant stress

*σ*^{0}): 2.3(iii) mixed-orthogonal (or displacement–traction) BC: 2.4

where **u** is the displacement vector and **t** is the traction vector on the specimen boundary ∂*B*_{δ}. The above boundary conditions can be generalized to elastic–plastic materials in an incremental setting (Hazanov 1998; Ostoja-Starzewski 2008). The microstructures in our study are made of perfectly bonded, homogeneous, isotropic grains of linear elastic/linear hardening plastic (*J*_{2}) type with an associated flow rule (Simo & Hughes 1998). Specifically, the constitutive response of any grain (i.e. a piecewise-constant region in a deterministic microstructure *B*(*ω*)) is described by
2.5a
2.5b
2.5c
where primes indicate deviatoric tensor parts, *E* is the elastic (Young’s) modulus, *ν* is Poisson’s ratio and *f* is the yield function. Equation (2.5b) indicates elastic deformation or plastic unloading and equation (2.5c) refers to the plastic stage. We consider the von Mises–Huber isotropic yield criterion and *f* is defined as
2.6
Here, *c* is the yield constant. In the one-dimensional case, the constitutive relations (2.5*a*–*c*) reduce to
2.7a
and
2.7b
where *σ*_{s} is the yield stress in uniaxial tension (), *ε*_{p} is the plastic strain and *E*_{p} is the plastic modulus. We can see that equations (2.7a,*b*) constitute a simple nonlinear model (piecewise-linear) characterizing the elastoplastic response. Figure 1 shows such a homogeneous body transitioning instantaneously from elastic (white) to plastic (black) state, accompanied by a kink in the stress–strain curve.

## 3. Simulations of fractal patterns of plastic grains

### (a) Evolutions of fractal patterns

We consider a two-dimensional plane-strain problem and apply shear loading through one of the three types of uniform BCs consistent with equations (2.2)–(2.4): 3.1

A computational mechanics study is conducted with a finite-element method commercial software (Abaqus 2008). We take a sufficiently large domain comprising 256×256 square-shaped grains. Each individual grain is homogeneous and isotropic, with its *E* and *E*_{p} being constant and *σ*_{s} being a field of independent identically distributed (i.i.d.) uniform random variables scattered up to ±2.5 per cent about the mean. Other kinds of randomness are studied in §3*c*. The material constants are taken from ABAQUS 6.8 Benchmarks 3.2.1 case 1,

Figure 2*a*–*d* shows elastic–plastic transition fields at different deformation stages under traction BCs. We follow here the binary format of figure 1 in the sense that elastic (plastic) grains are white (black). The plastic grains form plastic regions of various shapes and sizes, and we estimate the fractal dimension *D* of that entire plastic grain set by using a ‘box-counting method’ (Perrier *et al*. 2006),
3.2
where Nr denotes the number of boxes of size *r* needed to cover the object. Results of box counts for figure 2*a*–*d* are shown in figure 3*a*–*d*, respectively. With the correlation coefficients of log(Nr) versus log(*r*) close to 1 for all black–white patterns, we conclude that the elastic–plastic transition patterns are fractal. (Our computer code for the box-counting method was verified on the Sierpiński carpet: the theoretical value is and our estimation gives *D*=1.8911, with the correlation coefficient being 0.9999.) Note the orthotropic character of the plastic grain set, which reflects the propensity of shear bands and slip lines to form approximately at ±45^{°}.

Figure 4*a*,*b* shows response curves under these three BCs in terms of the averaged stress versus strain and the fractal dimension (*D*) versus strain, respectively. The responses of single-grain homogeneous phases are also given for reference. In both figures, the curves overlap, demonstrating that the (256×256) domain is the representative volume element (RVE). Here, we assign each grain with one finite element for convenience. To verify whether such a meshing is sufficient, we conduct numerical simulations using different element types (linear or quadric) and illustrate the results in figure 5*a*,*b*. Clearly, the curves are nearly identical. Thus, the 256×256 domain size is chosen so as to ensure the RVE-level response and computational accuracy, while having an acceptable spatial resolution of dependent fields to ensure a reliable assessment of fractal dimensions.

As expected from the hierarchies of scale-dependent bounds (Ostoja-Starzewski 2005, 2008), the response under mixed-orthogonal loading is always bounded from above and below by the kinematic and traction loadings, respectively. Here, we recall the hierarchies of bounds for elastic-hardening plastic composites (Jiang *et al*. 2001; Ostoja-Starzewski 2008), such as
3.3
In inequalities (3.3) and are the apparent tangent stiffness and compliance moduli, respectively, with the superscript d (or t) indicating the displacement (or traction) BC. Another type of hierarchy that applies is in terms of energies, as expressed by eqn (3.5) in Ostoja-Starzewski & Castro (2003).

### (b) Discussion of hardening effects

In order to investigate the influence of hardening in the transition, we conduct comparison studies among different materials with various plastic moduli and/or yield strain (*σ*_{s}/*E*). Their material parameters are listed in table 1, with Poisson’s ratio being 0.3 in all the cases. Here, we note that materials 1 and 2 are taken from Abaqus 6.8 Benchmarks 3.2.1 cases 1 and 2, respectively, while materials 3a and 3b are designated to be different in yield strain from materials 1 and 2; their parameters *E*,*σ*_{s} are those of 316 steel in ABAQUS 6.8 Example Problems 1.1.8. The series A materials (2a and 3a) have the same ratio of plastic-to-elastic moduli, *E*_{p}/*E*, as material 1, while series B materials (2b and 3b) are the same in their plastic moduli *E*_{p}. The compared responses are considered in terms of dimensionless quantities such as normalized stress or strain (rescaled by yield stress or yield strain) accordingly.

We begin the discussion by comparing responses of materials 1 and 2. As seen from table 1, these materials have the same yield strain but very different hardening parameters *E*_{p}. Note that the hardening constant of material 2 is much less than that of material 1, and we see from figure 6*a* that the tendency of approaching a homogeneous response in material 2 is slower than that in material 1. Accordingly, the fractal dimension grows slower in material 2 as observed in figure 6*b*, which means that the elastic–plastic transition develops more rapidly in materials with stronger hardening effects.

Note that materials 1 and 2 also differ in terms of elastic moduli. To see the hardening effects more clearly, we perform numerical simulations on two hypothetical materials 2a and 2b—the former with the same *E*_{p}/*E* of material 1, while, for the latter, the same *E*_{p}. The results are given in figure 7 for comparison. We observe that responses of materials 1 and 2a are nearly identical in the curves of both, the normalized stress–strain and the fractal dimension–strain. For material 2b with the same *E*_{p} of material 1 but a lower *E*_{p}/*E*, the tendency of approaching the homogeneous response appears slower, and the same conclusions can be drawn regarding the fractal-dimension curves.

While the study of materials 1 and 2 (figure 6*a*,*b*) has been restricted to cases of the same yield strain, we now turn to materials with different yield strains. This is done by comparisons of material 1 and two other hypothetical materials 3a and 3b—the same configurations of *E*_{p}/*E* or *E*_{p} but different yield strains. The results are illustrated in figure 8. We arrive at the same conclusions as for materials 2a and 2b. Based on these observations, we conclude that the elastic–plastic transition in random materials is characterized by their fractal-dimension–normalized plastic-strain responses. Specifically, the same response curves in fractal dimension–normalized plastic strain results in the same normalized stress–strain relations. The hardening effect is determined thoroughly by the ratio of plastic-to-elastic moduli, *E*_{p}/*E*. The same *E*_{p}/*E* leads to the same fractal transition patterns, and the larger the *E*_{p}/*E*, the faster the increase in fractal dimensions, i.e. the hardening facilitates the elastic–plastic transition in random materials.

While our study is restricted to linear hardening materials, a simple extension to nonlinear hardening models can also be drawn here—for those cases, the hardening parameter *E*_{p} is not a constant but some function with regard to the state of plastic strain, and the ratio *E*_{p}/*E* can be specified in each incremental step during the transition. We note that, in conventional stress–strain calibrations, the trends to approach homogeneous response curves are not easy to discern among different materials. On the other hand, the fractal dimension always increases from 0 to 2 during the transition, thus providing an optimal parameter to assess the transition process.

### (c) Other types of randomness

Now, we examine the influence of material-parameter randomness. This is again studied through comparisons of several cases. As suggested in §3*b*, we consider the response curves in terms of fractal dimension versus normalized plastic strain. First, the sensitivity of transition processes to the material’s model randomness is investigated in two scenarios.

—

*Scenario A*. Scalar random field of the yield limit, with three types of randomness,A1: yield limit is a uniform random variable up to ±2.5 per cent about the mean,

A2: yield limit is a uniform random variable up to ±0.5 per cent about the mean, and

A3: deterministic case: no randomness in the yield limit.

—

*Scenario B*. Random field of the yield limit or elastic moduli, with two cases,B1 = A1 and

B2: elastic modulus is a uniform random variable up to ±2.5 per cent about the mean.

Results for A1–A3 and B1–B2 are shown in figures 9 and 10, respectively. From figure 9, one can conclude that different random variables in the model configuration lead to different transition processes; overall, a lower randomness results in a faster elastic–plastic transition. Next, in figure 10 we observe the randomness in yield limits to have a stronger effect than that in elastic modulus.

Note that the hardening parameter (plastic modulus, *E*_{p}) can also be randomly perturbed. In fact, it is more realistic to combine the randomness in elastic and plastic moduli in applications of composites of various material constitutes. We conduct this study through a comparison of three cases.

—

*Case 1*. B2, the plastic modulus is fixed.—

*Case 2*. Elastic modulus is a uniform random variable up to ±2.5 per cent about the mean, while the plastic modulus is a corresponding linearly dependent random variable:*E*_{p}/*E*is fixed.—

*Case 3*. Elastic and plastic moduli are i.i.d. uniform random variables up to ±2.5 per cent about the mean.

Figure 11 illustrates the results for these cases. We find that the three curves are practically overlapping, which means that the randomness in the plastic modulus has a very small effect in the transition.

### (d) Further discussions

In §3*a*–*c*, we demonstrated that the set of all plastic grains at the elastic–plastic transition is an evolving fractal, and we studied hardening effects and material model randomness appealing to the fractal dimension. While the fractal dimension versus normalized plastic strain curves prove effective in comparing transition patterns among different materials, we note that the fractal dimension always reaches 2 at the end of the transition, which hints at a universal property for all elastic–plastic materials. To demonstrate this clearly, we plot the curves of fractal dimension versus normalized von Mises stress for materials 1, 2b and 3b in figure 12*a*. We choose these three materials as in §3*b*, and their fractal-dimension–normalized plastic strain curves are different. Interestingly, from figure 12*a*, it appears that the curves are nearly identical for all these materials.

One may argue that a simple variable such as the volume fraction of plastic regions could also assess the transition, since it reaches 1 in the fully developed stage. Accordingly, we plot the results of plastic volume fraction versus normalized von Mises stress in figure 12*b*. Although the responses of materials 2b and 3b nearly overlap (in effect the *E*_{p}/*E* ratios are very close), they differ from material 1. This leads us to conclude that the relation between *D* and the normalized von Mises stress is a universal property (or at least varies very little) for all isotropic elastic-hardening plastic materials. Since the fractal dimension can readily be estimated from, say, image analysis, this provides an effective approach to infer the stress in the material at the transition.

## 4. Why are fractals observed in elastic–plastic materials?

The fractal character of evolving sets of plastic grains observed in the computational mechanics simulations reported above may be explained by a reference to fractals on Markov random fields (MRFs). To see this, we first introduce a binary random variable *S* describing the state of any grain as
4.1
where *e* means an elastic state, *p* means a plastic state and *f*_{p} refers to the yield function of grain *p*. Next, consider grain centres of our polycrystal as a Cartesian lattice of spacing *a* in , that is
4.2
where *m*_{1} and *m*_{2} are integers ranging from 1 through to *N* (=256 in our simulations above). Given that (i) the properties (yield stresses or elastic/plastic moduli) of each grain are random and (ii) the state of each grain is a result of all the interactions in the entire system of all grains, the state *S* on *L*_{a} is a random field,
4.3
In other words, for any *ω*∈*Ω* (a particular realization of the entire material system) and any location **x** on the lattice, the state *S* is *s* (i.e. either *e* or *p*). If we do not specify **x**, then *s* stands for a realization *S*(*ω*)=*s* on *L*_{a}, where some grains are elastic and others plastic. Note that the body is subjected to a macroscopic load: either *ε*^{0} or *σ*^{0}, according to BCs (2.2) or (2.3).

*Markov property*. Recognize that, given the macroscopic load such as *ε*^{0}, the conditional probability of a grain at **x** being plastic at any macroscopic load level depends not on the state of all other grains *L*_{a}−{**x**}, but only on the state of its *nearest neighbourhood**N*_{x},
4.4
This relation defines *S* of equation (4.4) as an MRF. Given the square lattice topology of our composites, *N*_{x} comprises four neighbouring grains,
4.5
whereby we consider the effect of four grains acting through the vertices onto **x** as being much weaker than those that contact through the sides. Since this accounts for neighbours just one step away, *S* is called a 1-Markov field. It is understood, but not written explicitly in equation (4.4), that the conditional probabilities on both sides of that equation depend on the macroscopic applied loading, i.e. *ε*^{0} or *σ*^{0} or some combination thereof, applied, respectively, through equations (2.2)–(2.4).

Note that the nearest neighbourhood may be replaced by a set of neighbours up to two steps away (in the Manhattan metric), whereby *S* would be a 2-Markov field. This would lead to a richer model, but we do not need to pursue it here because, already, the 1-Markov field possesses the salient features of a correlated percolation.

The formulation above is analogous to that of an MRF for an Ising magnet on a square lattice, where the state (spin up or down) of each site is a function of the spins at four neighbouring sites and of the overall temperature *T* (rather than that of a mechanical state)
4.6
The temperature also appears when one writes a Gibbs specification of the random field,
4.7
This is called a *Gibbs random field* (GRF), with *Z* being a partition function ensuring the probability measure is normalized to 1,
4.8

Going back to the elastoplastic composite, instead of equation (4.7), we can also write a Gibbs specification,
4.9
Now, the equations (4.7) and (4.9) also betray the second reason for the formulation above to be analogous to the MRF for an Ising magnet. Namely, the internal energy of the latter is *U*(*ω*,*T*), while the internal energy of the MRF of an elastic–plastic composite is *U*(*ω*,*ε*^{0}). This is consistent with a continuum thermomechanics picture where the temperature is a control parameter for a thermal problem, while strain is a control parameter for a mechanical problem.

It is well known that every MRF is equivalent to a GRF, and vice versa (Preston 1974). However, the key issue that arises is whether the description in terms of conditional probabilities is equivalent to the description in terms of absolute probabilities . That is, if we specify an MRF in terms of local interactions, do we also specify its probability measure *Π* (and, therefore, its state *ω*) in a unique way? The answer to this question depends on whether the interactions are weak or strong. When they are weak in the sense that the state *s*(**x**) depends weakly on the states *s*(*N*_{x}) of the neighbours, then there is a unique correspondence between and . On the other hand, when they are strong (in the sense of strong such dependence), for a given specification of interactions , there is more than one probability measure . The critical point where this non-uniqueness occurs is characterized by fractal patterns.

In the case of the Ising model, this critical point is the Curie point *T*_{C} on the temperature scale, below which we have a ferromagnet and above which there is a spatial disorder of spins so that no single dominant (and hence macroscopic) spin emerges. A wide range of binary patterns—i.e. white (W) versus black (B) vertices—have been analysed for the entire range of control parameters by Hammersley & Mazzarino (1983). The control parameters are *α* (the likelihood of any single vertex to coalign with the external magnetic field) and *β* (the strength of interaction), so that a canonical form of the internal energy reads
4.10
where |*V*_{B}| is the number of black vertices and |*V*_{BW}| is the number of pairs having one black and one white vertex.

Note that the first term in equation (4.10) here is responsible for a Bernoulli type (i.e. uncorrelated) percolation on the lattice. While this percolation—i.e. retaining only the first term in equation (4.10)—is well known to also exhibit fractal patterns, the mechanics and physics considerations suggest that the second term is non-zero, so the elastic–plastic transition is characterized by a correlated percolation. In fact, for a weakly random microstructure, the plastic state is likely to ‘spill over’ to a neighbouring elastic grain. On the other hand, in a strongly random microstructure, plasticity tends to go via weak grains around grains with high yield limits (Jeulin *et al.* 2008). Thus, in some cases, the interactions may be attractive, and in others repulsive. Note that (i) the correlated percolation involves an interplay of both terms in equation (4.10) and (ii) the problem being tensorial in nature and orthotropic, it is more complex than what happens in the Ising model, which is only scalar.

In the early 1980s, fractals were hardly known, and this is probably why Hammersley & Mazzarino (1983) did not estimate fractal dimensions from their computer simulations, although fractal patterns are clearly seen in their figures at (and around) *T*_{C}. The task of generating fractal patterns via MRF models and computing their fractal dimensions was accomplished by Onural (1991) and followed by others in the field of image analysis (e.g. Ghozi 2001).

Returning back to our elastic–plastic transition in a random composite, we recapitulate that

(i) there are elastic (

*e*) and plastic (*p*) vertices in analogy to W and B vertices in the Ising model;(ii) the increasing applied loading

*ε*^{0}tends to cause the*e*→*p*transition at any single vertex, while the local conditioning is attractive in the sense that*p*states on*N*_{x}tend to make*S*(*x*)=*p*(with the same cause–effect relation holding for*e*); and(iii)

*S*is the MRF so that the evolution of the entire*V*set from a predominantly*e*state to a predominantly*p*state exhibits fractal patterns.

Since the responses under equations (2.2) and (2.3) loadings have been shown to be almost the same (i.e. the RVE level), the above arguments could be restated with *ε*^{0} replaced by *σ*^{0}. While this section provides only a qualitative explanation of the morphogenesis of fractal patterns at elastic–plastic transitions, a quantitative determination of conditional probabilities of the MRF is outside the present study. An MRF has also been used to model damage phenomena in polycrystals, albeit on the network of grain boundaries instead of grains, and was introduced in Ostoja-Starzewski (1989*a*,*b*).

## 5. Conclusions

We study evolving sets of plastic grains at elastic–plastic transitions in several random two-dimensional model materials. In particular, these are composites of linear elastic/hardening plastic, isotropic grains having fluctuations in yield limits or elastic/plastic moduli, which are taken as strict-sense white-noise (non-fractal) random fields. It turns out that, in every model material subjected to increasing uniform loading (kinematic, traction or mixed), the sets of evolving plastic grains are fractal, with the fractal dimension growing from 0 in the globally elastic state to 2 in a globally plastic state. Several models and kinds of randomness are examined, always displaying the same essential features.

As mentioned in §1, the formation of fractal geometries was also studied in elastic–brittle materials in the 1980s and 1990s. It was found there that a fractal crack was approached only asymptotically, whereas at the elastic/plastic transitions, the set of plastic grains is a continually evolving and ultimately plane-filling fractal structure. That is, the fractal dimension of that set changes—indeed in both the elastic/hardening plastic studied here and the elastic/perfectly plastic materials studied in Li & Ostoja-Starzewski (in press), as well as in several variations thereof. Also, as the random field noise vanishes, in both cases, the transition from elastic to plastic occurs instantaneously (i.e. *D* goes from 0 to 2 discontinuously), whereby the smooth response curve of the random composite turns into a stress–strain curve with a kink. This suggests that elastic–plastic transitions in random materials are representations of one universality class. This research has a very practical application: the curves of fractal dimension versus applied stress—which indeed display a universal character for a range of different materials—offer a simple, image-based method of assessing the inelastic state of the material. Finally, we present a Markov-field model of stochastic and fractal evolution of plastic grains, where the Markov property is dictated by nearest-neighbour interactions between the contiguous grains.

## Acknowledgements

We benefited from comments by anonymous reviewers. This work was made possible by the NCSA at University of Illinois and the NSF support under the grant CMMI-0833070.

## Footnotes

- Received July 9, 2009.
- Accepted October 7, 2009.

- © 2009 The Royal Society