Variation of elastic energy shows reliable signal of upcoming catastrophic failure

We consider the Equal-Load-Sharing Fiber Bundle Model as a model for composite materials under stress and derive elastic energy and damage energy as a function of strain. With gradual increase of stress (or strain) the bundle approaches a catastrophic failure point where the elastic energy is always larger than the damage energy. We observe that elastic energy has a maximum that appears after the catastrophic failure point is passed, i.e., in the unstable phase of the system. However, the slope of elastic energy vs. strain curve has a maximum which always appears before the catastrophic failure point and therefore this can be used as a reliable signal of upcoming catastrophic failure. We study this behavior analytically for power-law type and Weibull type distributions of fiber thresholds and compare the results with numerical simulations on a single bundle with large number of fibers.


I. INTRODUCTION
Accurate prediction of upcoming catastrophic failure events has important and far-reaching consequences. It is a central problem in material science in connection with the durability of composite materials under external stress [1][2][3][4][5]. The same problem exists at a large scale (field-scale) associated with mine and cave collapses, landslides, snow avalanches and the onset of earthquakes due to plate movements [6,7]. In medical science, understanding fracturing of human bones exposed to a sudden stress is an important research area [8]. These phenomena belong to the class of phenomena called stressinduced fracturing, where initially micro-fractures are produced here and there in the system and at some point, due to gradual stress increase, a major fracture develops through coalescence of micro-fractures and the whole system collapses (catastrophic event). Such stress-induced failures occur also in very different domains -for example, in breakdown of social relationships and mental health [9,10].
The central question is -when does the catastrophic failure occur? Is there any prior signature that can tell us whether catastrophic failure is imminent? The inherent heterogeneities of the systems and the stress redistribution mechanisms (inhomogeneous in most cases) make things complicated and a concrete theory of the prediction schemes, even in model systems, is still lacking.
In this article, we address this problem (prediction of catastrophic events) in the Fiber Bundle Model (FBM) which has been used as a standard model [11][12][13][14] for fracturing in composite materials under external stress. We will show theoretically that in the Equal-Load-Sharing (ELS) model: 1) At the catastrophic failure point, the elastic energy is always larger than the damage energy.
2) The elastic energy variation shows a distinct peak before the catastrophic failure point and this is a universal feature, i.e., it does not depend on the threshold distribution of the elements in the system. 3) The energy release during final catastrophic event is much bigger than the elastic energy stored in the system at the failure point. Our numerical results show perfect agreement with the theoretical estimates.
We organize our article as follows: After the brief introduction (section I), we define the elastic energy and the damage energy in the Fiber Bundle Model in section II. In sections III and IV we calculate the elastic and damage energies of the model in terms of strain or extension. In several subsections of sections III and IV we explore the theoretical calculations for power-law type and Weibull type distribution of fiber thresholds. Simulation results are presented and numerical results are compared with the theoretical estimates in these sections. We present a general analysis of elastic energy variations and existence of an elastic-energy maximum in section V. In section VI we identify the warning sign of catatrophic failure by locating the inflection point. Here, in addition to uniform and Weibull distributions, we choose a mixed threshold distribution and present the numerical results, based on Monte Carlo simulation, to confirm the universality of the behavior in the ELS models. Finally, we keep some discussions in section VII.

II. THE FIBER BUNDLE MODEL
The fiber bundle model consists of N parallel fibers placed between two solid clamps (figure 1). Each fiber responds linearly with a force f to a stretch or extension ∆, where κ is the spring constant. κ is the same for all fibers. Each fiber has a threshold x assigned to it. If the stretch ∆ exceeds this threshold, the fiber fails irreversibly. When the clamps are stiff, load will be redistributed equally on the surviving fibers and this is called the equal-load-sharing (ELS) scheme. Throughout this article we work with ELS models only. The fiber thresholds are drawn from a probability density p(x). The corresponding cumulative probability is When the fiber bundle is loaded, the fibers fail according to their thresholds, the weaker before the stronger. Suppose that n fibers have failed. At a stretch ∆, the fiber bundle carries a force where we have defined the damage When N is large enough, d may be treated as a continuous parameter.
We will now assume that the stretch ∆ is our control parameter. We can construct the energy budget according to continuous damage mechanics [1,15]. Clearly, when we stretch the bundle with external force, work is done on the system. At a stretch ∆ and damage d, the elastic energy stored by the surviving fibers is The damage energy of the failed fibers is given by The total energy at stretch ∆ and damage d is then

III. ELASTIC ENERGY AND DAMAGE ENERGY AT THE FAILURE POINT
We are going to analyze the energy relations when the bundle is in equilibrium. We know that there is a certain value, ∆ = ∆ c , beyond which catastrophic failure occurs and the system collapses completely. We are particularly interested in what happens at the failure point. Is there a universal relation between elastic energy and damage energy at the failure point?
When N is large, we can reframe Eqns. (5,6) and express the energies in terms of external stretch (or extension) ∆ as and The force on the bundle at a stretch ∆ can be written as The force must have a maximum at the failure point ∆ c , therefore setting dF (∆)/d∆ = 0 we get A. Uniform threshold distribution We start with the simplest threshold distribution: the uniform distribution, which is well-known in fiber bundle research [14]. For a uniform fiber threshold distribution within the range (0, 1), p(x) = 1 and P (x) = x. Therefore we get, from Eqn. (11), Now putting ∆ c = 1/2 in Eqns. (8,9), we get and Therefore, the ratio between damage energy and elastic energy at the failure point (∆ c ) is B. Power-law type threshold distribution Now we move to a general power law type fiber threshold distributions within the range (0, 1), The cumulative distribution takes the form We insert the expressions for p(x) and P (x) into Eqn. (11) and find the critical extention We can calculate the elastic energy and damage energy at the failure point ∆ c : and Plugging in the value of ∆ c (Eqn. 18) into the above equations for elastic energy and damage energy we end up with the following relation: Clearly, the ratio depends on the power factor α ( Figures  2 and 3). When α = 0, the threshold distribution reduces to a uniform distribution and we immediately go back to Eqn. 15.

C. Energy balance
It is easy to show that the work done on the system up to the failure point ∆ c is equal to the sum of the energies E e and E d . The total work done on the system can be calculated as Inserting the expression for F (∆) into the integral we get, for a general power law type distribution, which is the total of elastic energy and damage energy, Eqns. 19,20). In fact, this is the energy conservation relation in thermodynamic sense.

D. Energy release during the final catastrophic avalanche
It is known that when the extension exceeds the critical value ∆ c , the whole bundle collapses via a single avalanche called the final or catastrophic avalanche [14]. Can we calculate how much energy will be released in this final avalanche? It must be equal to the total damage energy of the fibers between threshold values ∆ c and the upper cutoff level of the fiber thresholds for the distribution in question.
We calculate the damage energy of the final avalanche for power-law type distributions as It is important to find out whether the damage energy for the catastrophic avalanche has a universal relation with the elastic or damage energies at the failure point. we can write the following relation: As E e ∆c = (3+α)E d ∆c , we can easily get the other relation: And That means the energy release during final catastrophic avalanche is much bigger than the the elastic energy stored in the system just before failure (final stable state when ∆ = ∆ c ).
It is commonly believed that during catastrophic events like earthquakes, landslides, dam collapses etc., the accumulated elastic energy releases through avalanches [6,7]. We observe a different scenario in this simple fiber bundle model where the system is doing work during the catastrophic failure phase as the external force is still acting on the bundle. As a result, the energy release (during catastrophic failure event) becomes much bigger than the elastic energy stored at the final stable phase.

IV. ENERGY-ANALAYSIS FOR WEIBULL DISTRIBUTION OF THRESHOLDS
We now consider the Weibull distribution, which has been used widely in material science [14]. The cumulative Weibull distribution has a form: where k is the Weibull index. Therefore the probability density takes the form: As the force has a maximum at the failure point ∆ c , inserting P (x) and p(x) values in expression (Eqn. 11) we get From the above equation we can easily calculate the critical extension value as The elastic energy at the critical extension ∆ c is: and the damage energy is we get This integral is exactly calculable for k = 1 and k = 2.
A. Weibull distribution with k = 1 For Weibull index k = 1, ∆ c = 1 and the damage energy expression at the failure point takes the form Using integration by parts we arrive at the result We get the elastic energy at the failure point directly by putting k = 1 in Eqn. (33), Therefore, the ratio between damage and elastic energies at the failure point for Weibull distribution with k = 1 is In Figure 6, we have shown numerical results of the variation of elastic and damage energies with strain for Weibull distribution (with k = 1). The theoretical estimates of the ratio between damage and elastic energies at the failure point are compared with numerical results in Figure 8. For Weibull index k = 2, ∆ c = 1/ √ 2 and the damage energy expression at the failure point is Again, using integration by parts we arrive at the result  We get the elastic energy at the failure point directly by putting k = 2 in Eqn. (33): Therefore, the ratio between damage and elastic energies at the failure point for Weibull distribution with k = 2 is In Figure 7, we have shown numerical results of the variation of elastic and damage energies with strain for Weibull distribution (with k = 2). The theoretical estimates of the ratio between damage and elastic energies at the failure point is compared with numerical results in figure 8. In Appendix A, we give a general argument that elastic energy will be alaways bigger than damage energy at the critical (failure) point.

V. ELASTIC ENERGY MAXIMUM
There are two distinct phases of the system: A stable phase for 0 < ∆ ≤ ∆ c and an unstable phase for ∆ > ∆ c . If we plot the elastic energy and damage energy vs. ∆, we see ( Figure 9) that damage energy always increases with ∆ but elastic energy has a maximum at a particular value of ∆, let us call it ∆ m . Can we calculate the exact value of ∆ m for a given threshold distribution? Is it somehow connected to ∆ c ? In this section we are going to answer these questions.  We recall the elastic energy expression Eqn. (8). If we differentiate the elastic energy with respect to the extension ∆, we get Which is 0 at ∆ m , with If we consider a general power law type distribution p(x) = (1 + α)x α , within (0, 1), we can write Therefore we can conclude that ∆ m is bigger than ∆ c , i.e., elastic energy shows a maximum in the unstable phase (Figures 9 and 10). A more general treatment for the relation between ∆ m and ∆ c is given in the Appendix B.

VI. ELASTIC ENERGY INFLECTION POINT: THE WARNING SIGN OF CATASTROPHIC FAILURE
Are there any prior indications of the catastrophic failure (complete failure) of a bundle under stress? In the fiber bundle model, although the elastic energy has a maximum, it appears after the critical extenstion value, i.e., in the unstable phase of the system. Therefore it can not help us to predict the catastrophic failure point of the system. However, if we plot dE e /d∆, the change of elastic energy with the change of extension value ∆, we see that dE e /d∆ has a maximum and, most importantly, this maximum appears before the critical extension value ∆ c (Figures 11 and 12). In this section we calculate the particular value of ∆ at which dE e /d∆ has a maximum. Let us call this ∆ value ∆ max . We will also see whether there is a relation between ∆ max and ∆ c .

A. Theoretical analysis
We recall the expression for the derivative of elastic energy with respect to strain of extension Eqn. (45). Taking derivative of the equation, we get (49) Setting d 2 E e (∆)/d∆ 2 = 0 at ∆ = ∆ max we get for a general power law type distribution This expression confirms that ∆ max < ∆ c for α ≥ 0. For a Weibull distribution with index k, we can write  The solution (of d 2 E e (∆)/d∆ 2 = 0) with (−) sign is the acceptable solution for the maximum. Hence, A more general argument is given in Appendix C for the relation between ∆ max and ∆ c .

B. Comparison with simulation data
In Figures 11 and 12 we compare the simulation results with the theoretical estimates. The simulations are done for a single bundle with large number (N = 10 7 ) of fibers and the agreement is convincing. We have used Monte Carlo technique to generate uncorrelated fiber thresholds that follow a particular statistical distributions (uniform and Weibull distributions). It is obvious that in simulations we can measure energy values in the stable phase only.

C. Simulation results for a mixed threshold distribution
Now we choose a mixed fiber threshold distribution, i.e., a threshold distribution. Can we see similar signature (maximum of dE e /d∆ appears before ∆ c ) as we have seen in previous section? The chosen distribution is a mixture of uniform distribution and Weibull distribution (k = 1) which is shown in Figure 13. We assign strength thresholds to N/2 fibers from a uniform distribution and to N/2 fibers from a Weibull distribution. The simulation result (figure 14) reveals that dE e /d∆ has a maximum which appears before the failure point ∆ c and ∆ max is somewhere in between the respective ∆ max values for uniform and Weibull threshold distributions -as expected intuitively.

VII. DISCUSSIONS
The Fiber Bundle Model has been used as a standard model for studying stress-induced fracturing in composite materials. In the Equal-Load-Sharing version of the model, all intact fibers share the load equally. In this work we have chosen the ELS models and we have studied the energy budget of the model for the entire failure process, starting from intact bundle up to the catastrophic failure point where the bundle collapses completely. Following the standard definition of elastic and damage energies from continuous damage mechanics framework, we have calculated the energy relations at the failure points for different types of fiber threshold distributions (power law type and Weibull type). At the critical or catastrophic failure point, the elastic energy is always larger than the total damage energy. Another important observation is that the elastic energy variation has a dis- tinct peak before the catastrophic failure point. Also, the energy-release during final catastrophic event is much bigger than the elastic energy stored in the system at the failure point (see section III and section IV). Our simulation results on a single bundle with large numbers (10 7 ) of fibers, show perfect agreement with the theoretical estimates. We have chosen a single bundle, keeping in mind that for prediction purposes it is important and necessary that the warning sign can be seen in a single sample.
These observations can form the basis of a prediction scheme by finding the correlation between the position (strain or stretch level) of elastic energy variation peak and the actual failure point. Moreover, it is also possible to predict the size (energy release) of the final catastrophic event by measuring the stored elastic energy of the system at the failure point.
Our observations in this work have already opened up some scientific questions and challenges: what happenes for Local-Load-Sharing (LLS) models [16]? Does the elastic energy variation show similar peaks before the catastrophic failure point? Can we measure and analyze the elastic energy during a rock-fracturing test [17] in terms of the applied strain?
Our next article will resolve some of these issues -we are now working on energy budget of LLS models.

Acknowledgments
This work was partly supported by the Research Council of Norway through its Centers of Excellence funding scheme, project number 262644.

VIII. APPENDIX
As stated in section II, the elastic energy in the system at extension ∆ is where P (∆) is the cumulative probability distribution of the fiber thresholds. The force per fiber σ = F/N required to continue the breaking process at a given extension ∆ is The critical extension ∆ c where the bundle collapses is hence given by A. Elastic versus damage energy at the critical point Numerical data seems to suggest that E e > E d at the critical point for most threshold distributions. Let us try to prove this analytically by investigating the difference between elastic and damage energy: (57) The derivative of this energy difference is (58) We can now express the energy difference in terms of the forces acting on the fiber bundle. We integrate this expression to find by partial integration. In particular, this gives the result at the critical point. Since σ c = max σ(∆), we see that E dif f (∆ c ) > 0 for all threshold distributions. The only exception possible is a threshold distribution with a constant force σ(∆) = σ c . But this results in a lower cut-off ∆ 0 = ∆ c > 0 (for the threshold distribution to be normalizable), and then E e (∆ c ) > E d (∆ c ) = 0.

B. Elastic energy maximum point
First, rewrite Eqn. (56) as This definition of g(∆) will be useful in the following derivations. The maximum of the elastic energy is found at an extension ∆ m , which is given by i.e., Comparing this expression to Eqn. (61) allows us to find a relation between ∆ c and ∆ m . We investigate the function g(∆): It is clear from Eqn. (64) that for a threshold distribution with only a single maximum in the load curve, g(∆ m ) = 1/2 must occur in the unstable phase, i.e., ∆ m > ∆ c .

C. Elastic energy inflection point
The elastic energy maximum occurs after the critical point and is hence unsuitable as a predictor for failure. But what about the maximum of the derivative of the elastic energy, the inflection point ∆ max ? As stated in the section VI: Setting this second derivative to zero and rearranging terms gives the equation 2(1 − P (∆ max )) − 4∆ max p(∆ max ) = ∆ 2 max p ′ (∆ max ).
(66) To find a relation between ∆ max and ∆ c we substitute ∆ max p ′ (∆ max ) = −2p(∆ max ) − d 2 σ/d∆ 2 ∆max and then divide by 2∆ max p(∆ max ). After rearranging terms the result is From this expression and Eqn. (64) we can see that at ∆ max , dσ/d∆ and d 2 σ/d∆ 2 must have different signs.
Let's once again assume that we are working with a threshold distribution that has only a single maximum in its load curve. Then dσ/d∆ > 0 corresponds to the stable phase and dσ/d∆ < 0 corresponds to the unstable phase. We see that any threshold distribution with d 2 σ/d∆ 2 < 0 everywhere in the unstable phase must have ∆ max in the stable phase, i.e., ∆ max < ∆ c . This is a general (but weak) condition that is sufficient, but not necessary, for ∆ max < ∆ c .