# ‘Damage cascade’ in a softening interface

Arnaud Delaplace, Stéphane Roux, Gilles Pijaudier-Cabot

► To cite this version:

Arnaud Delaplace, Stéphane Roux, Gilles Pijaudier-Cabot. ‘Damage cascade’ in a softening interface. International Journal of Solids and Structures, 1999, 36 (10), pp.1403-1426. 10.1016/S0020-7683(98)00054-7. hal-01006714

HAL Id: hal-01006714

https://hal.science/hal-01006714

Submitted on 31 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ‘Damage cascade’ in a softening interface

Arnaud Delaplace<sup>a\b</sup> Stéphane Roux<sup>b</sup> Gilles Pijaudier-Cabot<sup>c</sup>

<sup>a</sup> LMT, ENS Cachan/CNRS/Université P. et M. Curie, 61, Avenue du Président Wilson, F-94235 Cachan, France

<sup>b</sup> Laboratoire de Physique et Mécanique des Milieux Hétérogènes, Ecole Supérieure de Physique et de Chimie Industrielles, 10 rue Vauquelin, F-75231 Paris Cedex 05, France

<sup>c</sup> LMT, ENS Cachan/CNRS/Université P. et M. Curie, 61, Avenue du Président Wilson, F-94235 Cachan, France and Institut Universitaire de France

A model describing the damage at an interface which is coupled to an elastic homogeneous block is introduced. Resorting to a real-space renormalization analysis, we show that in the absence of heterogeneity localization proceeds through a cascade of bifurcations which progressively concentrates the damage from the global interface to a narrow region leading to a crack nucleation. The equivalent homogeneous interface behaviour is obtained through this entire cascade, allowing for the analysis of size effects. When random heterogeneities are introduced in the interface, prior to the onset of localization damage proceeds by a sequence of avalanches whose mean size diverges at the first bifurcation point of the homogeneous interface. The large scale features of the bifurcation cascade are preserved, while the details of the late stage are smeared out by the randomness.

## 1. Introduction

Progressive failure of quasi-brittle materials can be separated in three different phases: first, the material response is elastic, then microcracking appears and these microcracks coalesce eventually in order to form a macro-crack which propagates suddenly. From the theoretical point of view, the difficulties involved in the transition between the last two phases are quite important. In the second phase, the strain field is quasi-homogeneous at a macroscopic scale. Then, the strain field becomes more and more heterogeneous, and the strain grows only inside a narrow region. The subsequent apparition of a discontinuity is often called strain localization in a general sense. There are in the literature different approaches to the description of this transition. One of them is the continuum approach, e.g. with continuous damage models. It is based on the description of the average behaviour of the material (see e.g. Krajcinovic and Lemaitre, 1987; Laws and Brock-
enbrough, 1987; Lemaitre, 1992). For such continuum models, the transition is viewed as a bifurcation problem. When strain localization is due to strain softening, the tangent stiffness operator ceases to be definite positive. The partial differential equations of equilibrium lose their ellipticity which authorizes discontinuous rate of deformation fields to develop suddenly. The inception of strain localization might be for instance depicted under some restrictive assumptions by Hill's criterion (Hill, 1959):

$$ \det [\mathbf{n} \cdot \mathbf{H} \cdot \mathbf{n}] = 0 \qquad (1) $$

where **H** is the tangent stiffness operator at the continuum point level and **n** is the orientation of the localized band. Another one is the loss of stability at the material level in the sense of the Drucker postulate. Note that in some well-defined cases (associative constitutive laws), the loss of uniqueness coincides with the loss of stability in the rheological sense.

The second approach is discrete random modelling. It is directed towards the description of the study of the material heterogeneities (i.e. at a scale lower than the representative volume of the material) (see e.g. Delaplace et al., 1996; Fokwa, 1992). Because of heterogeneity, the usual employed localization criteria in continuum models cannot be used for two reasons mainly: the first one is that the solution is always unique. To some extent, the situation for discrete models is the same as the situation for some rate dependent models where bifurcation is not possible and strain localization cannot be viewed as a loss of uniqueness problem anymore (see e.g. Dudzinski and Molinari, 1991; Leroy, 1991). The second one is that no tangent operator can be calculated because of the discrete characteristic of the response, and because of the fluctuations that appear all along the curve. There lies a subtle difficulty: consider the dimensionless ratio $\epsilon = a/L$, of the microstructure units $a$ of a discrete model over the system size $L$, which characterizes the discreteness of the medium. When $\epsilon$ tends to 0, a continuum description is expected to hold. As we will see later, the stress-strain response of the system converges towards a smooth law whose tangent operator **H** can be defined. The latter does provide information on the stability of the structure. However, when stability is analysed using actual responses for non-zero $\epsilon$, it can be shown that the fluctuations in the stress-strain responses give rise to a non-differentiable law, and this feature brings some useful additional informations on the approach to the loss of stability. Elaborating over these notions leads to the useful concept of 'avalanches'.

The aim of this paper is to propose a tool to characterize the transition from a homogeneous state of microcracking to a localised one for the discrete models. This tool should also be applied to any response with fluctuations, like those met in experiments where dispersions and fluctuations due to material heterogeneity are unavoidable. Because one cannot deal in this case with a loss of uniqueness, this tool should be based on stability considerations in the broad sense. Therefore, we study the fluctuations that are encountered all along the response of the system. More precisely, the avalanche statistics of the fluctuations are analysed. This formalism is used in many kinds of model (Bak, 1996; Paczuski et al., 1995), from biological diffusion up to earthquake response. All these models have at least one thing in common: their evolutions are structured around a critical point, as for the fibre bundle model that we used as a basis for our different proposed models.

For the sake of simplicity, we will consider a model problem: the case of a band made of a strain softening (discrete) material assembled in series to an elastic block (Fig. 1). This system is loaded by uniaxial tension, perpendicular to the direction of the band. Thus, the problem of localization will be strictly one directional as the orientation of the localised band is fixed. Since
<figure><img src="image_2.png" /><figcaption>Fig. 1. The model problem.</figcaption></figure>

the band has a finite width which might become small with respect to the block dimensions, its response may also be regarded to be the same as that of a softening interface located in between a rigid substrate and an elastic body.

In the first part of the paper, we will recall the analytical results obtained on the fibre bundle model, also called Daniels model. We will particularly present the properties of the avalanches distribution in the presence of fluctuations due to the variability of the fibres strength and deal with a simple derived model, that is a Daniels model and a spring connected in series. We will look at the evolution of the avalanche properties, and apply them for detecting the loss of stability. In order to have a realistic representation of the local mechanical redistribution of the stress field when a micro-crack appears, we will use in the second part a hierarchical model that takes into account redistribution, i.e. a non-local load sharing on the surviving fibres when a bond breaks. Again, we will look at the evolution of the avalanche properties, and we will carry out a complete study in terms of stability. In order to better understand the mechanism of rupture and the apparition of successive bifurcation points, this model will be compared to the equivalent continuum one. This model also allows to have direct access to the damage profile at the onset of unstable propagation of a macrocrack.

## 2. The Daniels model and avalanche statistics

The Daniels model (Daniels, 1945) albeit simple, displays an amazingly rich behaviour which is—at least partly—representative of the role of heterogeneity in the mechanical behaviour of some materials. It is commonly called the fibre bundle model. N parallel fibres are equally stretched between two rigid beams. The fibre behaviour is elastic up to a threshold force where the fibre breaks irreversibly. The stiffness is the same for all fibres—and thus can be chosen to be unity—but the threshold force *t* is a random variable characterized by its probability distribution function *p*(*t*), or its cumulative distribution *P*(*t*) = ∫_{0}^{*t} *p*(*t*) *dt*'. The advantage of this model is that it is completely solvable analytically. For instance, the mean force *F*, that is the applied force divided by the total number of fibres *N*, vs displacement *u* is easily obtained as

$$F(u) = (1 - P(u))u \qquad (2)$$

Then, for a stiffness of 1 (*t* = *u* for a single fibre at failure), the displacement *u* varies between 0
<figure><img src="image_3.png" /><figcaption>Fig. 2. The response of a 500-fibre Daniels' model, with a uniform distribution between 0 and 1 for the thresholds. The thin curve is the mean response, that is a parabola. A sketch of the model is included inside.</figcaption></figure>

(free interface) and 1 (interface failure). As can be seen on Fig. 2, the response of a 500-fibre bundle follows closely the theoretical expectation, i.e. a parabola for a uniform distribution $p(t) = 1$ of threshold forces between 0 and 1, but all along the curve are superimposed fluctuations. Note that similar fluctuations are also encountered experimentally on quasi-brittle heterogeneous materials, like fibre-reinforced concrete, but they are usually not described by continuum models. It is to be noted right away that the amplitude of these fluctuations vanishes as $N^{-1/2}$, and thus considering the limit of an infinite system size, $N \to \infty$, the response of the system converges to the above given mean behaviour. It is our aim to show that these fluctuations, albeit of modest amplitude, are of interest both from an experimental and a theoretical standpoint, and that some care has to be taken when considering the infinite size limit.

Since we are interested in stability, it is important to incorporate in the analysis the boundary and load conditions. In the following study, we will consider that the bundle is loaded with a testing machine of known stiffness $k$. Hence, the analysis will be applied to a bundle connected in series with a spring whose stiffness is that of the testing machine as shown in Fig. 3. This simple system can also be seen as a rough model for the mechanical behaviour of an elastic body (the spring) attached to a rigid substrate through a damageable interface (the fibre bundle). The overall displacement (bundle plus spring) will be controlled during the loading sequence. If the bundle is loaded with a stiff enough testing machine, only one fibre may break for a constant loading. However, if the stiffness is reduced, one single failure may induce catastrophically a sequence of failures before reaching a new stable position. This sequence is what we call an 'avalanche'. In the numerical simulations, one can easily solve for the response of the bundle with an ideally stiff control, that is, imposing strictly a prescribed displacement. From such a response, one can also compute the response of the bundle under any boundary conditions (including a finite stiffness $k$, and large viscous damping to avoid interial overshoot), as a succession of equilibrium positions.

Let us first consider the limit of an infinite system size, and substitute a deterministic damageable
<figure><img src="image_4.png" /><figcaption>Fig. 3. The response of a 500-fibre Daniels' model connected in series with a spring, with again an uniform distribution. The thin curve is the mean response. A sketch of the system is included inside.</figcaption></figure>

interface instead of the bundle. The damage law of the interface is chosen to be the asymptotic mean force-displacement response of the bundle, i.e.

$$F = u(1-u) \tag{3}$$

as obtained above for a uniform distribution of fibre strengths Daniels (1945). The stability analysis of this system is quite straightforward (see Bažant and Cedolin, 1991): under a small enough prescribed displacement $U$ of the entire system (interface $u$ plus elastic body $v$), the system has a unique solution, i.e.

$$ \begin{aligned} u &= [(1+k) - \sqrt{(1+k)^2 - 4kU}]/2 \\ v &= u(1-u)/k \end{aligned} \tag{4} $$

The second-order work of the system is:

$$ \delta^2 \mathcal{W} = \frac{1}{2}(k+K)\delta u^2 \tag{5} $$

where $K$ is the tangent stiffness of the bundle, i.e. $dF/du$. $k$ is always positive, as $K$ is positive in the prepeak part of the bundle behaviour, and negative in the postpeak part. The state of the system is stable if $\delta^2 \mathcal{W} > 0$, that is equivalent to $K > -k$.

At the maximum displacement $U^* = (1+k)^2/4k$, we find a bifurcation point, where the solution is no longer unique. At this point, the displacement is $u^* = (1+k)/2$ and the tangent stiffness in the bundle is

$$ \frac{dF}{du} = -k \tag{6} $$

It is exactly opposite to that of the elastic body.
<figure><img src="image_5.png" /><figcaption>Fig. 4. The response of the continuous system, with the stable path and the unstable one.</figcaption></figure>

It corresponds also to the loss of stability, because the second-order work is zero. This is a trivial example of a localization point. If one tries to increase the displacement past $U^*$, then the entire interface fails catastrophically. Under an idealized controlled displacement, the system response follows a snap-back branch, that represents instability ($\delta^2 W' < 0$). Note that prior to the critical equilibrium, strain softening develops in a stable fashion as also pointed out in Bažant's analysis (Fig. 4).

Let us now come back to the finite size fibre bundle as the interface. The response of the bundle is no longer a differentiable law such as eqn (3), but rather a sequence of linear elastic responses limited by end-points where a fibre breaks. Because of randomness of the fibre strength, the response of the system is always unique, and no bifurcation point could be defined. Let us call {$u_i, F_i$} the sequence of failure displacements and forces, respectively, where i indicates the number of broken fibres. In a first approach, we can repeat the same analysis as previously. One cannot of course consider that eqn (6) holds at bifurcation as in the previous example because the response of the bundle is no longer differentiable. Nevertheless, tracing the maximum of the function $F(u) + ku$ yields a criterion for the onset of unstable behaviour, and if $F(u)$ is differentiable one recovers the previous analysis.

As we saw before, the response of such a system is just a succession of fluctuations. It can be analysed through 'avalanches'. With the introduced variables, we can define an avalanche in the fibre bundle as follows: an avalanche of size $\Delta$ and direction $k$, starting at {$u_i, F_i$}, is such that:

$$ \begin{cases} F_i' < F_i - k(u_i' - u_i) & \text{for } i < i' < i + \Delta \\ F_i' \ge F_i - k(u_i' - u_i) & \text{for } i' = i + \Delta \end{cases} \quad (7) $$
<figure><img src="image_6.png" /><figcaption>Fig. 5. A scale of the {$u, F$} response of a discrete model. Each point is the rupture of one fibre. The arrow represents the beginning of an avalanche of size $\Delta$, where $\Delta$ is the number of fibres that break under the arrow. For this example, it is an eight-size avalanche, represented by the blank points.</figcaption></figure>

This means that if the system is loaded up to a point where the ith fibre fails, a series of $\Delta$ fibres will fail simultaneously for the same overall displacement $U$ (Fig. 5). These avalanches are expected to be experimentally observable through e.g. acoustic emissions.

In the continuous case, prior to the point $u^*$, damage in the bundle is controlled and we can say that the avalanche size is 0. At $u = u^*$, the critical equilibrium state is reached and a single avalanche of size equal to the remaining number of bonds in the bundle is observed. In the thermodynamic limit, $N \to \infty$, this avalanche has a size which diverges to infinity. We observe that avalanches do reproduce the result of a standard stability analysis, with a simple 'binary' ($0-\infty$) avalanche size distribution. A crucial point is that this analysis is not entirely correct, in the sense that we have first considered the continuum limit for the force displacement curve, and analysed the avalanches on this mean response. Most of the information which can be derived from the concept of avalanches has been lost in this procedure.

Taking into account the full random and discrete nature of the model, Hemmer and Hansen (Hemmer and Hansen, 1992; Hansen and Hemmer, 1994) succeeded in determining the analytical solution of the probability distribution $n$ of observing an avalanche of size $\Delta$, for any stiffness $k$, starting at any prescribed displacement $u$. They found:

$$ \frac{n(\Delta, u, k)}{N} = \Delta^{-3/2} \Phi\left(\frac{\Delta}{\Delta^*}\right) \qquad (8) $$

with

$$ \Delta^* \propto (u^*(k) - u)^{-2} \qquad (9) $$

where $\Phi(x)$ is a scaling function which is constant for small arguments $x \ll 1$, and drops to zero rapidly for $x > 1$. $\Delta^*$ is the maximum avalanche size. Moreover, the exponent $-3/2$ and $-2$ appearing in eqn (8) are universal, in the sense that they do not depend on the chosen distribution
<figure><img src="image_7.png" /><figcaption>Fig. 6. The log-log graph of the probability distribution n of avalanche of size Δ at a given displacement u, with $u \ll u^*$ [eqn (8)].</figcaption></figure>

of threshold force p. In simple terms, they showed that n is a power-law distribution with an exponent $-3/2$, truncated at a maximum avalanche size, $Δ^*$, which diverges as one approaches the displacement $u^*(k)$ (Fig. 6).

Using the expression of n, one deduces easily that the mean avalanche size behaves as the upper distribution cut-off $Δ^*$, i.e.

$$ \langle \Delta \rangle \propto (u^*(k) - u)^{-2} \qquad (10) $$

Thus, Hemmer and Hansen (1992) showed that the presence of the small amplitude fluctuations around the mean behaviour induced a well-defined statistical distribution of avalanches instead of having the binary variation of the avalanche sizes (Fig. 7).

<figure><img src="image_8.png" /><figcaption>Fig. 7. The evolution of the avalanche mean size (Δ) vs the evolution of the displacement. The line is the binary response of the mean behaviour, as the dotted curve is the discrete response [eqn (10)].</figcaption></figure>
This statistical distribution of avalanches can potentially be used as a precursor of the macroscopic failure point $u^*$. According to Hemmer and Hansen (1992) analysis and in the course of loading the bundle, we would observe a series of avalanches, which can be referred to microinstabilities, which progressively becomes larger and larger, up to the point where they diverge and become ‘macroscopic’. This signals in particular that some care has to be taken when taking the thermodynamic limit. Taking first the continuous limit of the force-displacement response, and analysing its stability, erases the progressive development of the avalanches, and hence misses an important feature of the model.

To demonstrate the utility of such an analysis, we will analyse the fluctuation in the global failure displacement for a finite size system. Let us call $\delta u$ the distance to $u^*(k)$ where the final avalanche is initiated. $\delta u$ can be estimated by writing that the maximum avalanche size at this displacement, $\Delta^*$, allows to increase the displacement $u$ up to $u^*(k)$. Hence

$$ \frac{\Delta^*}{N} \propto \delta u \qquad (11) $$

Using eqn (8), $\Delta^*(\delta u) \propto \delta u^{-2}$, we deduce the scaling

$$ \delta u \propto N^{-1/3} \qquad (12) $$

This estimate gives the typical fluctuation of $\delta u$, i.e. $\langle \delta u^2 \rangle^{1/2}$ as $N$ varies (Fig. 8). It seems difficult if not impossible to derive this result without considering the notion of avalanches.

Finally, let us also note as a side-remark that in the continuous limit, the response of the system is continuous but not differentiable. In this case, one can show that locally, the fluctuating part of the force-displacement response becomes self-affine with a Hurst exponent of 1/2. Thus it belongs to the realm of $C_{1/2}$ functions, rather than $C_1$ as would be needed to apply a criterion such as

<figure><img src="image_9.png" /><figcaption>Fig. 8. The scaling form of the fluctuation. The points are numerical results obtained for different sizes $N$, as the line $y = -2x/3$ is a guide for the eyes.</figcaption></figure>
$dF/du = -k$. This lack of regularity of the response can be seen as the origin of the statistical avalanche distribution.

Because the displacement is constrained to be the same in each fibre, the response tends to a well-defined behaviour, and thus finite size effects play only a marginal role in the present example. For instance, the peak stress, $F$, tends to well defined value with fluctuations of order $1/\sqrt{N}$. However, other systems display a much more significant size effect, which can be tracked back to the cumulative effect of the avalanches. Hence, in this case, the concept of avalanches and of their statistical distribution is unavoidable.

A principal flaw of this simple model is the load redistribution when a fibre breaks: the external force is shared equally between all surviving fibres. On the other hand, when a micro-crack appears in a quasi-brittle material, the stress is redistributed mainly around it, and the local interaction decreases fast as the distance with the micro-crack increases (typically as a power law of the distance). To take into account this effect, we need to introduce this redistribution in the model. Note that a simple local redistribution, where the load of a failed fibre is redistributed equally on the nearest surviving fibres, has been already studied (see e.g. Harlow and Phoenix, 1991).

### 3. Continuous interface model with redistribution

#### 3.1. *Presentation and properties*

For the sake of simplicity, we are going to focus on the situation where the band of strain softening material is small with respect to the size of the elastic block modelled above by a spring and a rigid bar. Hence, we will deal with a softening interface embedded in between a rigid and an elastic substrate. This description is suited to adhesion, and can also be seen as a simplification of a 2-D medium since the redistribution process will be constrained to develop in the direction of the interface. Nevertheless, this example contains the basic features involved in the transitional behaviour between diffuse and localised cracking. At variance with the previously discussed case, we would like to incorporate an elastic coupling between the fibres as mediated directly by the elastic body, i.e. without the rigid bar which redistributed equally the displacement among the surviving fibres. In continuum mechanics, this effect could be represented by the elastic Green function of a semi-infinite plane. For the convenience of the analysis, we resort to a different choice based on a hierarchical decomposition of the elastic body. The structure is the following one: a block is split up in three sub-blocks (Fig. 9). The two lower blocks are then subdivided in three, and this recursively down to the lowest level chosen in the discretization. Each block is described using an elastic uniaxial behaviour. Thus, the elastic body can be seen as composed of springs, connected in parallel and alternatively in series. At the lowest level, each finer block is connected

<figure><img src="image_10.png" /><figcaption>Fig. 9. The two first decompositions of an elastic block.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10. A hierarchical system of generation 3, with 10 fibres for each bundle.</figcaption></figure>

to a fibre bundle of small size. These bundles represent the interface. A system of generation one (number of level in the hierarchical decomposition) is thus similar to the previous model.

Figure 10 shows an example of a 3-generation model. With 10 fibres for each bundle, a system of generation 12 is for instance made up of $10 \times 2^{12-1} = 20,480$ fibres. The simplicity of the construction allows for the numerical simulation of extremely large sizes, while still preserving the long range nature of the elastic couplings. In two dimensions, all springs have the same stiffness $k$ at all generations, but their initial lengths is divided by two as the generation is decreased by one. At the first generation, however, the aspect ratio of the element is twice that of all other generations, and thus the first springs in contact with the interface have a stiffness $k/2$. This allows to obtain a global stiffness for the entire elastic system which is independent of the discretization level as expected (see the Appendix). An important point which will be used later concerns the interaction between an intermediate level (say index $i$), with the rest of the medium. We can compute recursively the stiffness $\Lambda_i$ of this structure deprived from one subblock $i$, if a force is applied at this level. We find that $\Lambda_i$ is exactly equal to $k$, i.e. just as if the subsystem was simply connected to the exterior world by a single block (see the Appendix).

As the generation increases it can be shown that the elastic coupling will be the same as for a continuous model, i.e. with an influence function scaling with the same power-law as the Green function in an elastic continuum. In order to reach this result, one needs to introduce a distance suited to our discretization. Considering two points along the interface, we search from the smallest block which contains both points. If $j$ is the block generation (which could take value between 1 and $N$, which is the generation of the entire system), the distance is then defined as

$$d = 2^{j-1} - 1 \qquad (13)$$

This distance has the special feature of being ultrametric (i.e. the same properties than an usual distance, except for the triangular inequality where $d(A, C) \le \min(d(A, B), d(B, C))$. With this definition, one can show that under an applied force $F$ on a point of the interface, the induced displacement $v$ of an other point is

$$v(j) = (N - (j-1)) \frac{F}{k} + v_0 \qquad (14)$$
where $v_0$ is the initial displacement of the considered point, and $N-(j-1)$ is nothing but the number of springs that separate the two points. By introducing the distance $d$, and for large $j$, it leads to:

$$v(d) = A \log \left(\frac{B}{d}\right) \frac{F}{k} + v_0 \qquad (15)$$

where $A = 1/\log 2$ and $B = 2^N$ are constants. It is exactly the same form that the Green function of a semi-infinite plane. The only difference is that this influence function consists in constant plateaus whose size increases in geometric series. This is a residual effect of the two-fold splitting of each level.

## 3.2. Continuous interface

We now proceed by considering first the case where each fibre bundle is changed into a damageable element, with a behaviour law derived from the mean fibre-bundle response, eqn (3). We use the hierarchical construction to relate the interface law to the global (interface plus elastic body) response.

Let us construct a system at generation $(n+1)$, starting from two generation-$n$ subsystems. These last subsystems are supposed to be described by two force-displacement relations $F_1^{(n)}(U^{(n)})$ and $F_2^{(n)}(U^{(n)})$. We wish to find the global $F^{(n+1)}(U^{(n+1)})$ response. The two subsystems are subjected to the same displacement, hence the force is

$$F^{(n+1)} = F_1^{(n)}(U^{(n)}) + F_2^{(n)}(U^{(n)}) \qquad (16)$$

The same force also stretches a spring of stiffness $k$ in series with the blocks, and thus the global displacement is

$$U^{(n+1)} = U^{(n)} + \frac{F_1^{(n)}(U^{(n)}) + F_2^{(n)}(U^{(n)})}{k} \qquad (17)$$

These equations provide a parametric representation of the $(n+1)$th generation system as a function of the $n$th generation.

Let us first assume that the interface is homogeneous, thus in the previous analysis $F_1 = F_2$. Hence

$$U^{(n+1)}(F) = U^{(n)}(F/2) + F/k = U^{(1)}(F/2^n) + (1-2^{-n})F/k \qquad (18)$$

The final equation simply relates the interface displacement to the global one $u(F) = U^{(1)}(F)$. We note that as $n$ increases, the global behaviour is nothing but that of the elastic medium because the displacement in the interface represents a vanishing contribution. We can invert the previous relation to obtain the homogeneous interface law from the global response. This is a practical tool to compute the equivalent homogeneous interface law when some inhomogeneity exists locally.

From the previous equation and because the interface response is continuous, the tangent (subscript tg) and secant (subscript sc) stiffnesses of the entire system at generation $n$ can be computed:
$$ 
\begin{aligned}
K_{\text{tg}}^{(n)} &= \left( \frac{1 - 2^{1-n}}{k} + \frac{2^{1-n}}{1 - 2u} \right)^{-1} \\
K_{\text{sc}}^{(n)} &= \left( \frac{1 - 2^{1-n}}{k} + \frac{2^{1-n}}{1 - u} \right)^{-1}
\end{aligned} 
\qquad (19) $$

where we have used the local interface displacement $u$ to characterize the loading, $F = 2^{n-1}u(1-u)$, assuming in this formula a homogeneous displacement all along the interface.

### 3.3. Bifurcation analysis

For moderate displacements, and in the absence of randomness in the interface, every point of the interface undergoes the same damage. The system response however may cease to be unique at a particular displacement for which the ‘global’ displacement $U^{(n)}$ is maximum.

Let us first assume that one half of the interface is subjected to an increasing damage while the other half is elastically unloaded. Using eqn (17), we see that this bifurcation condition is reached when

$$ k + K_{\text{sc}}^{(n-1)} + K_{\text{tg}}^{(n-1)} = 0 \qquad (20) $$

From the expression for the secant and tangent stiffnesses eqn (19), we obtain an equation for the displacement, $u = u_1^*$ at the interface level for which a first bifurcation is encountered:

$$ (3 - 2\epsilon)(1 - 2\epsilon)(1 - u)(1 - 2u) + 2k\epsilon(2 - 2\epsilon)(2 - 3u) + 4k^2\epsilon^2 = 0 \qquad (21) $$

where $\epsilon = 2^{1-n}$. Focusing on the large system size limit, we can expand the solution in order of $\epsilon$ and obtain the solution as

$$ u_1^* = \frac{1}{2} + \frac{2k}{3l_1} + \left( \frac{10}{9k} - \frac{4}{9} \right) \frac{k^2}{l_1^2} + \mathcal{O}(l_1^{-3}) \qquad (22) $$

where we have introduced the size of the interface where the damage localises, $l_1 = 2^{n-1} = L/2$, to express the result in physical terms. $L$ refers here to the number of damageable elements (or fibre bundles in the discrete case).

At this stage, there is a bifurcation to three possible evolutions: either damage remains inhomogeneous (but this solution is unstable) or only one of the two subsystems continues to be damaged while the other is elastically unloaded. Due to the symmetry of the system these two solutions are identical.

For a large system size, $\epsilon \to 0$, we note that $u_1^*$ tends to $1/2$, i.e. the interface displacement at peak force. Bifurcation is however delayed to a larger displacement by a quantity proportional to $(k/L)$. The homogeneity in the latter expression can be restored if we consider that the stiffness of the interface fibres is not unity, and that the interface is in fact a band of softening material of width $h$. The offset of the first bifurcation point is then of order

$$ \frac{u_1^* - u_{\text{peak}}}{u_{\text{peak}}} \propto \frac{E_{\text{bulk}}(1 - v_i - 2v_i^2) h}{E_i(1 - v_i)} \qquad (23) $$
where $E_{\text{bulk}}$ is the Young modulus of the elastic block, and $L$ its size, $E_i$ and $v_i$ are the Young modulus and Poisson ratio of the band of width $h$. The occurrence of $v_i$ comes from the antiplane displacement in the layer.

In this analysis, we have postulated that the first bifurcation mode appeared at the macroscopic scale. One can perform the same computation for any intermediate level $1 \le i \le n$, keeping the boundary condition on $U^{(n)}$. The only variance with eqn (20) is that the stiffness $k$ has now to incorporate all the intermediate levels from $i$ to $n$. The hierarchical structure allows to compute this stiffness which remains simply equal to $k$ at all levels. Therefore, the localization at generation $i$ appears for a displacement $u_{n-i}$ given by eqn (22) where $l_{n-i} = 2^i = L/2^{n-i}$ is to be substituted to $l_1$.

$$u_{n-i}^* = \frac{1}{2} + \frac{2k}{3l_{n-i}} + \left(\frac{10}{9k} - \frac{4}{9}\right) \frac{k^2}{l_{n-i}^2} + \mathcal{O}(l_{n-i}^{-3}) \quad (24)$$

Thus, these modes will occur much later than the first one $l_1 = 2^{n-1}$. They will however be of interest if we proceed along one of the two symmetric stable branches past the first bifurcation point. The $n$th generation subsystem where the damage continues to progress will encounter a bifurcation point similar to the previous for $u = u_2^*$. Past this local displacement, the damage concentrates on one quarter of the system while the rest will be elastically unloaded. The same analysis can be carried out to any stage down the cascade of bifurcation always concentrating on a stable branch. The local displacement of the interface on the active part of the interface at the $i$th bifurcation is given by eqn (24).

We thus obtain a simple physical picture of the post-localisation regime (localisation is understood here as bifurcation) where the damage zone progressively condenses onto a smaller and smaller ('active') region, while the rest of the structure is elastically unloaded.

### 3.4. Post bifurcation response

An important feature which deserves a particular interest is the equivalent interface law which can be measured past the first bifurcation point. Indeed, as soon as the damage is no longer homogeneously distributed, the equivalent homogeneous law is no longer similar to the one of any of the constituents. However, if we were to perform the experiment, without any a priori knowledge of the cascade of bifurcations, the equivalent homogeneous law is the one we would extract from the load-displacement curve.

The easiest way to have access to such an equivalent law for a system at generation $n$ is to use the hierarchical nature of the decomposition of the elastic block. Let us assume that we know this equivalent law for a system of generation $n-1$, and express the equivalent law at the next generation. As discussed above, the first bifurcation point occurs first at the largest scale. Past this first point, one half of the lattice is simply elastically unloaded. The other half is described by the homogeneous equivalent law. Let $(1/2+x_1^{(n)}, 1/4-y_1^{(n)})$ be the displacement-force coordinate of the first bifurcation point. As observed above, we know that the $x_1^{(n)}$ form a geometric sequence of ratio $1/2$:

$$x_1^{(n)} = x_1^{(1)} 2^{-n}. \quad (25)$$
For the y coordinate, it suffices to observe that the first bifurcation point lies on the homogeneous characteristic, and hence $1/4 - y_1^{(n)} = (1/2 + x_1^{(n)})(1/2 - x_1^{(n)})$ where we have used the specific parabolic form of the ‘bare’ interface law. Thus,

$$y_1^{(n)} = (x_1^{(1)})^2 4^{-n}. \qquad (26)$$

Any point $(1/2+x, 1/4-y)$ of the $(n-1)$ generation equivalent homogeneous interface law is transformed into $(1/2+x', 1/4-y')$ such that

$$x' - x = \frac{2^n}{k}(y' - y)$$

$$y' - y = \frac{\frac{1}{2} + x + \frac{2^n}{k}\left(\frac{1}{4} - y\right)}{\frac{1}{2} + x_1^{(n)} + \frac{2^n}{k}\left(\frac{1}{4} + y_1^{(n)}\right)} \qquad (27)$$

We label the succession of bifurcation points by a subscript $j$ whereas the superscript $n$ refers to the system generation. Solving for the asymptotic (large $n$) behaviour of the $(x, y)$ variables provides

$$x_j^{(n)} = x_1^{(1)} 2^{-n} j$$

$$y_j^{(n)} = (x_1^{(1)})^2 4^{-n} \left( \left( \frac{1}{4} x_1^{(1)} - \frac{1}{18} k \right) 4^j + \frac{2}{3} k j - \frac{4}{9} k \right) \qquad (28)$$

Figure (11) shows the sequence of bifurcation points for $n = 7, n = 9$ and $n = 11$ computed exactly, together with the asymptotic expression shown as a curve. We observe that only the latest $j \approx n$ points are not well described by the asymptotic behaviour. However, as $n$ increases, most of the cascade is very accurately described.

To make the above result more explicit, we note that past the first bifurcation point, the force-displacement relation becomes size-dependent, but it can be cast in a simple scaling form using the system size $L = 2^n$:

$$\left(\frac{y}{y_1^{(1)}}\right) = L^{-2} \Psi\left(L \frac{x}{x_1^{(1)}}\right) \qquad (29)$$

where the scaling function is just:

$$\Psi(z) = \left(\frac{1}{4}x_1^{(1)} - \frac{1}{18}k\right)4^z + \frac{2}{3}kz - \frac{4}{9}k \qquad (30)$$

One important point to be noted here is the difference of exponents of $L$ which appears in $x$ and $y$. As a consequence, the equivalent homogeneous interface law shows a sudden decay of the force at constant displacement past the peak force.

It is also of interest to consider the scaling of other physical quantities. In particular, if we come back to the picture of the fibre bundle at the interface level, we may introduce another variable
<figure><img src="image_12.png" /><figcaption>Fig. 11. The sequence of the bifurcation points $(x_j^n, y_j^n)$. The continuous curve shows the asymptotic behaviour obtained from the recurrence relation. The dotted curves are the real sequence of bifurcation points for a 7-generation (crosses), 9-generation (black point) and 11-generation (triangle) system. Note that just the first points are well-described by this relation, and the snap-back part is not represented.</figcaption></figure>

which is the number *N* of broken fibres. The latter is simply related to the displacement as $N = 2^n u$. Therefore, we conclude that two consecutive bifurcations are separated by a fixed number of broken fibres. In fact the displacement in the active region for the consecutive bifurcations increases exponentially fast, as $2^j$, but simultaneously the active region shrinks also exponentially, as $2^{-j}$, so that the product of these two terms which gives the number of failed fibres remains constant.

It is now a simple matter to express the variations of *x* or *y* as a function of the number of broken fibres (past the peak force, where $N = N_p = L/2$): *x* is linear in $(N - N_p)$, whereas *y* grows exponentially fast. Simultaneously, the size of the active region is $L/2^j$, decreasing exponentially fast with *j* or equivalently $(N - N_p)$ or $x_j$.

### 3.5. From damage localization to crack nucleation

The physical picture which arises from this analytic solution is of particular interest, since it is one of the rare situations where some insight can be obtained past the first bifurcation.

We have seen that the interface degradation process consists in a progressive condensation of the damaging region from the structure scale down to the basic constitutive unit. At the end of this first cascade, exactly one of the smallest size interface element is totally broken. This naturally forms the initiation stage for a crack propagation regime. Unfortunately, following the crack propagation is of little interest in our model, which is then too much sensitive to the detail of the hierarchical decomposition to pretend any possible comparison with reality. However, up to the crack nucleation, we believe that the hierarchical interface model is a faithful description of continuum model, yet simple enough to be amenable to an analytic solution.
An interesting feature is to be noted at the crack nucleation stage: the progressive condensation of the damage can be read back from the damage profile along the interface. Indeed, the damage $D$ of the interface is a simple linear function of the maximum displacement ever encountered by a homogeneous domain, that varies between 0 and 1. Tracing backward the damage in the active region, we obtain that the damage $D$ at a distance $d$ from the crack nucleation point is

$$D(d) = \frac{1}{2} + \frac{k}{2d} \qquad (31)$$

where we have used the earlier defined distance [eqn (13)].

Thus, the cascade of bifurcation leads to a rather unusual damage profile ahead of the crack. If we define a 'process zone' as a damage zone ahead of a crack, we would conclude that the process zone is of infinite extent. However, this is important to note that the damage decreases very fast with the distance, i.e. as an inverse law. Then this process zone seems more similar than the classical one, that is a quasi-confined damage zone of finite length ahead of a crack.

## 4. Discrete interface model with redistribution

We have already underlined the importance of the notion of avalanches for a disordered fibre bundle. In the interface model, basically the results of Hemmer and Hansen still hold. The same statistics is expected in this case. The only variant comes from the boundary conditions. We have seen that an elastic coupling to the fibre bundle has the major effect of moving the interface displacement [defined in eqn (4)] at which the avalanche size diverges. In the interface case, the elastic coupling is a little more complex, and thus the point of divergence for avalanches requires some discussion.

Let us consider a block at generation $n$. This block is subjected to an imposed displacement through a device of stiffness $\kappa^{(n)}$. The avalanches which are meaningful at this level are those which are constructed from the global force-displacement characteristic at generation $n$ with a slope $-\kappa$. We would like to relate those avalanches to the one computed at the previous generation. We have seen above how to relate the force-displacement relations from one generation to the next. This provides a simple equivalent stiffness of the loading device $\kappa^{(n-1)}$ to be considered at the $(n-1)$th generation.

$$\kappa^{(n-1)} \equiv H(\kappa^{(n)}) = \frac{\kappa^{(n)}}{2(1 - \kappa^{(n)}/k)} \qquad (32)$$

where the function $H$ is shown in Fig. 12. Iterating the previous transformation allows to compute the elastic coupling to be considered directly at the interface level. The function $H$ has two fixed points, $\kappa = 0$ and $\kappa = k/2$. The first one is attractive, whereas the second one is repulsive.

In order to better understand what is the practical measuring of the slope, let us consider a system consisting in a few elements. If we are looking for the first bifurcation point, that is equivalent to the divergence of the avalanche sizes, we have to consider avalanches with a stiffness equal to $\kappa = k/2$ for just one bundle response. Because we are far from the fixed points, we use the relation eqn (32) to obtain this slope. Figure 13 illustrates this point for an 11-generation system.
<figure><img src="image_13.png" /><figcaption>Fig. 12. The function $H$, for $k = 0.5$. The thin line is the equation $y = x$, so the two intersections are the fixed points, $x = 0$ and $x = k/2$.</figcaption></figure>

Therefore, for most values of $\kappa^n$, the equivalent stiffness to be considered at the interface level $\kappa^{(1)}$ tends to 0 as the system size tends to infinity. This means in practice that for most boundary conditions, the avalanches should be analysed at the interface level with an elastic coupling which tends to 0, i.e. under a constant force condition. This is precisely what has been shown in the previous analysis, where we considered $\kappa^{(n)} \to \infty$, i.e. a constant displacement imposed on the entire elastic domain, and we have retrieved that the first bifurcation occurred for a displacement at the interface level which approached the apex of the force-displacement curve ($u = 1/2$). The slight delay in this displacement resulted from the last iterations of the function $F$. Indeed, for $\kappa \to 0$, $H(\kappa) \approx \kappa/2$, and thus, $\kappa^{(i-1)} \approx \kappa^{(i)}/2$ for $i \ll n$. We observed that the first bifurcation in a generation $n$ system occurred at points $u^*(1) = 1/2 + B2^{-n}$ where $B$ is a constant, and thus the tangent stiffness $du/dF(u = u^*(1)) = B2^{-n}$ is indeed a geometric series of ratio $1/2$.

The existence of the unstable fixed point $\kappa = k/2$ can also easily be understood: if we invert the relation eqn (32), we can relate the larger scale stiffness to the lower one, through the function inverse of $H$. In this case the fixed points remain obviously identical but their attractive or repulsive character is turned to the opposite. This means that the stiffness of the entire system tends to $k/2$ as $n$ increases to infinity. The $k/2$ is nothing but the stiffness of the elastic body computed in the preceding section.

This shows that the conditions for bifurcation become independent of the global boundary conditions as the system size diverge. It also underlines the fact that in order to observe the cascade of bifurcation, one should use an active control on the loading conditions, with the ability to decrease the loading fast compared to the typical time needed to fracture the fibres, or to redistribute the load to the fibres. This imposes some severe constraints on the monitoring of the experiment. A possible way to build this control might be to use the acoustic emission during loading.

Let us note that the notion of avalanche allows to understand naturally the cascading process.
<figure><img src="image_14.png" /><figcaption>Fig. 13. The interface response of a 11-generation model. The thin curve is the parabola, and the line corresponds to the slope of the first bifurcation point. The dotted curve is the response of the equivalent homogeneous system. Note the good agreement with the bifurcation of the interface response from the parabola.</figcaption></figure>

Indeed, we have seen that for a subblock embedded inside the entire structure the effective stiffness of the surrounding medium amounts to *k* (instead of *k*/2 if all other subblocks of the same generation are subjected to the same displacement). Therefore, at the bifurcation point where the damage is localized in a subblock of generation *i*, the two subblocks at generation *i*−1 are still
stable, i.e. the maximum avalanche size in each of these two blocks is finite. Hence the damage
will be shared between the two subblocks up to the next bifurcation point.

This argument also allows to estimate the validity of the cascade once the fluctuations due to
the random nature of the fibre bundles are taken into account. As the size of the active region, $l$, decreases, the force fluctuation increases as $l^{-1/2}$, and the proportion of broken fibres displays a fluctuation of order $l^{-1/3}$. Comparisons of the level of fluctuations with the increment of force, displacement or number of broken bonds, show that the late stage of the process ($l$ small enough) are dominated by the fluctuations, but in contrast, the early stage is well defined. Therefore, we anticipate that the first steps of the cascade may be correctly described by the above homogeneous situation, whereas the more mature stage may be scrambled by the presence of disorder. A representation of the location of the fibres that break under the loading gives a good physical idea of the cascade phenomena (Fig. 14).

The bifurcation cascade observed during the failure is not the usual idea of the failure of a joint:
such a failure generally occurs catastrophically, and then the first idea is to think that it is due to
a critical flaw. It is important to note that our model follows the same catastrophic behaviour if
we consider the global load-displacement response. Hence, Fig. 15 shows the global response of
the model for three different generation systems. As the generation increases, the behaviour
becomes more and more elastic brittle, as expected for a joint failure. But if we consider just the
interface response, we find effectively the previous behaviour with the bifurcation cascade.

Finally, in spite of their different description, we see that the heterogeneous discrete system (for
large enough sizes) and the homogeneous one have the same post-peak behaviour:

* The relation *F* vs *u* is similar until the apparition of the first crack.

• The onset of localisation appears at the same time (Fig. 13).

* The damage cascade is observed in both cases.

**5. Conclusion**

For continuous models, the localization is well-defined, using some criteria based on the loss of
uniqueness or the study of the tangent stiffness operator. For discrete models, the localization
could not be defined in such a manner. The solution is always unique, and no tangent could be
calculated on the response because of the fluctuations that are superimposed.

In some well-known cases, the loss of uniqueness coincides with the loss of stability, where a
bifurcation point is encountered. In the first part, we show that the study of avalanche statistics
allows to detect this point. Precisely, the divergence of the avalanche sizes could be directly
compared to the loss of stability in a continuous model.

After defining this equivalence, we propose, as an application, to study a damage interface
coupled with an elastic block. For the sake of simplicity, the interface is chosen to be thin, then
the damage propagates only in the interface direction. We are interested particularly in the unstable
path, that is very difficult to observe with continuous models. The discrete model that we use is a
hierarchical model, that has a good representation of the Green influence function in an elastic
continuum. Our conclusions are the following ones:
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. The representation of the condensation of the broken fibres for a 12-generation system, with 10-fibre bundle. The x-axis represents the fibre location (from 1–20,480), and the y-axis is the succession of the broken fibres under the loading. A scale is made from the first bifurcation to the initiation of the first crack.</figcaption></figure>

* The divergence of the avalanche sizes coincides effectively with the first bifurcation point on a continuous model.

* The damage in the interface is first homogeneous, and then condensate progressively into a
<figure><img src="image_17.png" /><figcaption>Fig. 15. The global force-displacement response for a 5-generation (dotted line), 7-generation (dashed line) and 9-generation (continuous line) system. The behaviour becomes elastic brittle as the generation increases.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. The variables of a hierarchical structure.</figcaption></figure>

narrow region, down to a single point. Then a crack is initiated and propagates up to the total failure.

* The bifurcation points that lead to the crack formation are well defined, and are separated by a constant number of broken bonds.

* Just before the apparition of the crack, the damage profile is obtained as an inverse power law, and is spreading over all the interface.
<figure><img src="image_19.png" /><figcaption>Fig. 17. The representation of $\Lambda_i$.</figcaption></figure>

* The post-peak response could be accounted for by unusual scaling law (29) with affine transforms on force-displacement relations.

Finally, for heterogeneous systems, we would like to emphasize the usefulness of the avalanche statistics, as a practical tool to cope with non-smooth responses. We saw in the example studied in this article, that it constitutes a complementary approach to standard stability analyses and provide additional information, in particular through precursors of global loss of stability. This feature is certainly worth being further studied experimentally.

## Appendix: Value of stiffness in the hierarchical model

We first propose to establish the recurrence relation between the stiffness of a *i*-generation structure, $K_i$, and a (*i*+1)-generation one, $K_{i+1}$.

Per definition of the stiffness, the external force $F$ is:

$$F = K_{i+1} U_{i+1} \qquad (33)$$

We search the expression of $K_{i+1}$, as a function of $K_i$ and $k$. Using the hierarchical structure of the block, we can write

$$F = 2K_i U_i = F = k(U_{i+1} - U_i) \qquad (34)$$

and thus,

$$F = k \left( U_{i+1} - \frac{F}{2K_i} \right) = \frac{2kK_i}{k+K_i} U_{i+1} \qquad (35)$$

By identifying with eqn (33), we obtain the recurrence relation:
$$K_{i+1} = \frac{2kK_i}{k + 2K_i} \qquad (36)$$

The stable fixed point of this recurrence relation is

$$K_i = \frac{k}{2} \qquad (37)$$

Then, choosing the value $K_1 = k/2$ (i.e. the stiffness of the first springs in contact with the interface), the global stiffness of the hierarchical structure is independent of the discretization level as expected.

We now give the stiffness $\Lambda_i$ of a $n$-generation structure deprived from one subblock $i$, if a force is applied at this level. Again using the hierarchical structure of the block leads to a simple recurrence:

$$\Lambda_i = \frac{\Lambda_{i+1}k}{\Lambda_{i+1} + k} + \frac{k}{2} \qquad (38)$$

The stable fixed point is thus

$$\Lambda_i = k \qquad (39)$$

## References

* Bak, P., 1996. How Nature Works: The Science of Self-Organized Criticality. Springer-Verlag, New York.

* Bažant, Z.P., Cedolin, L., 1991. Stability of Structures. Oxford University Press.

* Daniels, H.E., 1945. The statistical theory of the strength of bundles of threads. Proceedings of the Royal Society, London. A183, 405-435, 1945.

* Delaplace, A., Pijaudier-Cabot, G., Roux, S., Progressive damage in discrete models and consequences on continuum modelling. J. Mech. Phys. Solids, 44, 99-136.

* Dudzinski, D., Molinari, A., 1991. Perturbation analysis of thermoviscosplastic instabilities in biaxial loading. Int. J. Solids Structures, 27, 601-628.

* Fokwa, D., 1992. Matériaux hétérogènes: analyse expérimentale et modélisation numérique par une approche hiérarchique. Thèse de doctorat de l'université de Paris 6, LMT.

* Hansen, A., Hemmer, P.C., 1994. Burst avalanches in bundles of fibers: local versus global load-shearing. Physics Letters A 184, 394-396.

* Harlow, D.G., Phoenix, S.L., 1991. Approximations for the strength distribution and size effect in an idealized lattice model of material breakdown. J. Mech. Phys. Solids, 39, 173-200.

* Hemmer, P.C., Hansen, A., 1992. The distribution of simultaneous fiber failures in fiber bundles. Journal of Applied Mechanics, 59, 909-914.

* Hill, R., 1959. Some basic principles in the mechanics of solids without a natural time. J. Mech. Phys. Solids, 7, 209-225.

* Krajcinovic, D., Lemaitre, J., 1987. Continuum damage mechanics theory and applications. CISM Course, Springer-Verlag, Berlin.

* Laws, N., Brockenbrough, J.R., 1987. The effect of microcrack systems on the loss of stiffness of brittle solids. Int. J. Solids. Struct., pp. 1247-1268.

* Lemaitre, J., 1992. A Course on Damage Mechanics. Springer.

* Leroy, Y., 1991. Linear stability analysis of rate-dependent discrete systems. Int. J. Solids Structures, 27, 783-808.

* Paczuski, M., Maslov, S., Bak, P., 1995. Avalanche dynamics in evolution, growth and depinning models. Phys. Rev. E, 53, 414.