Avalanche Statistics of Interface Crack Propagation in Fiber Bundle Model : Characterisation of Cohesive Crack

A BSTRACT : This paper considers a model of crack propagation taking place at the interface between a rigid support and an elastic plate. The interface is modeled using a ﬁber bundle model (i.e., describing a damage behavior using a discrete set of elastic brittle elements having a random strength). This paper studies the ﬂuc-tuations of the force required to propagate the crack along the interface. The statistics of avalanches, deﬁned as a series of elements that are broken simultaneously under a load that decreases with the crack advance, are studied numerically and analytically. Local ﬁber breakage kinetics is related to a correlation length, which sets the size of a fracture process zone.


INTRODUCTION
Damage in quasi-brittle materials such as concrete is due to the accumulation of microscopic crack nucleation, propagation, and most importantly, arrest.Arrest of microcracks is an essential feature to produce progressive damage; otherwise, a brittle behavior would result.To allow for microcrack arrest, it is therefore important to deal with material heterogeneities.In spite of the importance of this factor, most damage approaches in continuum mechanics treat the solid as homogeneous and little attention is paid to the detailed material microstructure producing the progressive damage behavior.Only the mean effects of the microcrack nucleation, propagation, and arrest are considered.The effects that are neglected (i.e., the intrinsic heterogeneity of the solid) may, however, play a role, and this study illustrates some of these effects in terms of the fluctuations of the macroscopic response due to the existence of an underlying heterogeneity.
An experimental example illustrates this discussion.Fig. 1 shows the measured force-displacement relation for a fourpoint bending test performed with a steel-fiber reinforcedconcrete specimen.Different curves corresponding to different samples of nominally identical material and geometry are shown.One observes a rather large variability of the response.This noise is not an artifact of the measurement device but really reflects the material behavior.It renders difficult the identification of a constitutive behavior law from such a test.Moreover, considering one single curve, one observes a noisy response consisting in abrupt drops of the force at fixed strain.Elucidating the statistical feature of such a ''noise'' is one of the objectives of this work.
Another effect that may result from this heterogeneity is the size effect (i.e., size dependence of the macroscopic strength), a crucial property when scaling of structural response has to be performed.However, in the particular model considered hereafter, such a size effect does not exist.
This paper will mainly focus on the fluctuations of the macroscopic response, in a geometry such that a usual macroscopic modeling will produce a steady response.This particular choice will allow one to study the fluctuations directly and to emphasize the effective noise, which would be obtained experimentally, in the form of ''avalanches'' (i.e., a series of microcrack events that are produced quasi-simultaneously under constant conditions of loading to be specified more precisely below).The notion of avalanches is here borrowed from a large body of works in statistical physics where this notion was proven to be useful in the analysis of various critical systems (Hemmer and Hansen 1992;Paczuski et al. 1995;Delaplace et al. 1999a).Those avalanches depend not only on the total statistical distribution of critical loading but, more importantly, on time (or displacement) correlations.
There are few models that can address such questions without having to rely heavily on numerical simulations.Among these, one of the most intensively studied is the fiber bundle model introduced by Daniels (1945).It consists of a series of elastic brittle ''fibers'' loaded in parallel between two rigid bars.The advantage of this simple model is that almost all of its properties can be solved for analytically.Moreover, it has been extended to account for a wide variety of more complex phenomena, such as time-dependent failure.The fact that the fibers are loaded through two rigid parallel bars implies that all fibers are subjected to the same displacement; thus, the load is always equally shared by all the surviving fibers.This is called the global load sharing rule.The opposite limit of a very compliant interface has been introduced by Harlow and Phoenix (1991) and gives rise to a very different behavior where the system behaves as macroscopically brittle.An intermediate case, taking into account a coupling of the fibers to a semi-infinite elastic medium, has more recently been studied by Delaplace et al. (1999a).This last case gives rise to a homogeneous damage up to a localization point where the most unstable case has a wavelength equal to the system size, provided the elastic block is not too compliant.A major conclusion from this study was that, at the inception of interface cracking, there is no such thing as a finite size fracture process zone.If a process zone of finite size has to be expected, it must develop after the inception of cracking, in the propagation regime that is of concern in this paper.
Studying the original Daniels' model, Hemmer and Hansen (1992) have worked out analytically the statistics of avalanches for this system.This is one of the few works that really extracted information about the fluctuation of the force-displacement characteristics for this system.The same results can immediately be transferred to the case of a more realistic load transfer mechanism, such as the case of an elastic block.The statistics of avalanches has also been extended to the local load sharing rule (limit of extremely compliant coupling) in Kloster et al. (1997).The aim of the present study is to extend the study of these avalanches to a system derived from the fiber bundle model, where fibers are connected to a deformable media and where the geometry is suited to an inhomogeneous, but steady, state of loading.The geometry mimics the propagation of a crack front at the interface between a rigid interface and a deformable plate (or rather, a beam in two dimensions).It can be looked at as a schematic model for mode I crack propagation in a quasi-brittle material (the rigid side of the interface representing an axis of symmetry).The purpose of introducing a nonlinear interface is to concentrate cracking onto a specific region, a line in the present approach, the same as in a cohesive crack model.
In the context of a steady-state crack propagation, this problem has been studied with different approaches, motivated by an experimental study of the roughness of crack fronts by Schmittbuhl and Ma ˚løy (1997) and Delaplace et al. (1999b).Discrete fuse networks have been studied by Zapperi et al. (2000) to mimic this situation, considering the steady-state regime.The problem of nucleation of such interfacial cracks using a variety of load sharing rules has also been developed.However these studies were focused on either the global average behavior or the morphology of the crack front.The focus here will be different in the sense that one concentrates on the characterization of the fluctuations of the force-displacement curve.
This paper is organized as follows.First, the model is defined and a mean-field version of it is introduced in which the profile of imposed displacement along the interface is fixed.The statistical distribution of the loading is obtained simply in this limit and is shown to converge to a Gaussian as the size of the ''process zone,'' as compared to the microstructural size, diverges.Then, ''avalanches'' are introduced and their statistics studied numerically and analytically in the framework of the simplified model.Finally, the numerical study of the complete model is reported, in which the displacement at the interface results from beam deformation.The results are shown to be similar to those obtained for the simplified model.

General Case of Elastic Coupling
First recall the definition of the Daniels' model (1945).A collection of elastic brittle fibers are loaded in parallel between two rigid bars.Each fiber has a perfectly brittle behavior with the same stiffness .However, the strength of the fibers are random uncorrelated variables.The loading is an imposed displacement, and the response is the global force-displacement curve.
In the model, the fibers are connected to a semi-infinite 1D beam on one side and a rigid substratum on the opposite side.One could also consider another elastic beam instead of the rigid part (even having a different stiffness) without any change in the formulation of the problem.A normal displacement is imposed at one point of the elastic beam, which may move along the interface, as if a wedge was pushed in a double cantilever geometry.Fig. 2 illustrates the model.In the notation, the position of fiber i is x i and its extensional displacement is denoted by y i .The spacing between fibers is set to 1 (i.e., x iϩ1 Ϫ x i = 1), thus defining a fixed microstructural size.The maximum extension a fiber can stand before breaking y c is chosen randomly, as mentioned above.This distribution is uniform in the interval [0:1] in the analytical part of this study, whereas other distributions could also be considered in the simulations.A fiber that is strained up to its critical extension is broken instantaneously.The normal displacement imposed on the beam is such that all fibers under this extension are broken with a probability of 1.
The interface consists basically in three parts.Under the wedge, all fibers are broken.This constitutes the crack by itself.Ahead there is a ''process'' zone where a finite fraction of fibers are broken.Finally, moving away from the crack tip, broken fibers become more scarce and their influence can be neglected.There is no clear-cut separation in the last two regions; however, the probability of fiber rupture ahead of the crack tip decays exponentially and thus can be used to define the extension of the process zone.The terminology for the fracture process zone is borrowed from Hillerborg et al. (1976), who studied fracture in concrete with the help of a cohesive crack model.They assumed that the fracture process zone was a region of intense damage with the material, although aggregates bridge the crack and allow for stress transfer across the crack.

Imposed Displacement Profile
A simplified model where the displacement profile along the upper beam is imposed will be implemented in the analytical study.The interface opening has an exponential shape, which captures some of the features of a beam deflection (as will be seen later).For any abscissa x, the profile y is given by where the length scale is considered as a fixed parameter; and U = time-dependent horizontal displacement of the edge.
The major advantage of this simplifying assumption is that it allows derivations of closed-form expressions for the loading while it is not too far away from the exact beam deflection.

SIMPLIFIED ANALYTICAL MODEL Mean Behavior and Fluctuations
The total force exerted on the wedge is simply the sum of all fiber contributions where each individual force f i is In the case of a uniform distribution in the interval [0:1] of the breaking extension of individual fibers, one may compute the average force in fiber i for the ith fiber at a location such that i Ն U (the fiber spacing is set to 1) and ͗ f i ͘ = 0 otherwise.Hence, the mean total force is where d(U ) = U Ϫ int(U ) (i.e., the noninteger part of U); and ␣ ϵ exp(Ϫ1/).Fig. 3 shows the resulting force as a function of U. One observes a small amplitude component due to the discrete nature of the fiber distribution.However, the amplitude of the latter modulation is very small and typically it is negligible when compared to the fluctuation of the force, when the loading length scale is large enough.For integer valued displacement d(U ) = 0, and for large , the mean force takes its minimum value and can be approximated by the following expression [from ( 5 The maximum force is obtained for displacement such that Substituting the last expression in (5) results in the peak mean force and thus relative peak-to-valley difference vanishes as 1/(4 2 ) for large .Therefore in the following, one may disregard this systematic variation and only consider the case of integer displacement U.
The variation of the total force can also be obtained in a similar way, thanks to the fact that the threshold forces for different fibers are uncorrelated.The variance of the total force is simply the sum of the variances of the f i , and thus using one obtains for the total force variance 2 (F ) Considering only integer value In the limit of a large , the above expression reduces to Therefore the relative standard deviation of the total force will approach (F )/͗F͘ → for >> 1. 1/ 3

͙
The entire distribution of the total force can be obtained analytically in the context of this simplified model.In the limit of a large length , one can simply observe that F is given by a sum of statistically independent random variables and the law of large numbers applies.Consequently, F will have a Gaussian distribution.Thus the above computed average and variance are sufficient to specify entirely the distribution.Keeping only the leading terms, one can write

Correlations
The above section gave some information concerning the total force distribution.To characterize the evolution of the signal F(U), one needs to work out the correlations in this function.More precisely, knowing the force F(U ), the question is, What can be inferred for F(U ϩ ⌬U )?
One can get some insight into this question by computing the expectation value of dF(⌬U ) = F(U ϩ ⌬U ) Ϫ F(U ).The force difference can again be written as a sum of independent statistical variables where j extends over the unbroken fibers after the crack tip U.
The random variables df j assumes the following values: where the last condition has been written so that the sum can extend over all fibers for j > U ϩ ⌬U.One has to consider also the fibers in the range U < j < U ϩ ⌬U, which are broken with probability 1, but were surviving under the displacement U with probability 1 Ϫ y j (U ).At the steady-state value of the force, ͗dF͘ = 0 by definition.The expectation value of ͗dF 2 ͘ can thus be obtained by summing up the variances of the df j Thus the squared force difference increases first linearly with ⌬U and saturates to a constant value equal to twice the variance of the force.The interpretation of this property is straightforward: the total force is the sum of the order of uncorrelated random variables.Thus over this length scale, F(u) behaves as a random walk.However, for larger distances, the fluctuations of F become uncorrelated.As a consequence of this observation, note that the typical variation of the force over a short displacement ⌬U scales as Consequently, the fluctuating part of the signal be-⌬U.

͙
comes nondifferentiable when the microstructural size goes to 0, keeping fixed (note that this construction implies a redefinition of the physical scale because one chose here to measure distances in terms of the microscopic distance, the fiber separation distance).It is important to stress that such a conclusion has already been shown for the fiber bundle model in the case where a global load sharing rule was used and thus for a homogeneous macroscopic loading.
If the limit of a bundle with an infinite number of fibers is expected to represent a continuum response, the fact that the constitutive response is not differentiable is a striking departure from traditional assumptions.In fact, it means that when this continuum limit is considered, a smoothing of the constitutive response is performed at the same time so that the response becomes differentiable.It follows that upon taking this limit, the information contained in the fluctuations of the response for large-size systems is lost.As will be seen in the next section, this information enlightens the crack propagation regime and yields a parallel between the fiber bundle model and cohesive crack models in the context of this study.

Avalanches
There exists various definitions of avalanches.The first one is the most natural, consisting of fixing a level of force and computing the distance ⌬U over which the crack can propagate.The avalanches are characterized by their statistical distribution, p 1 (⌬U, F ).
To give a global characterization of the signal without considering a specific value of the force, one may consider the above avalanches for any crack length U such that a fiber would break and then average over all breaking events.These avalanches are the ''forward'' avalanches.A similar construction can be performed for avalanches ending at each fiber breaking event.They correspond to a similar construction as the previous ones after reversing the arrow of time.They are called ''backward'' avalanches.In some cases such as invasive percolation, these two statistics have very different properties.Fig. 4 shows the computed forward avalanche distribution for the present model where a fixed displacement profile is set.Three values of have been used.The first one is = 100 (⅜) with 10 8 broken fibers, the second one is = 1,000 (ϫ) with 2.10 7 broken fibers, and the third one is = 10,000 (᭟) with 5.10 6 broke fibers.
The main observation is that the distribution exhibits two distinct behaviors according to the value of the avalanche size ⌬ with respect to a crossover value ⌬*.The first regime, ⌬ < ⌬*, is a power law p 1 (⌬) ϰ with an exponent The second regime, ⌬ > ⌬*, is also a power law but with a different exponent 2 = 2.05 Ϯ 0.10 ( 19) 2 Fits to both of these power laws are plotted in Fig. 4. Finally, the crossover scale ⌬* scales as : the graph used the scaled variables ⌬/ and scaled distribution to show that the 1 p (⌬) 1 three curves collapse onto a single master curve.This data collapse shows that, indeed, the crossover scale ⌬* is proportional to .
One can easily understand the value of the two exponents.The first one corresponds to a regime where the force versus crack length U displays correlations similar to a random walk.The forward avalanche, in this case, can be interpreted as the time required for a random walk to return to the origin.This well-known statistical problem is indeed a power law of exponent 3/2 in agreement with the first regime.Note that this behavior is exactly the one that has been established for the global load sharing fiber bundle with rigid boundaries (Hemmer and Hansen 1992;Hansen and Hemmer 1994), using also a mapping onto a random walk problem.
For large avalanches, ⌬ > ⌬*, the forces are uncorrelated.One thus may resort to this simple case to work out the avalanche statistics: let (t) be a random uncorrelated noise, with a distribution p() and cumulative distribution P() =  ͐ p(x) dx.Starting at a given value of = 0 , the probability that an avalanche is >⌬ is Q(⌬, 0 ) = P( 0 ) ⌬ because the different values are uncorrelated.The cumulative distribution P 1 of forward avalanches is obtained from the integration of the above Q distribution over all starting points of distribution p(), hence where one has used ⌬ >> 1.The avalanche distribution p 1 is obtained from the derivative of the cumulative distribution and leads to the power law p 1 (⌬) = ⌬ Ϫ2 for all distributions p().
In this problem, the large avalanches correspond to this case and, indeed, one observes an exponent 2 = 2.One sees in this particular example that a simple statistical analysis performed on the force signal allows one to extract the correlation length ϰ ⌬* without knowing it beforehand.This correlation length defines the size of the fracture process zone.
There is another definition of avalanches that might be easier to handle in experimental studies.In the above construction, one needs to know the entire signal F(U ).This requires an infinitely stiff mechanical testing device.In contrast, if a finite stiffness K is used, the test would not be able to follow the characteristic of the system unconditionally.After a local maximum (U 0 , F 0 ), U would jump directly to a new position . Such a construction defines biased avalanches characterized by the value of the stiffness K. K , K 2 .The avalanches are power-law distributed up to ⌬ c with an exponent 3 Ϸ 1.5.One has seen above that the fluctuating part of the force signal had a self-affine character with a typical variation over a scale ⌬U proportional to dF = A⌬U 1/2 .Thus, at small scales, the fluctuations dominate whereas at a large avalanche size, the bias takes over the fluctuations.Thus there is a characteristic scale for the largest avalanche size, such that K⌬ c = A(⌬ c ) 1/2 ; hence For smaller sizes, the bias can be neglected; hence, the biased statistics are again similar to a first return time distribution for a random walk.Hence a power law of exponent 3/2 is ex-pected consistently with the numerical results.This regime is, however, valid only for the case where K Ϫ2 << .Otherwise, above , avalanches cannot be distributed in a similar way because the force signal becomes uncorrelated.Solving for the case of an uncorrelated signal with a Gaussian distribution, one finds that the avalanches are uniformly distributed up to a maximum scale depending on K as ⌬ c ϰ K Ϫ1/2 .Fig. 7 illustrates this scaling behavior for a small value of .
Finally, the distribution of backward avalanches has been checked.One finds the same distribution as for forward avalanches, with the same two power laws.Superposition of forward and backward avalanche distributions is shown in Fig. 8.This result may appear as natural; however, in other models (Paczuski et al. 1995) (e.g., invasive percolation), a very different behavior is obtained for the two definitions.In this case though, the random walk analogy as well as the case of uncorrelated noise signal both support such a symmetrical behavior.

MODEL WITH DEFORMABLE BEAM Mean Deflection of Beam
One now comes back to the complete model where the beam is not assumed to have a constant displacement profile.Rather one wishes to compute here the average shape of the deformed beam.Introducing E and I, the young modulus and transverse geometrical inertia of the beam, respectively, one can write, in the spirit of a continuum modeling (i.e., for length scales much larger than the fiber separation), an equation for the mean deflection of the beam y(x) This equation holds for a uniform distribution of critical fiber extension between 0 and 1, and for y < 1, whereas d 4 y/dx 4 = 0 for larger y.The boundary conditions are y(ϱ) = dy(ϱ)/dx = 0, y(0) = 1 (imposed displacement), and d 2 y(0)/dx 2 = 0 (no torque being applied at the loading point).The reason one has to distinguish between y and Y is that the deflection is not a monotonous function of x.Because damage is irreversible, one has to compute the maximum damage having been met by the corresponding section of fibers.One does not know the analytical solution to this problem; however, one can see that the quadratic nonlinear term becomes important at a large distance from the origin.Thus the asymptotic shape will have the following expression: ͙ ͩ ͪ The oscillatory component is the one that makes the deflection nonmonotonous and thus requires the distinction between y and Y. Omitting this oscillatory component and the nonlinearity, one retrieves the exponential decay that was imposed in the simplified version of the model.A numerical integration based on a relaxation method has been done, and the results are shown in Fig. 9.

Avalanches
The problem requires a much longer computation time compared to the simplified model because of the larger number of degrees of freedom.For numerical convenience, only the process zone is dealt with, thus neglecting the presence of a few broken fibers ahead of this region.However, this zone is taken into account in the computation by introducing a boundary condition at the end of the process zone that represents an infinite beam connected to the substratum through intact fibers.This involves two relations between the derivatives of order 0-3 of the deflection function y(x).The length of the domain considered numerically is set equal to .
Fig. 10 shows forward avalanche distributions for two stiffnesses ratios and hence two values of .In both cases the previous results are recovered (i.e., two power laws with exponents 1 = 1.50 Ϯ 0.05 and 2 = 2.1 Ϯ 0.15).One checks the reliability of the measure by plotting and evaluating the mean shape of the beam during crack propagation.The results are very close to the simplified model.Hence, the fluctuation analysis from the force-displacement response, an accessible experimental information, provides information on the existence of a well-defined length scale related to the fracture process zone in the interface.

CONCLUSIONS
The microstructure of heterogeneous materials plays a crucial role in the rupture process.This paper introduces a model of a 1D crack propagation that includes randomness.Fluctuations that are encountered on the force-displacement response are studied through an avalanche distribution analysis.Two cases are considered.First, a simplified model with a prescribed displacement opening profile is considered along the interface and it could be solved analytically.Second, a more complete model, where the coupling is done between a flexible beam and the fiber bundle, is considered.The latter was studied only through numerical simulations.Analytical predictions are supported by the numerical result of the simplified model.The main conclusions are • Although the response shows large fluctuations, an analysis of the avalanche distributions provides information on a redistribution length that scales the fracture process zone.• In the limit of an interface of an infinite size, the analytical model yields a response that is nondifferentiable.It cannot be compared with the smooth response expected in the context of a continuous interface.Therefore, the fluctuations contain important information on the interface fracture process, which is lost in a continuum approach.• Two regimes are found on the avalanche distributions.
The first one is a power law with an exponent 1.5, which can be mapped onto a random walk.It concerns avalanches whose sizes are smaller than the fracture process zone.The second one is also a power law with an exponent 2, which corresponds to an uncorrelated noise.It governs large avalanche sizes.• A ''biased avalanche'' has been introduced, which takes into account a finite stiffness of the loading device or of the system itself.Similar to the previous definition of avalanches, it exhibits two different regimes depending on the stiffness and on the size of avalanches.• The crossover avalanche size between the two regimes gives natural access to the extent of the fracture process zone.This information may prove to be a precious experimental tool to provide access to this a priori unknown scale.

FIG. 2 .
FIG. 2. Schematic Representation of ''Zip'' Model (Point Where Displacement Is Imposed May Move Along x-Axis) FIG. 3. Evolution of Mean Total Force (Bold Curve) Scaled by and Shifted by Ϫ1/2 as Function of Displacement U ( Has Been Set to 10, to 1; Long Dashed Lines Are Boundary of ͗F ͘-Values That Are Obtained from Eqs. (6) and (8); Thin Curve Represents Evolution of Mean Force for = 20 and Shows How Fast Variation with Decreases)

FIG. 7 .
FIG. 7. Statistical Distribution FIG. 5. Force Signal (Thin Line), When Followed with Loading Device Having Finite Stiffness K, Gives Access to Bold Curve Envelope, Where Decreasing Parts Have Slope ϪK Fig. 5 illustrates these avalanches.Fig. 6 shows the statistical distribution p b (⌬, K) obtained for different values of the stiffness K.The distribution has been scaled by whereas the avalanche size has been scaled by 2 1