S ep 2 01 9 Creep failure in a threshold activated dynamics : Role of temperature during a subcritical loading

Creep is time-dependent deformation of solids at relatively low stresses, leading to the breakdown with time. Here we propose a simple model for creep failure of disordered solids, in which temperature and stress are controllable. Despite its simplicity, this model can reproduce most experimental observations. Time dependence of the strain rate is well fitted with power laws resembling the Omori-Utsu and the inverse Omori laws in the primary and the tertiary creep regimes, respectively. Distribution of the creep lifetime obeys the log-normal distribution, and the average creep lifetime decays in a scale-free manner with the increasing stress. The above results are in good agreement with experiments. Additionally, the mean avalanche size as a function of temperature exhibits a series of jumps, and finite size scaling implies the existence of phase transitions.


I. INTRODUCTION
Failure point and the failure processes of disordered solids and composite materials depend on temperature, pressure, and driving conditions such as applied stress or strain rate [1]. Importantly, materials can deform and break with time even if the applied stress is below the critical value. This phenomenon is called as creep failure or creep rupture. Study of creep failure is essential in many contexts such as constructions, in which the building blocks are subject to the load for a considerable duration. The time elapsed until creep failure is known as the creep lifetime. From the practical point of view, estimate of creep lifetime leads to a direct forecast of the catastrophic failure event and therefore important in materials science. The creep lifetime depends on many ingredients, among which the applied stress and the temperature have prominent effects [2].
Phenomenology of creep includes several power laws. Suppose that the stress is applied to a specimen at t = 0 and kept constant thereafter. Then the specimen starts to deform, but the strain rate decays with time as t −p [3]. After a certain lapse of time, deformation accelerates toward breakdown with the strain rate increasing as (τ c −t) −p ′ , where τ c denotes the creep lifetime of the specimen. The two exponents, p and p ′ , are not generally the same, but typically range from 0.6 to 1.0 [3][4][5][6][7]. The latter power law describing acceleration toward breakdown has been of interest from practical point of view, since the time of catastrophic breakdown could be predicted by fitting the real-time data to the above formula [8][9][10]. Interestingly, similar power laws are also known for earth- * Electronic address: subhadeep.roy@ntnu.no † Electronic address: hatano@ess.sci.osaka-u.ac.jp quakes [11][12][13], which involve much larger scales than the experimental specimens. In this case, the creep lifetime may correspond to the time of earthquake. Therefore, to understand the physics behind these power laws may lead us to some novel principles that are common to fracture phenomena across different scales, as well as to develop useful tools for forecasting catastrophic failure events.
Since the physical mechanisms behind creep may be the accumulation of damage and plastic strain in a system, models for creep should take such aging processes into account. Structural randomness and heterogeneity in solids should be also considered appropriately. However, the fracture of heterogeneous solids is yet hard to handle with the elasticity theory for continuum, and therefore there have been some alternative approaches. A pioneering work by Main [14] is based on the subcritical crack growth dynamics and the phenomenological Voight's model for precursory strain [15]. Assuming a heterogeneous feedback system, the model can reproduce both the power laws with the exponents directly related to those assumed at the microcrack level. Another class of approaches makes use of a simple and intuitive model for disordered solids [5,[16][17][18][19][20][21][22][23][24], which is known as the fiber bundle model [25][26][27][28][29][30][31]. Since temperature plays an important role in creep rupture, thermal fluctuations are explicitly incorporated in these models in the form of noise or probabilistic time evolution. They appear to be successful in explaining the power law behaviors of the time-dependent strain rate as well as the temperature and the stress (or strain) dependences of the creep lifetime [16][17][18][19]. Instead of introducing stochasticity, Hidalgo et al. introduce the time evolution equation for the strain of each constituent (i.e., fiber) based on the Kelvin-Voigt rheology [21,22], and obtain the temperature and the stress dependences of creep lifetime that are somewhat different from those in the above-mentioned stochastic models [16][17][18][19]. Danku and Kun consider a damage accumulation process in each fiber by introducing the evolution equation for the damage variable [24]. Their model reproduces the inverse Omori law with exponent p ′ = 1, although it increases up to 5 if the smaller events are disregarded.
Contrastingly, much simpler models can exhibit creeplike behaviors even in the absence of any thermal fluctuations, damage variables, or rheological constitutive laws [32][33][34]. In particular, the present authors derived both the power laws for creep only by introducing the time evolution in a simple fiber bundle model [34]. However, the obtained exponents (p = 1.8 and p ′ = 2) are too large to be comparable to any experimental values.
Here we propose a slightly modified model, in which the effect of temperature is included. Probabilistic algorithm is adopted for the failure criterion of fibers, bringing thermal fluctuations in the model in addition to the time-independent (i.e., quenched) randomness. The interplay between these two kinds of randomness yields nontrivial time-dependent behaviors. We have numerically studied the model focusing on the effects of temperature and applied stress on the creep behaviors: statistics of the creep lifetime including the distribution function, the time evolution of strain rate, and the avalanche statistics. In the next section, we have provided a detailed description of the conventional fiber bundle model along with the modifications we have made for the present work. This is followed by the numerical results together with some brief comments and discussions on the future scope as a continuation of the present observations.

II. DESCRIPTION OF THE MODEL
After its introduction by Pierce in 1926 [25], the fiber bundle model has been proven to be a useful and yet arguably the simplest model to understand the failure process in disordered solids [26,[29][30][31]. The fiber bundle model consists of L vertical fibers in between two parallel bars. A load F is applied to the bars to create a stress σ = F/L per fiber. Each fiber has an individual strength chosen from a threshold distribution randomly. This heterogeneity in the strength may be regarded as a quenched randomness in the model, and the dispersion of threshold distribution measures the strength of disorder in the model. When the applied stress exceeds a threshold value, the corresponding fiber breaks irreversibly, and the load borne by that fiber is redistributed within the model; either among all the surviving fibers (mean field model) or among the surviving nearest neighbors only (local stress concentration). Due to such redistribution, the local stress values of some fibers increase and that can lead to further breaking and redistribution. This is an avalanche. After a certain number of avalanches, the model breaks completely or relaxes to a stable state with nonzero surviving fibers. In the latter case, the external load F needs to increase to break the next weakest fiber leading to further avalanches. This process continues un-til all the fibers are broken. The applied load just before the global failure is referred to as the critical load F c .
Fiber bundle model has some variations to include thermal effects in fracture [16-18, 23, 35-38]. Here we also adopt a probabilistic rule that is dominated by temperature. Let us assume that any fiber breaks with a certain probability that depends on the temperature T . In this study, the probability of the ith fiber to break at time t is given by Here σ(t, i) and σ th (t, i) are, respectively, the local stress and the stress threshold of ith fiber at time t, and P ′ is a constant chosen to be unity hereafter. Note that any fibers can rupture even if the local stress is less than its threshold value. If σ(t, i) ≥ σ th (t, i), P r reaches unity and the fibers break with probability 1. Note also that this probabilistic model reduces to a conventional model in the limit of T → 0, in which a fiber can break only when σ(i) ≥ σ th (i).
The time evolution of the system is as follows. At t = 0, the load is applied to the system, leading to a nonzero fracture probability of P r (0, i). This P r (0, i) is then compared with a random number P * (0, i) generated uniformly within the interval of [0, 1]. If P r (0, i) > P * (0, i), the ith fiber breaks and the borne load is redistributed. At the next time step (t = 1), a different rupture probability is calculated based on the new local stress profile σ(1, i). The model keeps evolving at each time by comparing the two probabilities, P r (t, i) and P * (t, i). The present authors adopt essentially the same time evolution rule in a fiber bundle model with a deterministic rupture rule (i.e., at T = 0) and find creeplike time evolution [34]. Here we extend this model to investigate the effect of stochasticity, the amplitude of which may be proportional to temperature.
Equation (1) suggests that P r increases with temperature T and the local stress σ(t, i). Then the nature of breakdown of the system will be determined by the interplay between the local stress profile, the individual strength of the fibers, and the temperature. In the present study, we concentrate on the mean field model, in which σ(t, i) are uniform (independent of i) as a result of democratic load redistribution.
The stress threshold of each fiber, σ th (i), is a random quantity to be sampled from a certain probability distribution, ρ(σ th ). Unless otherwise indicated, we adopt a uniform distribution with the mean of 1/2 and the halfwidth of δ:

III. NUMERICAL RESULTS
Numerical results are produced for the mean field model with system sizes ranging between 10 3 and 10 5 . To compute average values of physical quantities such as creep lifetime, 10 4 to 10 5 configurations are sampled at each system size. Here we have mainly studied the dispersion in creep lifetime and the evolution strain rate with time. Additionally, a series of dynamical transitions observed during the creep failure will be discussed in detail. In this section, we discuss the distributions of creep lifetime in the parameter space of σ − T to understand whether any extreme statistics is associated with it. Hereafter the creep lifetime is denoted as τ c , which varies from sample to sample. To compute the distributions of creep lifetime, 10 5 samples are computed for the system size of 10 3 . Figure 1 shows that the peak of the distribution shifts to lower values as either of the parameters, temperature or applied stress, is increased. This happens as the probability of rupture increases with increasing σ or T . A simultaneous increment in both parameters shifts the peak even faster.
The distribution of creep lifetime develops a tail, suggesting the existence of extreme statistics. As is shown in Fig.2(a) and (b), a closer look on the distribution functions reveals that they are well fitted with the log-normal distribution given below.
where v is the variance and µ is the mean of log τ c . The values of (µ, v) are shown in the legends of individual panels of Fig. 2. They depend on the applied stress and the temperature as shown in Fig.2 (b), (c), (e), and (f). We observe that the mean value decreases gradually as either the temperature or the applied stress increases. On the other hand, the variance is insensitive to the applied stress but decreases with increasing temperature. The lifetime distributions are also compared with the three-parameter Weibull distribution [39]. The comparison is shown in Fig.2(g) for T = 0.08 and σ = 0.14. Apparently, the log-normal distribution fits the simulation result better, particularly at the tails. This behavior remains at any temperature and stress as well as in the limit T → 0, at which the rupture events are deterministic [33].
At a constant temperature and varying applied stress, we observe the following scaling: with β = 4.0. The scaling is demonstrated in Fig.3(a) for different applied stress but at a constant temperature T = 0.08. The inset of the same figure shows the unscaled behavior. The above scaling also tells that, at constant temperature T , the average lifetime decreases with the applied stress in a scale-free manner with a temperature dependent exponent, β.
The behavior is examined numerically and shown in Fig.  3(b). The exponent β is a decreasing function of temperature.

B. Time evolution of strain rate
Next, we have studied the time evolution of the strain rate,ǫ. Since the load F is constant and the stiffness of the fibers are uniform, the strain (ǫ) in the bundle is proportional to F/L t , where L t denotes the number of intact fibers at time t. Then the strain rate is given bẏ where the stiffness of the fibers are assumed to be unity. Fig.4 shows the time evolution of strain rate for several sets of parameters: disorder δ, applied stress σ and temperature T . Fig.4(a) shows the strain rate as a function of time at several values of temperature, where the stress and the disorder strength are common (σ = 0.05 and δ = 0.5). We observe the following power-law behavior for the strain rateǫ.ǫ where t is the time elapsed after the loading and c is the time constant for the onset of scale-free behavior. This is essentially identical to Andrade creep, t −p , if the time constant c is negligible. However, note that the time constant c is nonzero in the present model as shown in Fig.4(c). The time constant c increases as the extent of disorder δ decreases as shown in the inset of Fig.4(c), but insensitive to the applied stress. The exponent p depends on temperature, and decreases gradually as the temperature increases as shown in Fig.4(a). The inset shows the variation of p with respect to the temperature. In a wide temperature range (0.02 ≤ T ≤ 0.1), the exponent p changes gradually from 1 to 0.67. This range of exponent values is in good agreement with experiments. The exponent p is insensitive to the applied stress σ as shown in Fig.4(b), and depends only gradually on the dispersion in the threshold distribution as shown in the inset of Fig.4(c).
In the tertiary creep, in which deformation accelerates toward breakdown, the strain rate can be fitted with a power law again as shown in Fig. 4: which is known as the inverse Omori law. Importantly, the characteristic time c is nonzero again. As shown in Fig. 5 (a), the c-value increases as the extent of disorder δ declines. Figure 5 (b) and (c) show that the cvalue also increases as the temperature increases (b), or the stress decreases (c). Namely, the acceleration lasts longer towards breakdown at lower temperature, lower stress, and/or stronger disorder. In terms of predictability of the creep lifetime, the number of parameters should be as small as possible. In this respect, the failure may be more predictable for vanishing c (i.e., large disorder). This tendency of predictability is consistent with an ex-periment with controllable disorder [10]. On the other hand, the exponent p′ appears to remain for the parameter range investigated here, although the power law itself is less clear if c is large. The above two power-law behaviors, the Omori and the inverse Omori laws, are common to a simple deterministic model at T = 0 [34]. However, the parameters behave differently in most respects. For T = 0 case, the exponent p ′ is approximately 2 irrespective of the randomness and the stress, whereas it is a decreasing function of disorder and temperature in the present case. For the primary creep, the characteristic time c in Eq. (7) is insensitive to the stress in the present model, whereas it increases with the stress at T = 0. For the tertiary creep, the characteristic time c in Eq. (8) vanishes at T = 0, whereas it is nonzero in the present model. Contrastingly, the temperature dependence of the characteristic time c in Eq. (8) extrapolates to the zero temperature case.

C. Transitions in abruptness
In this section, we have studied the abruptness of failure. Denoting the number of unbroken fibers as L t , we define the decreasing rate of fibers as s(t) = L t − L t+1 . Larger s thus implies that the breakdown is more abrupt. Since this s varies from sample to sample, the ensemble average, which is denoted by s , is taken over 10 4 con- figurations. Note that s /L = τ −1 c . Figure 6 shows the behavior of s /L. The applied stress σ is varied from 0.1 to 0.4, while the disorder strength δ ranges in between 0.1 and 0.3. At high temperatures, s ∼ L, suggesting that the total bundle breaks in a single time step. As T decreases, s /L decreases abruptly at a certain temperature. As the temperature decreases further, similar drops of s /L occur repeatedly but the amplitude of abrupt change declines gradually at lower temperatures. The jump heights gradually decreases as we go to lower T values. The first jump is from 1.0 to 0.5, and the second jump is from 0.5 to 1/3, and so on. The height of nth jump, denoted by ∆h n , is ∆h n /L ≃ 1 n(n + 1) .
This is understandable since s /L = τ c and τ c goes from n to n + 1 for the nth jump. Figure 6 also shows that the jumps gradually disappear as the stress σ increases. For example, there are a number of jumps with decreasing jump heights for σ = 0.1. On the other hand, for σ = 0.4, we observe only one jump from s /L = 1 to 0.5. This is because the failure occurs very abruptly at higher stresses, taking only a few time steps. Thus, there are only a few large avalanches at higher stress, and subsequent transitions to smaller avalanches are absent. This behavior remains almost unaltered with increasing disorder strength, whereas the transition temperatures are affected.
To study these jumps more closely, we have carried out the finite size scaling. Figure 7 shows the behavior of s /L at several system sizes ranging from L = 5 × 10 2 to 5 × 10 4 . The system size scaling is given as follows: where T c n is the transition temperature for the nth jump, and the exponent α is a constant. The inset of Fig.7 shows the system size scaling for n = 1 with T c 1 = 2.12. The exponent α is estimated as 0.45 irrespective of the stress and the temperature. We get a series of T c n values from the same scaling for different jumps. This scaling relation suggests that jumps in s /L with varying temperature are indeed phase transitions. Fig.8(a) and Fig.8(b) shows the variation of T c n with n for different values of σ and δ. Both behaviors suggest a decrease in T c n with n. In addition, the n values itself decreases as we go to higher σ values, suggesting a lesser number of jumps and at the same time big avalanches. On the other hand, n values remain invariant of the disorder strength. For a particular jump (n-value), the critical temperature decreases with increasing applied stress and decreasing strength of disorder.
Finally, Fig.9 shows the variation of T c n for the first and second transitions, when both applied stress σ and strength of disorder δ varies simultaneously. The figure shows two planes corresponding to the first (n = 1) and second (n = 2) transitions producing three separate regions. The region above the plane n = 1 corresponds to instantaneous failure where the total bundle breaks in a single avalanche producing a unit creep lifetime (τ c = 1). In this region, the distribution of τ c is a delta function at 1.

IV. DISCUSSIONS
In the present work, the effect of temperature in creep failure is probed by a fiber bundle model with a proba-bilistic algorithm. As a result of the stochasticity, the model can fail even if the applied stress is below the critical value at T = 0. The dynamical properties of the model, including the exponents for power-law behaviors, are comparable to those of real materials. The model exhibits a large scatter in the creep lifetime, which can be described by a log-normal distribution with the temperature-dependent mean and variance. The mean lifetime depends on the applied stress in a scale-free manner. A detailed study of the time-dependent strain rate enables us to detect the characteristic power laws with a nonzero time constant. The exponents depends slightly on the temperature, the range is comparable with the phenomenology known as the Andrade creep.
For a practical purpose, the average lifetime alone does not provide us with sufficient information on creep failure since the lifetime distribution is skewed [41]. Both the log-normal and the Weibull distributions are observed for the lifetime of some alloys [42][43][44]. The Weibull distribution was also found for STS304 stainless steel [49]. In our study, the creep lifetime distribution is fitted more convincingly with the log normal distribution than the Weibull distribution. Another major difference is in our case the Weibull distribution contains three parameters while the above experimental creep lifetimes are fitted well with two-parameter Weibull distribution. In addition, the stress dependence of the lifetime in our model agrees well with experimental behaviors [40,44,45]. Overall, the present numerical results are mostly consistent with experimental results on materials [40][41][42][43][44][45][46][47][48][49].
Understanding of creep behavior is a fundamental problem in material science and engineering. Prior knowledge on the creep lifetime for a particular material reduces the threat of sudden catastrophic failure and increases predictability in the failure process. At the same time, sufficient knowledge of creep dynamics makes us aware if the system is approaching the global failure. In this paper, we have presented a detailed study of the creep dynamics in a disordered system. Our next step is to explore the model with other modifications such as the inclusion of stress concentration around a broken fiber. Application of time-dependent loading allows us to explore the possibility of fatigue within the model.