Random fuse model in the presence of self-healing

We now suppose that the broken resistors can be restored by some self-healing process. We are addressing the self-healing process from a generic theoretical point of view, without considering any specific mechanism. The focus is on the effect of the self-healing process, i.e. the repair of a local microcrack, modelled as the reattachment of a broken link of the lattice, and on the time delay of this repair. These two features are common to all self-healing processes, with different parameter values, and can be included in the framework of the RFM, independently of the specific features of the practical method adopted to achieve self-healing in applications.

We define the self-healing rate η as the ratio between the number of restored resistors, namely n_sh, and of broken ones n_b, $\eta ={n}_{{\rm{sh}}}/{n}_{{\rm{b}}}$. Thus, for η = 0 the case without self-healing is recovered, and for η = 1, all fuses are repaired. In our model, the self-healing rate η corresponds to the probability that a broken resistor will be repaired.

In [23], link restoration was assumed instantaneous, while in this paper we assume that there can be a delay in the regeneration of a broken resistor, i.e. broken links are restored after a certain amount of time. Since the RFM is a quasistatic model, there is no real time evolution, but the applied voltage can be considered equivalent to time: if a time unit is associated with a voltage step, simulations with a linearly increasing voltage correspond to a linearly increasing time (in the RFM, any rupture and current redistribution is considered instantaneous with respect to the voltage step). Thus, a delay time can be expressed in voltage units of the simulation. If a resistor is broken at voltage V⁰ and its regeneration delay is V_sh, it is restored when the applied voltage is V⁰ + V_sh.

In a real life phenomenon, this is connected to a characteristic time of the self-healing process. For example, in some applications, microcracks are repaired by a healing agent contained in small capsules inside the material itself. In this case, the delay time would be the characteristic time scale for the healing agent to polymerise and close the microcraks. This time scale is present in principle in all self-healing processes, so it can be added phenomenologically in the model from a generic point of view.

We assign the delay of each resistor by extracting it from a Weibull distribution W(x, Δ_sh, k_sh), where the scale parameter Δ_sh is the characteristic delay and k_sh is the shape parameter. When a resistor is restored, its conductivity is set again to c = 1 and a new breaking threshold must be extracted from the Weibull distribution of the thresholds W(x, λ^', k), but the scale parameters must be rescaled so that the resistor is not instantaneously burnt. Thus, we set ${\lambda }^{{\prime} }=\lambda +{I}_{j}^{0}$, where I_j⁰ is the current flowing through it at the moment of the previous breakage. In this way, the global statistical properties of the resistors are unchanged by self-healing. A resistor can be broken and restored any number of times. The rest of the algorithm works as in the case without self-healing.

In this paper, we focus on the effects of the delay in the case η = 1, i.e. all resistors can potentially be restored. If not explicitly stated, all simulations are therefore performed with η = 1. The self-healing delay, expressed in voltage units, can be compared with the average value of V_max of the system. For ${\lambda }_{{\rm{sh}}}\ll \langle {V}_{\max }\rangle$ the case of 'fast' self-healing occurs, i.e. a broken resistor has a high probability of being restored before the failure of the system, and vice versa, for ${\lambda }_{{\rm{sh}}}\gtrsim \langle {V}_{\max }\rangle$ the case of 'slow' self-healing occurs. We will explore both of these limits. It is convenient to define the healing rapidity ${\rm{\Gamma }}\equiv 1/{{\rm{\Delta }}}_{{\rm{sh}}}$, so that the case without self-healing is obtained for Γ = 0.

In [23], an analytical calculation was performed to describe the evolution of damage in a FBM with self-healing in the limit of the Daniel's theory [19, 20], and to provide an approximate expression for its strength. The resulting expressions are approximated, since no account for damage or healing localisation is present, but they still provide a useful analytical benchmark. According to [23], the relation between stress σ (current I) and strain (voltage V) can be written as:

$$\begin{eqnarray}&&\sigma =E\ \epsilon \ {{\rm{e}}}^{(\eta -1){\left(\epsilon /{\epsilon }_{0}\right)}^{m}},\end{eqnarray} \tag{ 1 }$$

where ₀ and m are the scale and the shape parameter of the Weibull threshold distribution, respectively, and E is the Young's modulus. In the presence of a delay, a modified version of equation (1) can be provided if an effective strain-dependent self-healing rate η_eff is adopted in place of η, accounting for an average characteristic time delay Γ⁻¹ in the healing process. In particular, results are fitted with a good approximation by the phenomenological law:

$$\begin{eqnarray}&&{\eta }_{{\rm{eff}}}(\epsilon )\sim \eta \displaystyle \frac{{\left({\rm{\Gamma }}\epsilon \right)}^{m}}{1+{\left({\rm{\Gamma }}\epsilon \right)}^{m}}\end{eqnarray} \tag{ 2 }$$

which provides a simple expression with the correct asymptotic behaviour: the case Γ = 0, i.e. an infinite time delay, is equivalent to the case without self-healing, and the case ${\rm{\Gamma }}\to \infty$ is equivalent to an instantaneous self-healing as in [23]. For intermediate values of the delay, the maximum strength is reduced with respect to the instantaneous self-healing.

The maximum voltage and the peak current both increase with self-healing. This is shown in figure 2 by the characteristic I(V) curve for various self-healing rates. In particular, by increasing the value of Γ, the curves approach the ideally linear case I ∝ V, in which the last catastrophic fracture event, leading to the failure of the system, occurs immediately after the peak current and involves a large avalanche of ruptures. Thus, a signature of the self-healing process in the RFM is an increase in maximum strength and ultimate strain, but also a vanishing plastic phase and perfectly brittle behaviour with a catastrophic final event.

Zoom In Zoom Out Reset image size

To be more quantitative, in figure 3(a) we report the behaviour of the mean total damage $\langle d\rangle$ of the system as a function of applied voltage. In standard RFM simulations, the damage approximately increases linearly up to the final larger event. For slow self-healing values, i.e. ${{\rm{\Delta }}}_{{\rm{sh}}}\gtrsim \langle {V}_{\max }\rangle$, the curve lies over the standard case. For fast self-healing the slope of the curve progressively decreases and the spike due to the last catastrophic event increases, and for the fastest self-healing case the spread of damage is inhibited up to failure. In figure 3(b), we report the statistics of the number of broken fuses and restored ones at the end of the test. By increasing Γ, the number of broken fuses decreases and most of them are fractured at the last rupture event. These results suggest that with a fast self-healing rate the system damage may appear limited despite a catastrophic rupture avalanche is imminent.

Zoom In Zoom Out Reset image size

The abrupt failure of the system is not due to a weakness of the restored links, since on average only about ten percent of the final broken links have been previously restored. Links are repaired in an unstressed state, so that they have larger thresholds. Thus, we must conclude that with a self-healing delay Γ ≳ 2, the system may display a limited damage, but be in a critical state in which a single rupture can trigger the total failure.

We observe also that an imminent failure cannot be easily predicted from the statistical distribution of the currents. We define the effective current as sum over all unbroken links of the local currents normalised by their threshold values, i.e with the notation of section 2 ${I}_{{\rm{eff}}}={({n}_{{\rm{tot}}}-{n}_{{\rm{b}}})}^{-1}{\sum }_{j}I(j)/x(j)$. This is an average measure of the imminence to the failure of the links. Another observable is also the fraction of links whose ratio I(j)/x(j) is above the percentage p, namely r_p, so that, for example, r_0.95 indicates a link current of 95% of its threshold. This is an estimate of how many links are very close to failure. For both observables, results before the failure in the presence of self-healing display only a small percentage increase with respect to the standard case, as shown in figure 4, so that apparently these kinds of observables are not effective in indicating that the system is in a critical state close to failure. Average global quantities cannot capture this phenomenon because, in this case, the fracture propagates due to local current spikes (i.e. stress concentrations) around links that are progressively broken by the final avalanche. We have verified that this behaviour is not an artifact due to the choice of the elementary step ΔV. Thus, we must conclude that this is a critical feature of the system for large self-healing values.

Zoom In Zoom Out Reset image size

The scaling law of the peak current with the total size of the lattice in RFM is consistent with $\langle {I}_{\max }\rangle ={c}_{1}{L}^{\alpha }+{c}_{2}$, where c₁ and c₂ are constants and the exponent estimate is α ≃ 0.96 [46], obtained both with diamond and triangular lattices. Our data in the standard case, obtained using a square lattice in the range 12 ≤ L ≤ 160, leads to α = 0.954(1), which is consistent with the value of the literature. The other constants are c₁ = 0.308(2) and c₂ = 0.487(10).

A similar scaling law holds for the maximum voltage, which may be associated with the average lifetime of the system. For 30 ≤ L ≤ 160, we find a scaling law $\langle {V}_{\max }\rangle ={c}_{3}{L}^{\beta }+{c}_{4}$, with c₃ = 0.39(1), c₄ = 2.2(1) and exponent β = 0.941(8). Tables of datasets used for these fits are reported in the tables A1–A3. These results are stable, i.e. removing the larger lattices from the fitted range does not modify the results. This demonstrates that finite size effects are negligible.

In the presence of self-healing, these scaling laws are modified and cannot be fitted with a simple power-law behaviour. In figure 5, we report the maximum current increase as a function of Γ for different lattice sides L. For slow self-healing values, the increase of the maximum current with respect to the case without self-healing follows a power-law ${\rm{\Delta }}I\equiv \langle {I}_{\max }{(\eta )\rangle -\langle {I}_{\max }(0)\rangle \sim ({\rm{\Gamma }}L)}^{a}$, where the exponent value is a ≃ 2. This power-law behaviour in the presence of self-healing is found in the approximated theory [23, 47], and the quadratic behaviour is consistent with equations (1), (2) for small Γ. The collapse on the same curve for all L is clear from Γ L ≲ 2 . This can be derived by considering that V_max ≈ 0.38L and the condition of slow self-healing is Γ V_max ≲ 1, so that combining them we find the same condition.

Zoom In Zoom Out Reset image size

For fast self-healing values, data depends on the lattice side L, which explains the modification of the scaling law of $\langle {I}_{\max }\rangle$. However, the slope of the curve is decreasing with respect to the smallest values. This implies that by reducing the delay ${{\rm{\Delta }}}_{{\rm{sh}}}$, the self-healing process is progressively less effective.

In figure 6, we report for a comparison of the maximum current as a function of Γ for various η values and L = 30. These results are qualitatively consistent with a fit derived from inserting η_eff of equation (2) into equation (1) and calculating the maximum strength. A quantitative good agreement is found for the smaller self-healing rate and Γ ≲ λ. Moreover, results for various η collapse on the same curve, at least for slow self-healing values, by plotting ΔI/η, which is captured by equation (2).

Zoom In Zoom Out Reset image size

When a resistor is burned, the currents are redistributed inside the lattice and rupture avalanches can be triggered. The number s of resistors burned at the same voltage is the size of a single avalanche event. By iterating the simulations, the distribution p(s) of the avalanche sizes can be sampled, representing the probability of having an avalanche of size s. The avalanche distribution displays a decreasing power-law behaviour with a final peak due to the last catastrophic event involving many resistors (figure 7(a)). If the final avalanche is removed, the distribution can be written in the form of a decreasing power-law with an upper cut-off:

$$\begin{eqnarray}&&p(s)={s}^{-\tau }g(s/{s}_{0}),\end{eqnarray} \tag{ 3 }$$

where τ is the power-law exponent, and s₀ is the cut-off, which can be written as a function of the lattice size as s₀ = L^D, D being the fractal dimension of the avalanche. The exponent τ can be analytically calculated in the FBM, where it assumes exactly the value τ = 5/2. In the RFM, it has been evaluated in [41], resulting in τ = 2.75 for a diamond lattice and τ = 3.05 for a triangular lattice, with fractal dimension of about D = 1.18 for both.

Zoom In Zoom Out Reset image size

A robust method to extract these results from data is suggested in [41, 48]: once the last catastrophic event is removed, the qth moment of the distribution ${M}_{q}\equiv \langle {s}^{q}\rangle$ is evaluated as a function of the lattice size. These quantities have a power-law behaviour ${M}_{q}\sim {L}^{{\sigma }_{q}}$, from which an exponent σ_q can be evaluated. If the distribution coincides with equation (3), then the following conditions should be verified: if q < τ − 1, σ_q = 0, while for q > τ − 1 σ_q = qD + D(1 − τ). Thus, from the plot of σ_q as a function of q, both τ and D can be evaluated.

Our data are consistent with these results: we obtain D = 1.07(4) and τ = 3.00(8). This is close to the value found for a triangular lattice, which is analogous to our case with a non-uniform starting current distribution on the lattice. With these results, the collapse of the avalanche distributions on the same scaling function can be obtained, as shown in the inset of figure 7(a).

In the presence of self-healing, the avalanche distribution is modified, as shown in figure 7(b) for various self-healing delay values. As observed in the previous sections, the peak due to the last catastrophic events increases with Γ, and for fast self-healing values, there are substantial modifications. If we repeat the procedure to extract τ and D, for fast self-healing values, the data rescaled do not collapse on the same curve. On the contrary, for slow self-healing there is a perfect collapse and, for example, for Γ = 0.05 we obtain D = 0.95(2) and τ = 2.86(8). In the inset of figure 7(b) we report the comparison of the avalanche distribution once the last events are removed. The slope of the distributions decreases with fast self-healing values. An estimate can be obtained by fitting the first points of the avalanche distribution, where the cut-off function is less influential, and averaging it over the lattice sides. Results are shown in table 1.

We can conclude that in general the exponent of the avalanche distribution decreases with the addition of self-healing, i.e. there is a higher probability of larger rupture avalanches. For slow self-healing values, equation (3) is still valid, while for fast self-healing values there are substantial modifications to the cut-off expression.

At the end of the simulation, the lattice is split into two regions by a winding crack that is roughly perpendicular to the flow of current, i.e. along the vertical direction. It is possible to study the 'fractality' of the crack in a quantitative way by means of different observables [49], which should give similar results provided that the crack is indeed self-affine. For comparison between different self-healing rates, we choose the crack width whose scaling law has been investigated in [41].

The first step is to obtain a single-valued function f(y) to describe the crack, where y ∈ [0, L] is the vertical coordinate of the lattice and f(y) ∈ [0, L] is the coordinate of the crack point along the horizontal direction. In our case, we identify the links belonging to the crack as those connecting points of the two split regions of the lattice. In order to obtain a single value for f(y), from this set we take only the vertical link that has minimum f(y) for each value of y. This is a conventional choice that does not affect the scaling laws.

Hence, we can define the crack width as ${w}_{l}^{2}\equiv \langle \overline{{(f(y)-\overline{f(y)})}^{2}}\rangle$, where the average $\overline{f}$ is performed over an interval of size l ≤ L and $\langle ...\rangle$ is the statistical average over all the samples. The global width $W\equiv w(L)$ is expected to scale as a power-law W ∼ L^ζ, where ζ is the fractal exponent of the crack.

In the RFM, an anomalous scaling has been observed [41], i.e. there is another local exponent ζ_loc so that for regions of size l ≪ L, the scaling law is $w(l)\sim {l}^{{\zeta }_{{\rm{loc}}}}{L}^{\zeta -{\zeta }_{{\rm{loc}}}}$, with ζ = 0.80 and ζ_loc = 0.7 for a diamond lattice. In the case without self-healing, our data are in good agreement with this. In figure 8(a), the crack width is plotted as a function of l for various lattice sizes. From the fit of the global width W we obtain ζ = 0.791(1). The exponent ζ_loc can be estimated from a fit on the region l ≪ L of the case L = 100 and we obtain ζ_loc ≃ 0.70.

Zoom In Zoom Out Reset image size

In table 2, results introducing self-healing are reported. The power-law behaviour of the crack width is confirmed except for the largest values of Γ, whose goodness of fit decays in the same range of L ∈ [12, 100]. This can be explained by the larger rupture avalanche at the system failure, so that the crack does not propagate as in the standard case, rather it is formed almost simultaneously. In any case, the ζ exponent increases, since the lattice regions with restored links have smaller rupture probability, so that the final crack propagates around them with a winding path. Thus, on large lattices there is also a larger probability that the winding path amplitudes increase and, in other words, that the crack width increases with L at greater rate.

On the contrary, we do not observe significant variations of the local exponent ζ_loc, since the previous effect is less influential on local scales. The invariance of the local scaling properties is confirmed by another observable, i.e. the total crack length l_c, reported in figure 8(b). This can be calculated as the number of links connecting points belonging to the two split parts of the lattice. In the standard case, we find a power law-behaviour l_c = l₁L^δ + l₀ , with δ = 1.16(1). This is another measure of crack tortuosity, since it indicates super-linear length increase of a one-dimensional curve with system size. In the self-healing case, the value of total crack length increases, but the scaling exponent δ does not exhibit substantial modifications. This implies that the tortuosity of the crack on local scales is not affected significantly by self-healing.