Fast physically accurate rendering of multimodal signaturesof distributed fracture in heterogeneous materials
Yon Visell, Member, IEEE

Fig. 1. Fast, interactive rendering of fracture in heterogeneous materials.
Abstract— This paper proposes a fast, physically accurate method for synthesizing multimodal, acoustic and haptic, signaturesof distributed fracture in quasi-brittle heterogeneous materials, such as wood, granular media, or other fiber composites. Fracture
processes in these materials are challenging to simulate with existing methods, due to the prevalence of large numbers of disordered,quasi-random spatial degrees of freedom, representing the complex physical state of a sample over the geometric volume of interest.
Here, I develop an algorithm for simulating such processes, building on a class of statistical lattice models of fracture that have been
widely investigated in the physics literature. This algorithm is enabled through a recently published mathematical construction basedon the inverse transform method of random number sampling. It yields a purely time domain stochastic jump process representing
stress fluctuations in the medium. The latter can be readily extended by a mean field approximation that captures the averaged
constitutive (stress-strain) behavior of the material. Numerical simulations and interactive examples demonstrate the ability of thesealgorithms to generate physically plausible acoustic and haptic signatures of fracture in complex, natural materials interactively at
audio sampling rates.
Index Terms—Physical simulation, multimodal rendering, virtual reality

1

I NTRODUCTION

It has been a longstanding goal in virtual reality and computer animation to simulate the interactive, multimodal behavior of natural objects
and surroundings in plausibly realistic ways. However, the mechanics of many common interactions and materials are highly complex.
For example, when subjected to sufficiently large stresses, a diverse
variety of natural media, ranging from wood to gravel, give rise toquasi-random, fracture-like processes yielding fast and unpredictable
force fluctuations [31, 10, 41, 6, 43, 28, 42, 30]. Some of these include
fracture in natural fiber composites [14, 2, 1] and stress fluctuations
during shear sliding on granular media [4, 8, 34, 33]. Such processes
have attracted interest in scientific and engineering research becauseof the unusual variety of phenomena they exhibit, and because of theirrelevance to many natural materials and physical interactions [14, 2].
A full account of the physics of these processes remains elusive, due tothe microscopically disordered nature of the materials and their large
number of constituent elements. Nonetheless, some of the most salient
aspects of their behavior have been revealed in statistical physics analyses undertaken during the last few decades [2, 14, 31].
A basic problem is to predict the macroscopic patterns of failurefrom micromechanical properties, and to simulate their behavior without explicitly integrating the full micromechanical state, which is of-• Author is with Drexel University. E-mail: yon.visell@drexel.edu.
Manuscript received 18 September 2014; accepted 10 January 2015; posted
online 23 March 2014; mailed on 1 May 2014.
For information on obtaining reprints of this article, please sende-mail to: tvcg@computer.org.

ten computationally prohibitive. Stress distributions in granular media
also share many common features with fiber composites. Remarkably, both material categories can be considered to bear stress predominantly along arrays of brittle, one-dimensional structures withrandomly distributed strengths. In prior literature, as in the present
contribution, this has been exploited to facilitate both analysis andmodeling.
Rendering brittle fracture in heterogeneous materials
Methods for the computational rendering and animation of fracturelike processes have advanced in recent decades. The notion of simulating material deformation, fracture, and related phenomena was introduced by Terzopoulos [47]. Numerous subsequent studies have introduced different methods of fracture rendering [20, 35, 36, 39, 49, 32,
45, 49, 12, 27, 29]. A review would be beyond the scope of the present
contribution. Fewer studies have investigated rendering fracture in microscopically heterogeneous materials, such as wood. To date, thosethat have done so are based on physically inaccurate heuristics [38], or
have been otherwise impractical for rendering damage in macroscopically large samples, due to the high density of constituent elements
involved [46]. In contrast, the methods we introduce here are strongly
grounded in the relevant physics, and can efficiently render fracture inlarge ensembles of fibers in real time.
Particle-based rendering methods have been used to animate granular media [5, 50, 25]. However, known algorithms are unable to capture the characteristic spatial and temporal stress fluctuations that are
observed in real granular materials [14]. Multimodal (haptic and/or
auditory) rendering of such phenomena is even more prohibitive, due
to the high sampling rates that are required for perceptual fidelity, forFig. 2. Fiber bundle models, the foundation of the present work, havebeen successfully used to model brittle fracture in materials that effectively bear stress along arrays of one-dimensional structures loaded inparallel. Examples include fiber composites such as wood (left), andgranular media, such as sand or gravel (right).

stable, artifact-free rendering [26], and due to the rapid temporal fluctuations involved.
The mechanics of quasi-brittle, heterogeneous materials is complicated by the predominance of highly disordered micromechanicalstructures, which lead to large random spatial and temporal stress fluctuations. Physically accurate simulation of such materials has generally only been possible using numerically intensive methods from material science and physics [14, 2, 31], where lattice models with large
numbers of spatial degrees of freedom predominate. These models arecomputationally costly, but phenomenologically accurate.
The present contribution is enabled through recent work in the statistical physics of fracture [3]. It identified a purely temporal fracture process, reviewed below, which captures the stochastic temporal dynamics without representing spatial degrees of freedom. In thatwork, we identified a stochastic jump process exactly capturing thefluctuating pattern of failure in a fiber bundle model of fracture andobtained a factorization, separating out a mean constitutive responsethat closely matches the nonlinear stress-strain behavior of the exact
model. The latter aspect allowed the (smooth) stress-strain dynamics
to be simulated deterministically, without tracking the (rapidly fluctuating) random failure history in the material. The former provides aniterative stochastic description of instants of failure in the material as it
is loaded. The work presented here makes use, and builds upon, bothinnovations in order to enable real-time rendering of interactions withthese complex materials.
Summary of Contributions and Content
The main contributions of the present paper are:
1. A novel algorithm for the efficient, real-time rendering of multimodal signatures of fracture in natural and man-made fiber composite materials, based on a statistical model of fracture developed by the author [3].
2. A derivation of a stochastic point process model of strained or
shear granular media, based on the inverse Fiber Bundle Model.
This construction is inspired by the aforementioned one, and ispresented here for the first time.
3. An algorithm for fast multimodal rendering of fracture-like processes that arise in strained or sheared granular media.
4. Novel applications of the foregoing algorithms to the interactive,
multimodal rendering of fracture in fiber composite materials.
Interactions are enabled by a force-feedback device and (in thecase of audio) a loudspeaker.
Efficiency is guaranteed by a construction which allows a distributed physical medium to be modeled as a lumped-parameter point
process, as in item (1.) above. This, in turn, effectively reduces therendering computations to the integration of a 2nd-order dynamical(nonlinear) dynamical system coupled with a sequential random number sampling procedure. For this reason, most of the material below

is devoted to exposition of this method and its application to the material simulations of interest. Further details are provided in a separate
publication of the authors [3].
The organization of the remainder of this paper is as follows. The
next section provides a brief review of a well studied lattice statisticalphysical model of fracture processes, the Fiber Bundle Model, that has
been found to capture the salient physical phenomena in materials andprocesses like those mentioned above. Then, we present the construction of a low-dimensional stochastic process that exactly captures the
behavior of the aforementioned model in the limit of long-range loadtransfer, as demonstrated theoretically and numerically.
2

F IBER B UNDLE M ODEL OF M ATERIAL F RACTURE

Fiber Bundle Models (FBMs) are lattice statistical physical models offracture for failure processes in quasi-brittle, heterogeneous materials[2, 40]. Fracture processes in these materials unfold in unpredictableways, particularly when compared with fracture in homogeneous solidmaterials. They are accompanied by a host of distinctive phenomena that have been documented in prior literature, including statisticalstrength distributions, stress fluctuations, reorganization accompanying failure, acoustic emissions, and accumulated damage. These phenomena, while highly salient, are not well captured by deterministic
continuum mechanics models of fracture, which has led to the development of statistical approaches, including the FBM.
These models consist of lattices of N parallel elastic or viscoelastic
fibers each sharing a quantity σi , i = 1, 2, . . . , N, of the force, F, on thebundle (Figure 3). The strain xi (t) of each ith fiber is governed by thephysical equation
NN 
F(t) = ∑ σi (t) = ∑ φi σiF (t) + σiR (t) .
i=1

(1)

i=1

The strain-dependent, per-fiber load is σi = σi (xi , ẋi , ẍi , · · ·). It may
be expressed in terms of a part σiF born by intact fibers (for which theindicator variable φi = 1) and another, σiR , that is born by the surrounding (non-fibrous) material matrix, which persists after failure (φi = 0).
A minimal micromechanical model for viscoelastic and plastic effects(Fig. 3, modified Kelvin-Vogt model) may be given as follows:σiF (t) = (bF ẋi + kF xi ), σiR = (bR ẋi + mẍi )(2)

The dynamical equation depends upon an effective (per-fiber) masscoefficient m (proportional to the effective fiber mass), elastic constant
kF and damping constants bF and bR representing the pre- and postfailure relaxation of the fiber. [17]. A number of variations of thismicromechanical model are possible [40].Fibers fracture when xi (t)exceeds a fiber-specific breaking threshold ξ . The thresholds are random variables, with ξ ∼ p(ξ ), where p(ξ ) is a random distributionof fiber strengths. After fracture, the load is redistributed among thosethat survive. A common simplifying assumption (Equal Load Sharing,
ELS) is that upon the failure of an individual fiber its load is equallydistributed, so that the load on each intact fiber depends on the totalload and number NF of intact fibers σi = F/NF . We adopt this assumption here, and discuss extensions to local load sharing in Sec. 5,
below.
Each fracture event decreases the number NF (t) of intact fibers attime t, increasing the load on surviving fibers, and cascading in further failures. This continues until x(t) < ξ ∗ , where ξ ∗ is the thresholdof the weakest surviving fiber. When a critical value Fc of the applied stress is reached, the bundle is incapable of supporting the redistributed load, and all remaining fibers break. The number NF of fibers
surviving at a given load thus depends on the load history and random
assignment of thresholds.
Disorder in this model can be introduced through the initial assignment of random fiber thresholds, which model the random strengthsof micromechanical elements in the medium. These thresholds canbe viewed as random initial conditions, with the subsequent evolution
occurring deterministically. Instead, we can regard the sequence of

A

B

F(t )
x (t )i (t )x i (t )mx (t )bR

kF bF

Fig. 3. Illustration of the fiber Bundle Model of fracture. A. A lattice of N
viscoelastic-plastic fibers is loaded with external force F(t). The weakest surviving fiber strains under local stress σi . B. A mechanical analog
in the form of a modified Kelvin-Vogt element. The plastic element fails
when the force on it is greater than a random threshold ξ , a disconnecting spring kF and damper bF . Post-fracture relaxation is modeled via apersistent damping factor, bR .

failure points as a random point process whereby the threshold ξ representing the weakest fiber jumps in value at the time of each fracture.
In the next section, we provide an explicit formulation of the this process.

is surpassed. One can interpret the failure series as a sequential,
monotonically increasing random process that reproduces the specified strength distribution p(ξ ). A temporally local definition of theprocess ξ (t) that conforms to the desired strength distribution p(ξ )may obtained by sequentially sampling failure thresholds, rather thanintegrating the dynamical state of the lattice of all fibers. This canbe done using a method derived from the inverse transform samplingtechnique of random number generation [9]; we briefly review themethod, which is presented for the first time in a separate work of theauthor [3]. Let P(·) be the cumulative distribution function of p(ξ ),and let NF, j = N − j + 1 is the number of surviving fibers prior to thejth failure. Denote by ξ (t) = ξ j , where ξ j represents the strengthof the weakest surviving fiber at time t. When a fracture event occurs, ξ (t) jumps in value and the number of surviving fibers, NF (t),decreases for each failed fiber. This happens whenever the dynamic
strain x(t) exceeds ξ (t). A fracture event at time t yields a jump fromξ to a new value ξ 0 given by
s0 = s + (1 − s) 1 − u1/NF (t) ,
(3)
ξ0

=

P−1 (s0 ).

(4)Upon substituting and using ξ 0 = P−1 (s0 ), this provides a simple iterative expression for the fracture event:
ξ 0 = P−1 P(ξ ) + 1 − P(ξ ) 1 − u1/NF (t) .
(5)

ξ=P-1(x)The resulting sequence ξ j is equivalent to a set of N independent samples from p(ξ ) sorted in increasing magnitude.
The fracture size at time t depends on the state of the co-evolving
random process NF (t) and on the value of the strength ξ (t) prior to
failure. The fracture time depends on the values of ξ (t) and the strainx(t). This algorithm reproduces the ensemble of samples from p(·) by
sequentially sampling the conditional distributions p(ξ j | ξ j−1 , NF, j ),as in Eq. (5); see Fig 4. The resulting process is non-Markovian, because it depends on the sample history through NF (t), and through the
dynamical stress-strain equation (1).

ξ3
ξ2
ξ1

p(x)

p(x1)

p(x2|x1)

x1

p(x3|x2)x2

p(x4|x3)

x3

x=1

Fig. 4. The sequential sampling algorithm described in Eq. (??-??). The
next value ξ j of the threshold is determined from s j−1 , with the distribution represented by the shaded area to the right of s j , and the number ofremaining fibers, NF, j . ξ (t) is sampled in a way that depends on its valueξ (t − dt) at the preceding instant and the load history, σ (τ), τ ≤ t. The
process is non-Markovian, because the kth such distribution depends
on the number N − k − 1 of preceding fractures.

3

T IME -D OMAIN S TOCHASTIC P ROCESS F ORMULATION

The two key variables reflecting the instantaneous state of this model
are the number NF of intact fibers and the breaking threshold ξ ∗ of theweakest intact fiber. At the instant of failure, ξ ∗ increases by a random amount ∆ that is related to p(ξ ) and the number NF of surviving
fibers. The result is a stochastic jump process for a temporally fluctuating threshold ξ (t), equal at any instant t to ξ ∗ , whose jth piecewise
constant value ξ (t) = ξ j is reached at the instant ξ (t − dt) = ξ j−1

Mean damage approximation: In the foregoing construction, theglobal strain threshold ξ (t) depends on the random level of damage
at the time of fracture (represented by N − NF ), which depends on thehistory of the sample. This section describes an effective model withsmooth constitutive behavior, facilitating efficient simulation.
The instantaneous strain x(t) depends on the stochastic evolution ofdamage in the lattice. However, the real dynamics will tend to average over the fluctuations in stress. This suggests a modeling approachthat averages over fluctuations to approximate the nonlinear strain evolution deterministically, with random effects captured by the variableξ . One can achieve this by replacing NF (t), where it appears in thethreshold and evolution equations (5) and (1), by the expected survival number N F (t) given the load history. Under ELS, its value is
N F = N(1 − P(x)), with x the maximum strain achieved during loading. N F is the mean damage that would be expected for an ensembleof instances of the model subjected to the given load history.
The expected survival number N F (t) depends on the strain via
N F (t) = N (1 − P(x∗ (t))) , x∗ (t) = max
x(t 0 ).
0
t <t

(6)

Under monotonically increasing loading, this equals N F (t) = N(1 −P(x(t))). Upon replacing NF by N F one can factorize the model into
a deterministic part governing the nonlinear stress-strain response,F(t) = N F (t)σ F (t) + Nσ R (t)and a stochastic process describing the stress fluctuations:
ξ 0 = P−1 P(ξ ) + 1 − P(ξ ) 1 − u1/N F
while x(t) > ξ .

(7)

(8)
(9)

k,b,m
σ(t)

x

x(t)

D

σ

S

E

E, Δt

ξ'

t

Fig. 5. The mean damage model M is composed of a deterministic
process D describing the mean constitutive response and a stochastic
process S describing the fluctuations.

For now, refer to the original FBM model as M and the approximation
obtained through this “mean damage” replacement as M. The lattertakes on explicit form when a fiber strength distribution p(ξ ) is specified.
3.1

Statistical Distributions p(ξ ) of Fiber Breaking Thresholds

Uniform fiber strength distribution: In one commonly used fiber
strength distribution, p(ξ ) is assumed to be uniformly distributed on
the unit interval [0, 1]. Assuming a minimal Kelvin-Vogt type viscoelastic micromechanical model, the homogenized nonlinear stressstrain response becomes:
F(t) = N(1 − x) (bF ẋ + kF x) + N(bR ẋ + mẍ)(10)

while the failure event is described by sampling ξ 0 as:ξ 0 = ξ + (ξ − 1)(u1/(N(1−ξ )) − 1)(11)with u taken to be a sample of a random variable uniformly distributed
in [0, 1].
Weibull fiber strength distribution: The Weibull function is frequently used as a statistical distribution for solid strength in materials
science [48]. The Weibull CDF can be writtenP(ξ ) = 1 − exp(−ξ /λ )k Θ(ξ ).

(12)Here λ and k are scale and shape parameters, and Θ(·) is the Heaviside
step function, with Θ(ξ ) = 0 for ξ < 0 and Θ(ξ ) = 1 for ξ > 0. The
constitutive model is
σ (t) = N exp{(−x/λ )k } (bF ẋ + kF x) + N(bR ẋ + mẍ)(13)

and the fracture event is sampled as
 k


1/k
k
k
ξ
log 1 − e(−ξ /λ ) 1 − u1/(N exp{(−ξ /λ ) })(14)
ξ0 = λλConstitutive behavior and fluctuations: The value of the critical stress Fc and distributions of fluctuations as failure approaches are
known to depend weakly on the precise distribution of fiber strengths[2, 40]. Specifically, the constitutive evolution of M differs from thatof M due to fluctuations in the survival number NF about its mean.
This can be regarded as a source of high-frequency noise that should
approximately integrate to zero, so we reasoned that even if these fluctuations are significant, the model M would yield similar behavior to
M. Indeed, this turns out to be the case; detailed numerical comparisons are provided in a separate publication. Simulation runs are qualitatively indistinguishable for both models. Figure 7 provides represntative results (samples of model M), and a comparison of constitutive behavior and failure thresholds for both models is given in 6. In

Fig. 6. Left: The empirical stress-strain behavior of models M and M.
Both were obtained under quasistatic loading, at supra-critical values ofthe maximum stress, F > Fc , computed from 200 uniformly distributed
values of the terminal load, using simulations with 2000 fibers. Dashedblue line: critical stress values of Fc (M) = 2800, Fc (M) = 2850. Dashed
red lines: 95% confidence intervals. Right: Critical stress values obtained from the same simulations. Results from the study [3].

these simulations, the per-fiber load increases linearly until the current
threshold level ξ (t) is met, at which point one or more fibers fail anda new value of ξ (t) is sampled. The size and density of these failure
events increase with the applied load, until the model reaches a stable
point, or a critical load is achieved, causing the bundle to collapse.
Stress-strain relationships and critical load estimates are comparedfor both models in Figure 6. Qualitatively and quantitatively, both arenearly identical, with an error of less than 2% in the critical load Fc forall bundle sizes examined (size N = 1000 to 256000), indicating thatthe mean macroscopic behavior of the model is well-preserved by theapproximation.
Figure ?? provides evidence that this construction preserves themost salient physical phenomena captured by the FBM. It shows thatthe empirical distributions of failure event time intervals dt and of energy fluctuations dE exhibit power-law scaling with exponents consistent with prior literature [13, 40, 22],
4

S TRESS FLUCTUATIONS IN COMPRESSED GRANULAR ME DIA : I NVERSE F IBER B UNDLE M ODEL

Complex patterns of force fluctuations accompany the compressionof stiff granular materials, which, as in the case of fiber composite
materials, have been successfully modeled in terms of the inelastic
failure of one-dimensional structures loaded in parallel [15, 18]. The
latter arise due to the organization of stress in granular materials inhighly heterogeneous force networks, in the form of force chains percolating prevailingly along the direction of the principle axis of stress.
These chains are created and restructured due to relative displacementsof the grains, which rearrange their positions until a stronger, stableconfiguration is reached. Energy is rapidly dissipated in granular assemblies, and hence these rearrangements occur rapidly and highly
inelastically, giving rise to large transient stress fluctuations that provide valuable information about the evolution of force chains duringloading [23, 15].
4.1

Inverse FBM

A model for the evolution of force chains during compression was
proposed by Hidalgo et al. [15, 18], as an extension of the continuous
damage Fiber Bundle Model [24, 16]. In experiments on ensembles of
granular particles, the model was able to effectively capture both nonlinear constitutive behavior of the medium and predicted the amplitude
distribution of acoustic signals emitted by microscopic rearrangementsof the particles during compression.
In this model, granular force chains are represented by an arrayof fibers loaded in parallel with randomly distributed failure thresholds. However, in contrast to the standard FBM, chains restructure
after failure, becoming multiplicatively stronger and stiffer. If k is thepre-failure effective stiffness, and k0 is the value after failure, thenk0 = ak = k(1 + α), α > 0.

(15)

σ

σ'
k, ξ

k', ξ'

Fig. 8. Compression results in rapid events through which granular force
chains become stronger, stiffer and straighter.

small (in the simulations below, α = 0.005). Thus, one can approximate
NA(t) ≈ ∑ (1 + α n j (t)) = N + α IF (t),

(19)

j=1

Fig. 7. Response to a linear ramp load, F(t) = F0 r(t), where r(t) = t/τ0
for 0 ≤ t ≤ τ0 = 0.8. N = 2000 fibers, and Weibull-distributed thresholds(parameters k = 4 and λ = 2.5). Top row: Survival fraction NF (t)/N vs t.
Maximum stress, F0 , is shown relative to critica l stress, Fc ≈ 2800. Here,F0 = 1.2Fc . Failure occurs at t = 0.68. Bottom row: Threshold ξ (t) andstrain x(t) vs t.

For this reason, it has been called the inverse Fiber Bundle Model.
The bundle grows progressively stronger, rather than weaker (asin the FBM), corresponding to an increase in packing density of themedium. The constitutive dynamics of the bundle can be described,
using the parametrization of Section 3, by
σ (t) = A(t)(bF ẋ + kF x), σ (t) = (N + αIF )(bF ẋ + kF x)

(16)

where

where IF (t) is the total number of failure events until time t.
The breaking threshold of a fiber increases upon failure. Thus, because all threshold are assumed to be drawn from the same distribution, the sequential algorithm of Section 3 can be used for simulation.
After a failure, the next lowest threshold of the model is drawn according toξ 0 = P−1 P(ξ ) + 1 − P(ξ ) 1 − u1/NT.
(20)
In the mean damage approximation, the discrete random variables
IF (t) and NT (t) appearing in (16) and (20) may plausibly be replacedby their expected values, based on arguments similar to those used in
the foregoing, yielding
σ (t) = (N + α I F (t))(bF ẋ(t) + kF x(t)).

The expected number of failures I F (t) at a given strain level is givenby
I F (t) = nmax NP(x∗ (t)), x∗ (t) = max
x(t 0 ).
(22)
0
t <t

The expected survival number is, as for the Fiber Bundle Model, given
by N F (t) = N(1−P(x∗ (t))). Finally, the expected value of the numberof remaining failures in the bundle is N T = N F nmax . The constitutive
and fluctuating response of this model can be written as
σ (t) = N(1 + α nmax P(x∗ ))(bF ẋ + kF x).

N

A(t)

=

∑ (1 + α)n (t) ,
j

(17)

N

=

kF ∑ (1 + α)n j

(23)

and

j=1

K

(21)

(18)

0

−1

ξ =P




[Nnmax (1−P(x∗ ))]−1.(24)
P(ξ ) + 1 − P(ξ ) 1 − u

j=1

IFBM: Time-domain stochastic process formulation

Failure criteria For a material sample under a uniform, uniaxial load, either a maximum stress or maximum strain criterion may be
used (the latter being identical to that employed for the fiber composite
model of Section 3). In a local model of failure in a three-dimensional
volume of granular material, the breaking threshold depends on themagnitude and directions of stresses or strains, because the load bearing force chains in the material reorganize to align with the directionof maximum stress. For coarse-grained soils or other granular media, failure can often be described by a Mohr-Coulomb yield criterion.
This can be written as
s/n = tan(φ ),(25)

Although the inverse FBM is amenable to theoretical analysis, thepresent contribution focuses on applying an analogous modeling andsimulation method to that of Section 3. The value of α, the relative
increase in stiffness at each fracture, can generally be assumed to bewhere s/n is the ratio of maximum and minimum principle stresses
at the location of interest, and the parameter φ , the angle of repose(or angle of internal friction) of the granular material, describes thematerial strength.

and n j (t) is the number of failures of the jth fiber at time t. Here, it isassumed that the damping ratio is preserved during evolution, so thatb0 = b(1 + α).
When a given fiber has failed a maximum number nmax of times(nmax is a parameter of the model), it is assumed to be unable to further
restructure, attaining a maximum stiffness of kmax = kF (1 + α)nmax .
Subsequently, that fiber contributes a linear response. When all fibers
have failed nmax times, the bundle enters a linear response regime.
4.2

4.3 Numerical validation
The mean damage inverse Fiber Bundle Model presented above wassimulated under uniaxial, stress-controlled loading, in the same manner as in Section 3.1. A Weibull distribution of breaking thresholds
was used with the same values of the parameters as used in the previous simulations. Figure 9 shows typical response profiles of the model
under ramp loading conditions. As noted above, the nonlinear stressstrain relation is, opposite to the case of the Fiber Bundle Model, concave, corresponding to the increasing stiffness of the granular medium
as it is compressed.
Applied stress σ(t)2.5(1-Nf ) / N

2.0
σ(t)

1.5
1.0

0.4

t

σ(t)

0.0

0.8

0

4000

σ (x)

6
x(t)

1000

0.2

4

0.4

t

Stress vs. strain

0

3000

.020

2000

.015

1000

2
0.4

t

0.8

500

0.8

E(t)

0

Strain x(t)

0
0

0.3

0.1

0.5
0

Fractures per fiber

0.4

which is consistent with published results [15]. Theory predicts critical scaling in this model to occur when the maximum number of
permitted rearrangements nmax of each fiber is large, and theoretical
calculations of the values of the scaling exponents have been foundto agree very well with experimental observations of fluctuations in acompressed granular system [15].
An example showing the response of the inverse Fiber Bundle
Model to a dynamically variable applied stress (load profile recorded
from a human footstep), is shown in Figure 11. Although the behavior observed is plausibly realistic, with most fracture energy produced
during initial loading (heel strike) and in stress buildup prior to unloading (toe off), a more realistic response would be obtained if thespatially distributed nature of foot-ground contact were accounted for.

0
0

.01
0.0

2

x

0

4

Fig. 9. Response of the inverse Fiber Bundle Model under ramp loading conditions. Top Left: Applied stress vs time. Top Right: Numberof fracture events per fiber vs time. Bottom Left: Strain x(t). Bottom
right: Constitutive, stress-strain behavior. Data shown is based on 200
simulations of bundles with 2000 fibers subjected to ramp loading withmaximum stress σ = 2667. Other simulation parametesr were as given
in Figure 7.

Fig. 10. Left: Empirical distribution of times between failure events, ∆t,
for the inverse Fiber Bundle Model. Right: Distribution of failure energy fluctuations (computed in windows of duration 0.005), exhibiting
power-law scaling over more than one order of magnitude. Inset: Energy fluctuation size vs time during one simulation trial. Distributions
were computed from 500 simulations of bundles with 2000 fibers subjected to ramp loading with maximum stress σ0 = 105 and nmax = 20.
Other simulation details were as given in Figure 7.

t1

Fig. 11. The response of a granular medium, represented by a bundleof 2000 iFBM fibers, to a dynamic load, consisting of recorded data froma human footstep. Shown are the applied load σ (t), and the energy E(t)computed in windows of duration ∆t = 0.005.

Distributed and spatially varying media This paper has so far
treated interactions with heterogeneous materials with material properties that are constant, in a statistical sense, within the simulated volume. This work could also be extended to address to interactive simulation of fracture processes that vary in space, or that transfer load
preferentially to nearby fibers [19], as in large volume of distributed
granular material (Fig. 12). This is facilitated by the mean damage
approximation. In a distributed setting, the averaged nonlinear strain
response (Eq. (23) for the inverse Fiber Bundle Model) can be obtained
from a homogenized continuum mechanics model of the constitutive
behavior of the distributed medium, simulated using finite element,
elastic mesh, or point cloud methods. Using our model, fracture processes can be simulated by associating independent instances ξk (t) of
the fracture process effective volume elements (at points xk ) definedby finite elements or by point based radial basis functions Φi (x). Conditions for local fracture can then be obtained from a measure of maximal strain obtained from the strain tensor ε, similar to criteria used
to nucleate and grow cracks in continuous media for computer graphics rendering [39, 36, 37]. Since the constitutive behavior is smooth,it can be integrated at a lower sample rates than those used to render
fracture events.
5

Figure 10 presents empirical distributions of the times ∆t between
failure events and of the energy fluctuations of the model during loading. Both were computed as in Section 3.1. Energy fluctuations areobserved, for large enough values of nmax , to scale as p(E) ∼ E −2.0 ,0.5

I NTERACTIVE M ULTIMODAL
NEOUS F RACTURE

R ENDERING OF H ETEROGE -Interactive simulation of failure processes in quasi-brittle, heterogeneous materials is of interest for scientific visualization or immersive
virtual reality. By virtue of the stochastic formulation reviewed above,

Force3.0

x 10 4

2.0

F(t)

1.0
0.0

Strain

1.2

xk

0.8

x(t)

0.4

1.0

0.0

0.0

it was possible to implement the algorithms efficiently, at audio sampling rate ( fs = 48 kHz). This made it possible to reproduce stressfluctuations accompanying fracture evolution across the entire vibrotactile and audible bandwidths (Fig. 13). Figure 15 shows a still framefrom an interactive simulation, with control and force feedback mediated via a haptic interface (Falcon, Novint Inc., with libnifalcon open
source drivers), and with concurrent auditory feedback. The simulation runs in real time at high sampling rates (dynamics and audio: 48
kHz, haptic control loop: 1 kHz).
A variable stress σ (t) (either simulated or acquired via instrumentation) is provided as a real-time input to the dynamics, resulting ina time-varying strain x(t), computed as the solution of Equation (1).
Fracture events occur as described in Eq. (11)-(12). Each event evokes
a sequence of transient failures with energies E = ∆NkF x(t)2 , where∆N is the number of fibers that fail in a given time-step.
In the implementation, discretization of the dynamics was performed via Laplace transformation of the ordinary differential equations, followed by the bilinear transform mapping from the s- to the zplane. This allows the constitutive response of the system to be implemented via low-order infinite impulse response (IIR) discrete time filters. This explicit integration method provides a reasonable accuracyperformance tradeoff at high (audio bandwidth) sampling rates [44],far better than is achieved with a finite-difference discretization ofcomparable order. At low frequencies, both finite-difference and bilinear transform discretizations produce essentially equivalent results.
The fracture algorithm requires minimal computation. When fracture events are sparse, most computation time is spent on the integration; for model M, the stress-strain dynamics may be implemented
as a second order IIR filter requiring five multiply-accumulate operations per sample period. As the density of fracture events per sample
period increases, computation time is increasingly spent updating the
breaking threshold ξ after each failure, which requires calculating theinverse fiber strength CDF P−1 (·). However, by interpolating values
stored in a function table or by precomputing sampled thresholds, thiscost is reduced to a small, constant number of operations per fracture
event.
Random Modal Filter Model of Acoustic Transfer The qualityof mechanical vibrations or acoustic emissions from fracture processes
are strongly determined by the transfer properties of the medium. The
latter strongly influence the perception of these acoustic events, especially for natural materials such as wood [11, 21]. Thus, the rendering algorithm includes a linear filter model of acoustic transfer in themedium (Fig. 14). It is represented by a modal resonant filter with
impulse response [7]Mh(t) = ∑ Ai e−bi t sin(2π fi t).
i=1

(26)

0.8

-1.0

ξ(t)

0.4

0.99

t

1.0

0.0
8 kHz

Frequency, fFig. 12. Rendering high-frequency failure processes in distributed media at interactive rates is facilitated by the separation they provide between the slower, homogenized volumetric response and the rapid fluctuations, which can be represented through local, stress-driven fracture
processes associated with points xk .

Threshold

1.2

4 kHz

0 Hz
0

0.5
Time, t1

Fig. 13. Representative data from a real-time implementation of Fiber
Bundle Model M. Data shown corresponds to a model with 350 fibers.
Acoustic transfer is simulated by a resonant filter with N = 10 modes.
Left column, Top to Bottom: Force input F(t), strain x(t), threshold ξ (t),and time-frequency spectrogram of acoustic energy of fracture. Right:
Fracture event waveform.

Ampl, A

A∝f
i
i

X

X

-r

X

X

X

fi ~ p( f ) ∝ f

ɑ

X

X

freq, fFig. 14. The magnitude frequency response of a random modal filter
model used to model acoustic transfer in the simulated medium captures the physical phenomena of geometric growth of mode density anddamping-dependent decrease in amplitude with frequency, using justtwo key parameters.

Although the resonant filter parameters could, in principle, be precomputed by modeling the material geometry and constitution, a simpler approach is adopted here. Resonant frequencies fi are selected to
be randomly distributed as p( f ) ∝ f α , where the parameter α controls
the rate of growth of mode density with frequency. In real materials, this is related to the geometry, dimension, and wave velocity ofthe propagating medium. Decay factors bi describe the filter bandwidth, which is related to the intrinsic damping in the material. The
amplitudes Ai are chosen to provide a decaying spectrum envelopewith Ai = A0 fi−r . This rate of amplitude decrease with frequency
corresponds to a decay of 6r dB/octave [44], modeling frequencydependent damping in acoustic propagation. By choosing the number

tuations in hundreds to hundreds of thousands of distributed fibers inreal-time. This could aid a wide range of applications in virtual reality simulation of natural environments, and could enable others inscientific simulation and in interactive visualization.
R EFERENCES

Fig. 15. Interactive haptic and auditory rendering of fracture in a woodbeam, using the physical model presented here. Data shown corresponds to a model with 700 fibers. Acoustic transfer is simulated by aresonant filter with N = 10 modes. The mesh shown is a mass-springmodel to which we apply strain computed by the model for visualization purposes. The beam bending in the mesh-spring model is driven
by the lumped parameter stochastic fracture model, for the purpose ofvisualizing strain.

Nm of modes, their decay rate r with frequency, and their distribution
α appropriately, the acoustic response can be shaped to qualitatively
resemble the material of interest (e.g., wood, stone; see Supplementarymaterial).
6

C ONCLUSIONSThis paper presented algorithms for simulating failure in quasi-brittle
heterogeneous materials, such as natural fiber composites and granular materials, based on recently developed formulations of two wellestablished lattice physical models: the Fiber Bundle Model and Inverse Fiber Bundle Model. The constructions involved accurately reproduce the deterministic and stochastic response of the exact algorithms, in agreement with known physics. In short, these algorithms
can accurately reproduce the most salient features of rapid fracture infiber composite and granular materials. The main purpose of introducing these constructions here is that they facilitate the efficient simulation of these complex processes for the first time, by eliminating theneed to dynamically integrate large numbers of spatial degrees of freedom, and by homogenizing the constitutive (stress-strain) behavior.
The latter has the further effect of isolating the fast varying dynamical
degrees of freedom, preventing them from contributing transient fluctuations to the dynamical state, and reducing the computation requiredfor integrating the stress-strain dynamics to a small, constant numberof operations per sample period.
The method applied here considers failure due to a load appliedalong a single dimension, which could be viewed as a limitation. Forfiber composite materials (such as wood), stresses with arbitrary directions can readily be accommodated in the model presented here bydecomposing applied stress according to the fiber geometry. Granularmaterials rearrange load bearing structures in more complicated ways
as the direction of applied stress is varied. As such, more work isneeded in order to adapt the method presented here to such cases.
The model presented here is efficient enough to simulate stress fluc-[1] F. Adjémian and P. Evesque. Different regimes of stick-slip in granular matter: from quasi- periodicity to randomness. EGS XXVII General
Assembly, Nice, 21-26 April 2002, 2002.
[2] M. Alava, P. Nukala, and S. Zapperi. Statistical models of fracture. Advances in Physics, 55(3):349–476, 2006.
[3] Y. Visell and G. Millet. Fracture in disordered heterogeneous materials
as a stochastic process. arXiv:1408.530.3 (cond-mat.stat.mech), 2014.
[4] A. Baldassarri, F. Dalton, A. Petri, S. Zapperi, G. Pontuale, and
L. Pietronero. Brownian forces in sheared granular matter. Physical review letters, 96(11):118002, 2006.
[5] N. Bell, Y. Yu, and P. Mucha. Particle-based simulation of granular materials. In Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, page 86. ACM, 2005.
[6] R. Chudoba, M. Vořechovskỳ, and M. Konrad. Stochastic modeling ofmulti-filament yarns. i. random properties within the cross-section andsize effect. International Journal of Solids and Structures, 43(3):413–
434, 2006.
[7] P. Cook. Real Sound Synthesis for Interactive Applications. AK Peters,Ltd., 2002.
[8] F. Dalton, F. Farrelly, A. Petri, L. Pietronero, L. Pitolli, and G. Pontuale.
Shear stress fluctuations in the granular liquid and solid phases. Physicalreview letters, 95(13):138001, 2005.
[9] L. Devroye. Non-Uniform Random Variate Generation. Springer Verlag,
1986.
[10] G. Dill-Langer, R. C. Hidalgo, F. Kun, Y. Moreno, S. Aicher, and H. Herrmann. Size dependency of tension strength in natural fiber composites.
Physica A: Statistical Mechanics and its Applications, 325(3):547–560,
2003.
[11] B. L. Giordano and S. McAdams. Material identification of real impact
sounds: Effects of size variation in steel, glass, wood, and plexiglass
plates. The Journal of the Acoustical Society of America, 119(2):1171–
1181, 2006.
[12] L. Glondu, M. Marchal, and G. Dumont. Real-time simulation of brittle
fracture using modal analysis. Visualization and Computer Graphics,
IEEE Transactions on, 19(2):201–209, 2013.
[13] A. Hansen and P. Hemmer. Burst avalanches in bundles of fibers: Local
versus global load-sharing. Physics Letters A, 184(6):394–396, 1994.
[14] H. Herrmann and S. Roux. Statistical models for the fracture of disordered media. North Holland, 1990.
[15] R. Hidalgo, C. Grosse, F. Kun, H. Reinhardt, and H. Herrmann. Evolution of percolating force chains in compressed granular media. Physicalreview letters, 89(20):205501, 2002.
[16] R. Hidalgo, F. Kun, and H. Herrmann. Bursts in a fiber bundle modelwith continuous damage. Physical Review E, 64(6):66122, 2001.
[17] R. Hidalgo, F. Kun, and H. Herrmann. Slow relaxation of fiber composites, variable range of interaction approach. Physica A: Statistical
Mechanics and its Applications, 347:402–410, 2005.
[18] R. C. Hidalgo, H. J. Herrmann, F. Kun, and E. Partelli. Force chains in
granular packings. In F. Mallamace and H. Stanley, editors, The Physics
of Complex Systems. IOS Press, Amsterdam, 2004.
[19] R. C. Hidalgo, Y. Moreno, F. Kun, and H. J. Herrmann. Fracture modelwith variable range of interaction. Physical review E, 65(4):046148,
2002.
[20] K. Hirota, Y. Tanoue, and T. Kaneko. Generation of crack patterns with aphysical model. The Visual Computer, 14(3):126–137, 1998.
[21] R. L. Klatzky, D. K. Pai, and E. P. Krotkov. Perception of material
from contact sounds. Presence: Teleoperators and Virtual Environments,
9(4):399–410, 2000.
[22] M. Kloster, A. Hansen, and P. C. Hemmer. Burst avalanches in solvablemodels of fibrous materials. Physical Review E, 56(3):2615, 1997.
[23] R. Koerner, W. McCabe, and A. Lord Jr. Acoustic emission behavior andmonitoring of soils. In Acoustic emissions in geotechnical engineering
practice: a symposium, page 93. ASTM International, 1981.
[24] F. Kun, S. Zapperi, and H. J. Herrmann. Damage in fiber bundle models. The European Physical Journal B-Condensed Matter and Complex
Systems, 17(2):269–279, 2000.

[25] T. Lenaerts, B. Adams, and P. Dutré. Porous flow in particle-based fluid
simulations. ACM Transactions on Graphics (TOG), 27(3):49, 2008.
[26] M. Lin and M. Otaduy. Haptic rendering: foundations, algorithms, andapplications. AK Peters Ltd, 2008.
[27] N. Liu, X. He, S. Li, and G. Wang. Meshless simulation of brittle fracture.
Computer Animation and Virtual Worlds, 22(2-3):115–124, 2011.
[28] K. Loades, A. Bengough, M. Bransby, and P. Hallett. Planting density
influence on fibrous root reinforcement of soils. Ecological Engineering,36(3):276–284, 2010.
[29] S. Mérillou and D. Ghazanfarpour. A survey of aging and weathering
phenomena in computer graphics. Computers & Graphics, 32(2):159–
174, 2008.
[30] G. Michlmayr, D. Cohen, and D. Or. Sources and characteristics of acoustic emissions from mechanically stressed geologic granular media—a review. Earth-Science Reviews, 112(3):97–114, 2012.
[31] L. Mishnaevsky Jr and P. Brøndsted. Micromechanical modeling of damage and fracture of unidirectional fiber reinforced composites: A review.
Computational Materials Science, 44(4):1351–1359, 2009.
[32] M. Müller, L. McMillan, J. Dorsey, and R. Jagnow. Real-time simulationof deformation and fracture of stiff materials. Springer, 2001.
[33] S. Nasuno, A. Kudrolli, A. Bak, and J. Gollub. Time-resolved studies of stick-slip friction in sheared granular layers. Physical Review E,
58(2):2161–2171, 1998.
[34] S. Nasuno, A. Kudrolli, and J. Gollub. Friction in granular layers: hysteresis and precursors. Physical Review Letters, 79(5):949–952, 1997.
[35] A. Norton, G. Turk, B. Bacon, J. Gerth, and P. Sweeney. Animationof fracture by physical modeling. The Visual Computer, 7(4):210–219,
1991.
[36] J. O’Brien and J. Hodgins. Graphical modeling and animation of brittle
fracture. In Proc of the 26th annual conference on Computer graphics
and interactive techniques, pages 137–146. ACM Press/Addison-Wesley
Publishing Co. New York, NY, USA, 1999.
[37] E. Parker and J. O’Brien. Real-time deformation and fracture in a gameenvironment. In Proceedings of the 2009 ACM SIGGRAPH/Eurographics
Symposium on Computer Animation, pages 165–175. ACM, 2009.
[38] E. G. Parker and J. F. O’Brien. Real-time deformation and fracture in a game environment. In Proceedings of the 2009 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pages 165–
175. ACM, 2009.
[39] M. Pauly, R. Keiser, B. Adams, P. Dutré, M. Gross, and L. Guibas. Meshless animation of fracturing solids. Proceedings of the 2005 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 24(3):957–
964, 2005.
[40] S. Pradhan, A. Hansen, and B. Chakrabarti. Failure processes in elastic
fiber bundles. Reviews of Modern Physics, 82(1):499–555, 2010.
[41] N. M. Pugno, F. Bosia, and T. Abdalrahman. Hierarchical fiber bundlemodel to investigate the complex architectures of biological materials.
Physical Review E, 85(1):011903, 2012.
[42] I. Reiweger, J. Schweizer, J. Dual, and H. Herrmann. Modelling snow
failure with a fibre bundle model. Journal of Glaciology, 55(194):997–
1002, 2009.
[43] A. Rinaldi. Statistical model with two order parameters for ductile andsoft fiber bundles in nanoscience and biomaterials. Physical Review E,
83(4):046126, 2011.
[44] J. Smith. Physical audio signal processing. J. O. Smith, 2004.
[45] J. Smith, A. Witkin, and D. Baraff. Fast and controllable simulation of
the shattering of brittle objects. In Computer Graphics Forum, volume 20,
pages 81–91. Wiley Online Library, 2001.
[46] S. M. Sutherland. Fibre based modeling of wood dynamics and fracture.
PhD thesis, University of British Columbia, 2012.
[47] D. Terzopoulos and K. Fleischer. Deformable models. The visual computer, 4(6):306–331, 1988.
[48] W. Weibull. A statistical distribution function of wide applicability. Journal of applied mechanics, 1951.
[49] C. Zheng and D. L. James. Rigid-body fracture sound with precomputed
soundbanks. ACM Transactions on Graphics (TOG), 29(4):69, 2010.
[50] Y. Zhu and R. Bridson. Animating sand as a fluid. In ACM SIGGRAPH
2005 Papers, page 972. ACM, 2005.