Thermally activated intermittent dynamics of creeping crack fronts along disordered interfaces
Microstructural controls of anticrack nucleation in highly porous brittle solids
Jonas Ritter, Henning Löwe & Johan Gaume
Nucleation fronts ignite the interface rupture that initiates frictional motion
Shahar Gvirtzman & Jay Fineberg
Phase-field framework with constraints and its applications to ductile fracture in polycrystals and fatigue
Fei Xue, Tian-Le Cheng, … You-Hai Wen
Learning local, quenched disorder in plasticity and other crackling noise phenomena
Stefanos Papanikolaou
Variable-order fracture mechanics and its application to dynamic fracture
Sansit Patnaik & Fabio Semperlotti
Solitonic State in Microscopic Dynamic Failures
H. O. Ghaffari, W. A. Griffith & M. Pec
Velocity correlated crack front and surface marks in single crystalline silicon
Lv Zhao, Didier Bardel, … Daniel Nelias
Time limited self-organised criticality in the high rate deformation of face centred cubic metals
Lewis Lea, Lawrence Brown & Andrew Jardine
Plastic flow anisotropy drives shear fracture
A. Amine Benzerga, Nithin Thomas & Joshua S. Herrington
Tom Vincent-Dospital1,2,
Alain Cochard1,
Stéphane Santucci3,4,
Knut Jørgen Måløy2 &
Renaud Toussaint1,2
Scientific Reports volume 11, Article number: 20418 (2021) Cite this article
Condensed-matter physics
Statistical physics, thermodynamics and nonlinear dynamics
We present a subcritical fracture growth model, coupled with the elastic redistribution of the acting mechanical stress along rugous rupture fronts. We show the ability of this model to quantitatively reproduce the intermittent dynamics of cracks propagating along weak disordered interfaces. To this end, we assume that the fracture energy of such interfaces (in the sense of a critical energy release rate) follows a spatially correlated normal distribution. We compare various statistical features from the obtained fracture dynamics to that from cracks propagating in sintered polymethylmethacrylate (PMMA) interfaces. In previous works, it has been demonstrated that such an approach could reproduce the mean advance of fractures and their local front velocity distribution. Here, we go further by showing that the proposed model also quantitatively accounts for the complex self-affine scaling morphology of crack fronts and their temporal evolution, for the spatial and temporal correlations of the local velocity fields and for the avalanches size distribution of the intermittent growth dynamics. We thus provide new evidence that an Arrhenius-like subcritical growth is particularly suitable for the description of creeping cracks.
In the physics of rupture, understanding the effects that material disorder has on the propagation of cracks is of prime interest. For instance, the overall strength of large solids is believed to be ruled by the weakest locations in their structures, and notably by the voids in their bulk samples1,2. There, cracks tend to initiate as the mechanical stress is concentrated. A growing focus has been brought on models in which the description of the breaking matrix remains continuous (i.e., without pores). There, the material disorder resides in the heterogeneities of the matrix3,4,5,6,7,8,9. The propagation of a crack is partly governed by its spatial distribution in surface fracture energy, that is, the heterogeneity of the energy needed to generate two opposing free unit surfaces in the continuous matrix10, including the dissipation processes at the tip11. From this disorder, one can model a rupture dynamics which holds a strongly intermittent behaviour, with extremely slow (i.e., pinned) and fast (i.e., avalanching) propagation phases. In many physical processes, including12,13,14 but not limited15,16,17,18 to the physics of fracture, such intermittency is referred to as crackling noise19,20. In the rupture framework, this crackling noise is notably studied to better understand the complex dynamics of geological faults21,22,23,24,25, and their related seismicity.
Over the last decades, numerous experiments have been run on the interfacial rupture of oven-sintered acrylic glass bodies (PMMA)26,27,28. Random heterogeneities in the fracture energy were introduced by sand blasting the interface prior to the sintering process. An important aspect of such experiments concerns the samples preparation, which allows to constrain the crack to propagate along a weak (disordered) plane. It simplifies the fracture problem, leading to a negligible out-of plane motion of the crack front. This method has allowed to study the dynamics of rugous fronts, in particular because the transparent PMMA interface becomes more opaque when broken. Indeed, the generated rough air-PMMA interfaces reflect more light, and the growth of fronts can thus be monitored.
Different models have successfully described parts of the statistical features of the recorded crack propagation. Originally, continuous line models4,5,20,29 were derived from linear elastic fracture mechanics. While they could reproduce the morphology of slow rugous cracks and the size distribution of their avalanches, they fail to account for their complete dynamics and, in particular, for the distribution of local propagation velocity and for the mean velocity of fronts under different loading conditions. Later on, fiber bundle models were introduced6,30,31, where the fracture plane was discretized in elements that could rupture ahead of the main front line, allowing the crack to propagate by the nucleation and the percolation of damage. The local velocity distribution could then be reproduced, but not the long term mean dynamics of fronts at given loads. One of the most recent models (Cochard et al.8) is a thermally activated model, based on an Arrhenius law, where the fracture energy is exceeded at subcritical stresses due to the molecular agitation. It contrasts to other models that are strictly threshold based (the crack only advances when the stress reaches a local threshold, rather than its propagation being subcritical). A notable advantage of the subcritical framework is that its underlying processes are, physically, well understood, and Arrhenius-like laws have long shown to describe various features of slow fracturing processes26,32,33,34,35,36. In particular, this framework has proven to reproduce both the mean behaviour of experimental fronts37 (i.e., the average front velocity under a given load) and the actual distributions of propagation velocities along these fronts8, whose fat-tail is preserved when observing cracks at different scales38. It has recently been proposed39,40 that the same model might also explain the faster failure of brittle matter, that is, the dramatic propagation of cracks at velocities close to that of mechanical waves, when taking into account the energy dissipated as heat around a progressing crack tip. Indeed, if fronts creep fast enough, their local rise in temperature becomes significant compared to the background one, so that they can avalanche to a very fast phase, in a positive feedback loop39,40.
Here, we only consider slow fronts (i.e., fronts that creep slowly enough so that their temperature elevation is supposed to remain negligible). Building on the work of Cochard et al.8, we study various statistical features that can be generated with this Arrhenius-based model (re-introduced in the "Propagation model" section), when simulating the rupture of a disordered interface. By comparing these features to those reported for the PMMA experiment by Tallakstad et al.28,38, Santucci et al.27 and Maløy et al.26, we show a strong match to the experimental data for many of the scaling laws describing the fracture intermittent dynamics, including the growth of the fracture width ("Growth exponent and fracture energy correlation length" section), its distribution in local propagation velocity ("Local velocity distribution and fracture energy standard deviation" section), the correlation of this velocity in space and time ("Local velocities correlations" section), the size of the propagation avalanches ("Avalanches size and shape" section) and the front Hurst exponents ("Front morphology" section). We hence re-enforce the relevance of simple thermodynamics coupled with elasticity in the description of material failure.
Propagation model
We consider rugous crack that are characterised by a varying and heterogeneous advancement a(x, t) along their front, x being the coordinate perpendicular to the average crack propagation direction, a the coordinate along it, and t being the time variable (see Fig. 1). At a given time, the velocity profile along the rugous front is modelled to be dictated by an Arrhenius-like growth, as proposed by Cochard et al.8:
$$\begin{aligned} V(x,t) = V_0 \text { min}\left[ \exp \left( {-\frac{\alpha ^2[G_c(x,a)-G(x,t)]}{k_B T_0}} \right) , 1\right] , \end{aligned}$$
where \(V(x,t)=\partial a(x,t)/\partial t\) is the local propagation velocity of the front at a given time and \(V_0\) is a nominal velocity, related to the atomic collision frequency41, which is typically similar to the Rayleigh wave velocity of the medium in which the crack propagates42. The exponential term is a subcritical rupture probability (i.e., between 0 and 1). It is the probability for the rupture activation energy (i.e., the numerator term in the exponential) to be exceeded by the thermal bath energy \(k_B T_0\), that is following a Boltzmann distribution41. The Boltzmann constant is denoted \(k_B\) and the crack temperature is denoted \(T_0\) and is modelled to be equal to a constant room temperature (typically, \(T_0=298\,\) K). Using this constant temperature corresponds to the hypothesis that the crack is propagating slowly enough so that no significant thermal elevation occurs by Joule heating at its tip (i.e., as inferred by Refs.39,40). Such propagation without significant heating is notably believed to take place in the experiments by Tallakstad et al.28 that we here try to numerically reproduce, and whose geometry is shown in Fig. 1. Indeed, their reported local propagation velocities V did not exceed a few millimetres per second, whereas a significant heating in acrylic glass is only believed to arise for fractures faster than a few centimetres per second40,44. See the supplementary information for further discussion on the temperature elevation.
(Left): Separation of two rugous and sintered PMMA plates, as reported by Tallakstad et al.28 (side view). The rugosity of the (quasi-plane) interface is here massively exaggerated (the plates are centimetres thick while the standard deviation in the interface topography is less than a micrometer43). A local position of the front has an advancement a(x, t) and a velocity V(x, t). The out of frame coordinate is x and t is the time variable. (Right): top view, showing the crack font roughness, which arises from the disorder in the interface's fracture energy.
In Eq. (1), the rupture activation energy is proportional to the difference between an intrinsic material surface fracture Energy \(G_c\) (in J m\(^{-2}\)) and the energy release rate G at which the crack is mechanically loaded, which corresponds to the amount of energy that the crack dissipates to progress by a given fracture area. As the front growth is considered subcritical, we have \(G<G_c\). We here model the fracture energy \(G_c\) to hold some quenched disorder that is the root cause for any propagating crack front to be rugous and to display an intermittent avalanche dynamics. This disorder is hence dependent on two position variables along the rupture interface. For instance, at a given front advancement a(x, t), one gets \(G_c=G_c(x,a)\). The coefficient \(\alpha ^2\), in Eq. (1), is an area which relates the macroscopic G and \(G_c\) values to, respectively, the microscopic elastic energy \(U=\alpha ^2 G\) stored in the molecular bonds about to break, and to the critical energy \(U_c=\alpha ^2 G_c\) above which they actually break. See Vanel et al.36, Vincent-Dospital et al.40 or the supplementary information for more insight on the \(\alpha ^2\) parameter, which is an area in the order of \(d_0^3/l\), where \(d_0\) is the typical intra-molecular distance and l is the core length scale limiting the stress divergence at the crack tip.
Finally, the average mechanical load that is applied on the crack at a given time is redistributed along the evolving rugous front, so that \(G=G(x,t)\). To model such a redistribution, we here use the Gao and Rice3 formalism, which integrates the elastostatic kernel along the front:
$$\begin{aligned} G(x,t) = \overline{G}(t) \left[ 1 - \frac{1}{\pi } \text {PV} \int _{-\infty }^{+\infty } \frac{\partial a(x',t) / \partial {x'}}{x-x'} \mathrm {d}x'\right] . \end{aligned}$$
In this equation, \(\overline{G}\) is the mean energy release rate along the front and \(\text {PV}\) stands for the integral principal value. We, in addition, considered the crack front as spatially periodic, which is convenient to numerically implement a spectral version of Eq. (2)45 as explained by Cochard et al.8.
Equations (1) and (2) thus define a system of differential equations for the crack advancement a, which we have then solved with an adaptive time step Runge-Kutta algorithm46, as implemented by Hairer et al.47. The complete code for the crack simulation is available as a Software Heritage archive48. Further details on the code can be obtained by contacting the authors.
In this section, we discuss the main principles we have used in choosing the numerical accuracy of our solver. The related parameters are illustrated in Fig. 2.
Illustration of the discretization principles and of the solver and observation grids. Three crack fronts at three successive times are shown, over which the parameters discussed in the "Discretization" section are defined.
In attempting to correctly reproduce the experimental results of Tallakstad et al.28, this solver needs to use space and time steps, here denoted \(\Delta x_s\) and \(\Delta t_s\), at least smaller than those on which the experimental fronts were observed and analysed. Thus, \(\Delta x_s\) needs to be smaller than the experimental resolutions in space (the camera pixel size) \(\Delta x=\Delta a\) of about 2 to \(5\,\upmu\)m and \(1/\Delta t_s\) needs to be higher than the experimental camera frame rate \(1/\Delta t\). This frame rate was set by Tallakstad et al.28 to about \((100\overline{V})/\Delta x\), where \(\overline{V}\) is the average front velocity of a given fracture realisation. The propagation statistics of our simulated fronts, henceforward shown in this manuscript, have, for consistency, always been computed on scales comparable to the experimental \(\Delta x\), \(\Delta a\), \(\Delta t\) steps. Thus, as \(\Delta x_s<\Delta x\) and \(\Delta t_s<\Delta t\), we have first decimated the dense numerical outputs on the experimental observation grid, by discarding smaller time scales and by averaging smaller space scales to simulate the camera frame rate and pixel size.
As the camera resolution was 1024 pixels, the lengths L of the crack segments that Tallakstad et al.28 analysed were \(1024\Delta _x=3\) to 7 mm long, and we have then analysed our numerical simulations on similar front widths. Yet, these simulations were priorly run on longer front segments, \(L_s>L\), in order to avoid any possible edge effects in the simulated crack dynamics (for instance in the case where L would not be much bigger than the typical size of the \(G_c\) quenched disorder).
Overall, we have checked that the numerical results presented henceforward were obtained using a high enough time and space accuracy for them to be independent of the discretization (as is shown in the supplementary information).
Physical parameters values
For the model dynamics to be compared to the experiments28, one must also ensure that the \(V_0\), \(\alpha\), \(T_0\), G and \(G_c\) parameters are in likely orders of magnitude.
As \(V_0\) is to be comparable to the Rayleigh velocity of acrylic glass, we have here used \(1\,\hbox {km}\,\hbox {s}^{-1}\)49. Lengliné et al.37 furthermore estimated the ratio \(\alpha ^2/(k_B T_0)\) to be about \(0.15\,\hbox {m}^{2}\,\hbox {J}^{-1}\) and they could approximate the quantity \(V_0\exp (-\alpha ^2 \overline{G_c}/[k_B T_0])\) to about \(5\times 10^{-14}\,\hbox {m}\,\hbox {s}^{-1}\), where \(\overline{G_c}\) is the average value of \(G_c\) in the rupturing interface. With our choice on the value of \(V_0\), we then deduce \(\overline{G_c}\sim 250\,\hbox {J}\,\hbox {m}^{-2}\) (note that the trade-off between \(V_0\) and \(G_c\) should be kept in mind when comparing our results with those by Cochard et al.8, as both papers use a different \(V_0\)). The value thus inverted for the fracture energy (\(250\,\hbox {J}\,\hbox {m}^{-2}\)), that is to represent the sintered PMMA interfaces, is logically smaller but comparable to that inferred by Vincent-Dospital et al.40 for the rupture of bulk PMMA (about \(1300\,\hbox {J}\,\hbox {m}^{-2}\)). Qualitatively, the longer the sintering time, the closer one should get from such a bulk toughness, but the less likely an interfacial rupture will occur.
Experiments in two different regimes were run28: a forced one where the deflection of the lower plate (see Fig. 1) was driven at a constant speed, and a relaxation regime, where the deflection was maintained constant while the crack still advances. In both scenarii, the long term evolution of the average load \(\overline{G}(t)\) and front position \(\overline{a}(t)\) was shown8,37 to be reproduced by Eq. (1). In the case of the experiments of Tallakstad et al.28, the intermittent dynamics measured in the two loading regimes were virtually identical. Such similarity likely arises from the fact that the avearge load \(\overline{G}\) was, in both cases, computed to be almost constant over time, in regard to the spatial variation in G, described by Eq. (2) (see the supplementary information). Here, we will then consider that the crack is, in average along the front, always loaded with the same intensity (i.e., \(\overline{G}(t)=\overline{G}\)).
The actual value of \(\overline{G}\), together with the average surface fracture energy of the medium \(\overline{G_c}\), then mainly controls the average crack velocity \(\overline{V}\). This average velocity was investigated over five orders of magnitude in Ref.28, from 0.03 to \(140\,\upmu \hbox {m}\,\hbox {s}^{-1}\), which, in our formalism, shall correspond to values of \((\overline{G_c}-\overline{G}\)) between 145 and \(85\,\hbox {J}\,\hbox {m}^{-2}\), respectively, which is actually consistent with the values of \(\overline{G}\) measured by Lengliné et al.37 for cracks propagating at similar speeds. The intermittency of the crack motion was experimentally inferred to be independent on \(\overline{V}\) and we show, in the supplementary information, that it is also the case in our simulations. The velocity variation along the front shall then only arise from the disorder in \(G_c\) and from the related variations of G due to the roughness of the crack front. Further in this manuscript, we will use \(\overline{G}=120\,\hbox {J}\,\hbox {m}^{-2}\), which corresponds to an average propagation velocity of about \(1.5\,\upmu \hbox {m}\,\hbox {s}^{-1}\).
Heterogeneous fracture energy
Of course, the actual surface fracture energy field in which the rupture takes place will significantly impact the avalanches dynamics and the crack morphology. Such a field is yet a notable unknown in the experimental set-up of Tallakstad et al.28, as their interface final strength derived from complex sand blasting and sintering processes. Although these processes were well controlled, so that the rough rupture experiments were repeatable, and although the surfaces prior to sintering could be imaged43, the actual resulting distribution in the interface cohesion was not directly measurable. While this is, to some extent, a difficulty in assessing the validity of the model we present, we will here show that a simple statistical definition of \(G_c\) is enough to simulate most of the avalanches statistics.
We will indeed consider a normally distributed \(G_c\) field around the mean value \(\overline{G_c}\) with a standard deviation \(\delta G_c\) and a correlation length \(l_c\). Such a landscape in \(G_c\) is shown in Fig. 3a, and we proceed to discuss the chosen values of \(\delta G_c\) and \(l_c\) in the "Growth exponent and fracture energy correlation length" and "Local velocity distribution and fracture energy standard deviation" sections.
Growth exponent and fracture energy correlation length
Among the various statistical features studied by Tallakstad et al.28, was notably quantified the temporal evolution of their fracture fronts morphology. It was interestingly inferred that the standard deviation of the width evolution of a crack front h scales with the crack mean advancement:
$$\begin{aligned} \text {rms}\big (h(t)\big ) = \sqrt{{<}{h(t)}^2{>}_{x,t_0}} \propto \big (\overline{V}t\big )^{\beta _G}. \end{aligned}$$
In this equation, x is a given position along the front, t is a time delay from a given reference time \(t_0\), and h writes as
$$\begin{aligned} h_{x,t_0}(t) = \big [a(x,t_0+t)-\overline{a}(t_0+t)\big ]-\big [a(x,t_0)-\overline{a}(t_0)\big ], \end{aligned}$$
\(\overline{a}\) being the average crack advancement at a given time. To mitigate the effect of the limited resolution of the experiments and obtain a better characterization of the scaling of the interfacial fluctuations on the shorter times, we computed the subtracted width,
$$\begin{aligned} W(t) = \sqrt{\text {rms}(h(t))^2-\text {min}(\text {rms}(h(t)))^2}, \end{aligned}$$
as proposed in Barabasi and Stanley50, and done by Tallakstad et al.28 (whose experiments we here reproduce) and Soriano et al.51.
(a) Normal distribution of the fracture energy \(G_c\) considered for the simulations. The average value is \(\overline{G_c}=250\,\hbox {J}\,\hbox {m}^{-2}\), with a standard deviation \(\delta G_c=35\,\hbox {J}\,\hbox {m}^{-2}\) and a correlation length \(l_c=50\,\upmu\)m. The three lines are the modelled propagating front at three different times \(t_1<t_2<t_3\), using Eqs. (1) and (2). (b) A crack front reported by Tallakstad et al.28 (Fig. 3 of the experimental paper), plotted on the same spatial scales. (c,d) Local velocity maps V(x, a) in the space–space domain (c) and V(x, t) in the space-time domain (d) for a modelled crack propagating in this \(G_c\) landscape. Both maps are shown with the same color scale and they are computed on a resolution similar to that of the experiments by Tallakstad et al.28, using the waiting time matrix (see text for details). The velocity are plotted related to the average crack velocity \(\overline{V}=1.5\,\upmu \hbox {m}\,\hbox {s}^{-1}\). All parameters used to run the corresponding simulation are summarised in Table 1.
The scaling exponent \(\beta _G\) is referred to as the growth exponent, and we will here show how it allows to deduce a typical correlation length for the interface disorder. Indeed, \(\beta _G\) was measured to be \(0.55\pm 0.08\) by Tallakstad et al.28. This value is close to 1/2, that is, consistent with an uncorrelated growth process (e.g.,50), such as simple diffusion or Brownian motion. We thus get a first indication on the disorder correlation length scale \(l_c\). To display an uncorrelated growth when observed with the experimental resolution (\(\Delta x\sim 3\,\upmu\)m), the fronts likely encountered asperities whose size was somewhat comparable to this resolution. Indeed, if these asperities in \(G_c\) were much bigger, the growth would be perceived as correlated. By contrast, if they were much smaller (orders of magnitude smaller), the rugosity of the front would not be measurable, as only the average \(\overline{G_c}\) over an observation pixel would then be felt. Furthermore, and as shown in Fig. 4a, the exponent \(\beta _G\) was observed on scales \((\overline{V}t)\) up to \(100\,\upmu\)m, above which W stabilised to a plateau value of about \(30\,\upmu\)m. A common picture is here drawn, as both this plateau value and the typical crack propagation distance at which it is reached are likely to also be correlated with \(l_c\), as the front is to get pinned on the strongest asperities at this scale.
(a) Standard deviation of the width evolution of the crack front as a function of the mean crack advancement, as defined by Eqs. (3) to (5) for the chosen simulation (plain points) and for the experiments28 (hollow stars) (out of Fig. 8, Expt. 5 of the experimental paper). The continuous line has a slope 0.6, close to that of the experimental points: \(\beta _G\sim 0.55\). The numerical \(\beta _G\), obtained with a linear root mean square fit of the growth of W, is estimated as \(\beta _G=0.60\pm 0.05\). The dashed lines mark the observation scale \(\Delta x\), corresponding to the experimental camera pixel size, and the chosen correlation length for the simulation \(l_c=50\,\upmu\)m. (b) The same width function for simulations with different correlation lengths \(l_c\). The rest of the parameters are as defined in Table 1. The slope and plateau of the experimental data (shown in (a)) is marked by the dashed line for comparison.
From all these clues, we have considered, in our simulations, the correlation length of the disorder to be about \(l_c=50\,\upmu\)m, and we show in Fig. 4a that it allows an approximate reproduction of the front growth exponent and of the plateau at high \(\overline{V}t\). Note that the accuracy reported for the exponents in this manuscript is estimated by fitting various portions of the almost linear data points and reporting the dispersion of the thus inverted slopes. In Fig. 4b, we also show how varying \(l_c\) impacts W, and, in practice, we have chosen \(l_c\) by tuning it when comparing these curves to the experimental one. Noteworthily, the thus chosen \(l_c\) is in the lower range of the size of the blasting sand grains (\(50-300\,\upmu\)m) that were used28 to generate the interface disorder. It is also comparable to the correlation length of the blasting induced topographic anomalies \(\sim 18\,\upmu\)m on the post-blasting/pre-sintering PMMA surfaces, as measured by Santucci et al.43 by white light interferometry.
Local velocity distribution and fracture energy standard deviation
While the crack advances at an average velocity \(\overline{V}\), the local velocities along the front, described by Eq. (1), are, naturally, highly dependent on the material disorder: the more diverse the met values of \(G_c\) the more distributed shall these velocities be.
Maløy et al.26 and Tallakstad et al.28 inferred the local velocities of their cracks with the use of a so-called waiting time matrix. That is, they counted the number of discrete time steps a crack would stay on a given camera pixel before advancing to the next one. They then deduced an average velocity for this pixel by inverting this number and multiplying it by the ratio between the pixel size and the time between two pictures: \(\Delta a/\Delta t\). Such a method, that provides a spatial map V(x, a), was applied to our simulated fronts, and we show this V(x, a) map in Fig. 3c. As to any time t corresponds a front advancement a(x, t) (recorded with a resolution \(\Delta a\)), an equivalent space-time map V(x, t) can also be computed, and it is shown in Fig. 3b. The experimental report28 presented the probability density function of this latter (space-time) map \(P(V/\overline{V})\), and it was inferred that, for high values of V, the velocity distribution scaled with a particular exponent \(\eta = 2.6 \pm 0.15\)28,38 (see Fig. 5a). That is, it was observed that
$$\begin{aligned} P\left( {V}/{\overline{V}}\right) \propto \left( {V}/{\overline{V}}\right) ^{-\eta }. \end{aligned}$$
Cochard et al.8, who introduced the model that we here discuss, inferred that the \(\eta\) exponent was mainly depending on \(\alpha ^2(\delta G_c)^2/[k_B T_0(\overline{G_c})^2]\). Truly, a more comprehensive expression could also include other quantities, such as \(V_0\) or \(l_c\). Yet, as all other parameters have now been estimated, we can deduce \(\delta G_c\) by varying it to obtain \(\eta \sim 2.6\). We show, in Fig. 5b, how varying \(\delta G_c\) impacts \(P(V/\overline{V})\) and \(\eta\). We found \(\delta G_c\sim 35\,\hbox {J}\,\hbox {m}^{-2}\). In Fig. 5a, we show the corresponding velocity distribution for a simulation run with this parameter, together with that from Tallakstad et al.28, showing a good match. Note that the ability of the model to reproduce the local velocity distribution was already shown by Cochard et al.8, and this figure mainly aims at illustrating our calibration of the fracture energy field. The model we present is also slightly different to that of Cochard et al., as the interface fracture energy is, contrarily to this previous study, now described at scales below its correlation length, similarly to the observation scale of the experiments. We here verify that the reproduction of the local velocity distribution is still valid at these small scales. Satisfyingly, the inverted value of \(\delta G_c\) is not too far from the value found by Lengliné et al.37 for their fluctuation in the mean fracture energy \(\overline{G_c}\) along their sintered plates, when studying the mean front advancement (i.e., neglecting the crack rugosity) in similar PMMA interfaces, which was about \(25\,\hbox {J}\,\hbox {m}^{-2}\).
(a) Probability density function of the local propagation velocity along a simulated front (plain points), computed from the space-time map of Fig. 3d. The experimental probability28 (out of Fig. 5, Expt. 5 of the experimental paper) is shown for comparison (hollow stars). The continuous line has a slope \(-2.6\), close to that of the experimental points. This was achieved by setting the standard deviation for the disorder in fracture energy to \(35\,\hbox {J}\,\hbox {m}^{-2}\). The numerical \(\eta\), obtained with a linear root mean square fit of the distribution tail, is estimated as \(\eta =2.6\pm 0.1\). (b) The same distribution for three three simulations with different values of \(\delta G_c\). We chose the value of \(\delta G_c\) by tuning it and fitting the slope and maximum of the experimental data, which are illustrated by the dashed lines. The rest of the parameters used in these simulations are as defined in Table 1. Note that the ability of the model to reproduce the local velocity distribution was already shown by Cochard et al.8 (see text for explanation and discussion).
Further statistics
We have now estimated the orders of magnitude of all parameters in Eqs. (1) and (2), including a likely distribution for an interface fracture energy representative of the experiments28 we aim to simulate (i.e., including \(\overline{G_c}\), \(\delta G_c\) and \(l_c\)). For convenience, this information is summarised in Table 1.
Table 1 Summary of all parameters that are considered in this manuscript.
We will now pursue by computing additional statistics of the crack dynamics to compare them to those reported by Tallakstad et al.28.
Local velocities correlations
In particular, we here compute the space and time correlations of the velocities along the front. That is, four correlation functions that are calculated from the V(x, t) and V(x, a) matrices (shown in Fig. 3) and defined as:
$$\begin{aligned} C_{ij}(\delta i)=\left\langle \frac{\left[ V(i_0+\delta i,j)-\overline{V_j}\,\right] \left[ V(i_0,j)-\overline{V_j}\,\right] }{{(\delta V_j)}^2}\right\rangle _{i_0}, \end{aligned}$$
where i and j are the variables of either V(x, t) or V(x, a) and \(\delta i\) a given i increment. \(\overline{V_j}\) is the mean of V(i, j) taken along j at a given \(i_0\). The corresponding \(\delta V_j\) is the velocity standard deviation along the same direction and for the same \(i_0\). The correlation functions hence defined are the same as those used by Tallakstad et al.28 on their own data, allowing to display a direct comparison of them in Fig. 6. A good general match is obtained.
One can notice the comparable cut-offs along the x axis (Fig. 6a,c), indicating that our chosen correlation length for the interface disorder (\(l_c\) inferred in the "Growth exponent and fracture energy correlation length" section) is a good account of the experiment. Yet, one can notice that \(C_{xt}\) (the velocity correlation along the crack front shown in Fig. 6a) is higher in the numerical case than in the experimental one. It could translate the fact that the experimental disorder holds wavelengths that are smaller than the observation scale \(\Delta x\), and that our modelled \(G_c\) distribution, where \(l_c>\Delta x\), is rather simplified.
Local velocity correlation functions in space and time as defined by Eq. (7). The plain points were computed from the simulation whose parameters are presented in Table 1 and the hollow stars are some of the experimental data points extracted from Figs. 6 and 7 of Tallakstad et al.28. In plot (a), the line overlying the numerical data set corresponds to a fit using Eq. (8). (b) is plotted in a domain that allowed a good collapse of the experimental data for many experiments28. The parameters A, \(a^*\) and \(\tau _t\) were inverted from Eq. (9), and the related fit is shown by the dashed line overlying the numerical data set. Plots (a,c,d) hold two curves for the experiments, corresponding to two distinct sets of experiments done on two different sintered PMMA bodies. (b) Shows Expt. 5 of Tallakstad et al.28.
To go further, Tallakstad et al.28 modelled \(C_{xt}\) as
$$\begin{aligned} C_{xt}(\delta x)\propto {\delta x}^{-\tau _x}\exp \left( -\frac{\delta x}{x^*}\right) , \end{aligned}$$
and inverted the values of \(\tau _x\) and \(x^*\) to, respectively, 0.53 and about \(100\,\upmu\)m. Doing a similar fit on the simulated data, we found \(\tau _x\sim 0.13\) and \(x^*\sim 94\,\upmu\)m. The related function is displayed in Fig. 6a (plain line). Our small \(\tau _x\sim 0.13\) may derive, as discussed, from the better correlation that our simulation displays at small \(\delta x\) (\(\tau _x\) may in reality tend to zero for scales smaller than those we observe) compared to that of the experiments, while the matching \(x^*\) probably relates to a satisfying choice we made for \(l_c\). Overall, the existence of a clear scaling law at small offsets, as defined by Eq. (8), is rather uncertain (see the two experimental plots in Fig. 6a) so that mainly the cut-off scale is of interest.
On the time correlation \(C_{tx}\) (Fig. 6b), one can similarly define the parameters A, \(\tau _t\) and \(a^*\) to fit Eq. (7) with a function
$$\begin{aligned} C_{tx}(\delta t)\approx A{\delta t}^{-\tau _t}\exp \left( -\frac{\overline{V}\delta t}{a^*}\right) , \end{aligned}$$
where A is a constant of proportionality. Fitting this function to Eq. (7) with a least-squares method, we found \(\tau _t\sim 0.3\) and \(a^*\sim 4.3\,\upmu\)m, and this fit is represented by the dashed line in Fig. 6b. Tallakstad et al.28 reported \(\tau _t\sim 0.43\) and \(a^*\sim 7\,\upmu\)m. Figure 6b shows the experimental and simulated correlation functions in the \(\overline{V}\delta t/a^*\)—\(C_{tx}(a^*/\overline{V})^{\tau _t}/A\) domain, as this allowed a good collapse of the data from numerous experiments28. We show that it also allows an approximate collapse of our modelled correlation on a same trend. Finally, the derived value of \(a^*\) consistently matches the apparent cut-off length in the \(C_ax\) correlation function in Fig. 6d. This length being of a magnitude similar to that of the observation scale \(\Delta a\), the crack local velocities appear uncorrelated along the direction of propagation, which is consistent with the \(\beta _G\sim 1/2\) growth exponent (e.g.,50).
Avalanches size and shape
We pursue by characterising the intermittent, burst-like, dynamics of our crack fronts and, more specifically, the avalanche (or depinning) and pinning clusters shown by the local front velocity V(x, a). We define an avalanche when the front velocity locally exceeds the mean velocity \(\overline{V}\) by an arbitrary threshold that we denote c, that is, when
$$\begin{aligned} V(x,a)>c\overline{V}. \end{aligned}$$
Similarly, we state that a front is pinned when
$$\begin{aligned} V(x,a)<\frac{\overline{V}}{c}. \end{aligned}$$
We then map, in Fig. 7, the thus defined avalanching and pinned locations of the crack. Following the analysis of Tallakstad et al.28, we compute for each of these clusters the surface S, the crossline extent \(l_x\) (that is, the maximum of a cluster width in the x direction) and the inline extent \(l_a\). The definition chosen for \(l_a\) varies for the avalanche clusters, where the maximal extent along the a direction is regarded, or the pinned one, where the mean extent along the a direction is rather used. This choice was made28 because the pinning clusters tend to be more tortuous so that their maximum span along the crack direction of propagation is not really representative of their actual extent (see Fig. 7).
Positions of the avalanches (left) and pinning locations (right) in the front local velocity map V(x, a) shown in Fig. 3c, as per Eqs. (10) and (11). Two thresholds are here used to define these maps relatively to the mean velocity: \(c=3\) and \(c=6\). The white areas are the locations of interest, of surfaces S, crossline extents \(l_x\) and inline extents \(l_a\). Bottom images: Difference in definition of \(l_a\) for the avalanche (or depinning) and pinning clusters shown in Fig. 7. For the former, \(l_a\) is the maximum extent along the a direction. For the latter it is the average width in the same direction. In both cases, \(l_x\) is the maximum extent along the x direction and S the full surface (in white) of the cluster. The square pattern marks the pixel size (\(\Delta x=\Delta a = 3\,\upmu\)m).
In Fig. 8a, we show the probability density function of the cluster surface P(S) and compare it to the experimental one. One can notice that it behaves as
$$\begin{aligned} P(S)\propto S^{-\gamma }, \end{aligned}$$
with \(\gamma =1.44\pm 0.15\). This value is comparable to the exponent inverted experimentally28, that is, \(\gamma =1.56\pm 0.04\).
(a) Probability density function of the surface of the modelled avalanche clusters (plain points) and of the modelled pinning clusters (crosses), for a threshold \(c=3\). The straight dashed line has a slope \(\gamma =1.44\), as per Eq. (12). For comparison, the hollow stars show the experimental probability density function obtained by Tallakstad et al.28 for the avalanche and pinning clusters (both are overlapping, see Fig. 10 of their manuscript). (b) Same probability density function for various c values: \(c=1.5\) (squares), \(c=3\) (plain points), \(c=6\) (stars), \(c=12\) (circles). The straight line has a slope \(\gamma =1.4\), as per Eq. (12). (c) Variation of the mean avalanche size \(\overline{S}\) as a function of the threshold c for the simulation (plain points) and the experiments (hollow stars). The modelled \(\overline{S}\) is expressed in pixels (one pixel is \(9\,\upmu \hbox {m}^{2}\)) and the experimental \(\overline{S}\) reported by Tallakstad et al.28 (in their Fig. 13) is in an arbitrary unit, so that the magnitude of both should not here be compared. We have here shifted up this experimental data set for an easier comparison with the numerical one. The straight line has a slope 0.68.
Of course, the size of the avalanche (depinning) clusters highly depends on the chosen threshold c, but we verified, as experimentally reported, that the value of \(\gamma\) inverted from the simulated data is not dependent on c, as shown in Fig. 8b. We also show, in Fig 8c, that the mean cluster size \(\overline{S}\) varies with c approximately as \(\overline{S}\propto c^{-m}\), with \(m\sim {0.68}\). This value is comparable with the experimental scaling law28 measured to be \(\overline{S}\propto c^{-0.75}\).
We also computed the probability density function of \(l_x\) and \(l_a\), that are respectively compared to their experimental equivalent in Fig. 9. These functions can be fitted with
$$\begin{aligned}&P(l_x)\propto {l_x}^{-\beta _x}, \end{aligned}$$
$$\begin{aligned}&P(l_a)\propto {l_a}^{-\beta _a}, \end{aligned}$$
and we found \(\beta _x=1.7\pm 0.2\), close to the reported experimental value28 \(\beta _x\sim 1.93\). The value we found for \(\beta _a\sim 2.2=0.2\) is also inline with that of Tallakstad et al.28, who reported \(\beta _a\sim 2.36\).
It should be noted that, while we have here fitted P(S), \(P(l_x)\) and \(P(l_a)\) with plain scaling laws (i.e., with Eqs. (12) to (14)), Tallakstad et al.28 also studied the cut-off scales above which these scaling laws vanish in the experimental data, and the dependence of these cut-off scales with the arbitrary threshold c. In our case, such scales are challenging to define, as one can for instance notice in Figs. 8 and 9, where an exponential cut-off is not obvious. This may result from a limited statistical description of the larger avalanches in our simulations. Similar cut-off scales, decreasing with increasing c should however hold in our numerical data, in order to explain the decrease of average avalanche size with c, as shown in Fig. 8c.
(a) Probability density function of the crossline extent \(l_x\) of the modelled avalanche clusters (plain points) and of the modelled pinning clusters (crosses), for a threshold \(c=3\). The straight line has a slope \(\beta _x=1.7\), as per Eq. (13). The hollow stars shows the experimental probability density function obtained by Tallakstad et al.28 for the pinning and avalanche clusters (from their Fig. 16a, inset, c = 3). (b) Probability density function of the inline extent \(l_a\) of the modelled avalanche clusters (plain points) and of the modelled pinning clusters (crosses), for a threshold \(c=3\). The two straight dashed lines have a slope \(\beta _x=2.2\), inline with that of the experimental data from Tallakstad et al.28 for the pinning clusters (hollow stars, from their Fig. 16b, inset, c = 3).
Front morphology
Finally, we show, in Fig. 10a, the relations between the clusters surface S and their linear extent \(\overline{l_x}\) and \(\overline{l_a}\). Here, \(\overline{l_x}\) and \(\overline{l_a}\) are the mean extents for all the observed clusters sharing a same surface (with the given pixel size limiting the resolution). We could fit these relations with \(\overline{l_x}\propto S^{0.77}\) and \(\overline{l_a}\propto S^{0.25}\) for the pinning clusters, and with \(\overline{l_x}\propto S^{0.64}\) and \(\overline{l_a} \propto S^{0.47}\) for the avalanches clusters. It is in qualitative agreement with the laws observed by Tallakstad et al.28: \(\overline{l_x}\propto S^{0.63}\) and \(\overline{l_a}\propto S^{0.34}\) for the pinning clusters, and \(S\propto \overline{l_x}^{0.61}\) and \(\overline{l_a}\propto S^{0.41}\) for the avalanches clusters. These exponents were experimentally reported with a \(\pm 0.05\) accuracy, and we estimated comparable error bars for the numerically derived ones. Thus, the shape of our simulated avalanches and pinned locations is rather similar to the observed experimental ones. Note that, from all the previous exponents, one can easily define H such that \(l_a\propto {l_x}^H\), and we thus have \(H_p\sim 0.25/0.77=0.32\pm 0.1\) and \(H_d\sim 0.47/0.64=0.73\pm 0.01\) for, respectively, the simulated pinning and depinning clusters (see Fig. 10b).
(a) Mean linear extents of the simulated pinning and depinning clusters as a function of cluster size. The four data sets are, from top to bottom, \(\overline{l_x}\) for the pinning clusters (hollow stars), \(\overline{l_x}\) for the avalanche clusters (crosses), \(\overline{l_a}\) for the avalanche clusters (plain points), \(\overline{l_a}\) for the pinning clusters (hollow points). The straight lines correspond to the fits described in the inset. See text for the equivalent experimental exponents. (b) Mean inline extent \(\overline{l_a}\) as a function of the mean crossline extent \(\overline{l_x}\) for the pinning and depinning clusters. The straight lines have a slope of, respectively, \(H_p=0.32\) and \(H_d=0.73\).
It was suggested5,52 that H is a good indicator of the front morphology, as the front shape is to be highly dependent on the aspect ratio of its avalanches. To verify this hypothesis, we computed the advancement fluctuation along the front \(\sigma\), that is
$$\begin{aligned} \sigma (\delta x)=\sqrt{{<}\left( a(x_0+\delta x, t)-a(x_0, t)\right) ^2{>}_{x_0,t}}. \end{aligned}$$
While this quantity was not presented by Tallakstad et al.28, it was provided by other experimental works done on the same set-up26,27, and Fig. 11a shows \(\sigma\) as reported by these authors, together with that computed in the output of our simulation. One can notice that the numerical fronts are less rugous than the experimental ones. Such a mismatch is here due to the fact that the experiment from Santucci et al., shown in Fig. 11a, had more rugous crack fronts than the one from Tallakstad et al., to which the simulation is calibrated (as shown in Fig. 4). In both cases, the data sets seem to present two self-affine behaviours (e.g.,50) with a Hurst exponent \(\zeta\) that differs at low and high length scales. Noting \(\delta x^*\) the cut-off between these length scales we indeed have:
$$\begin{aligned}&\sigma \propto {\delta _x}^{\zeta ^-}\,\,\,\text {for}\,\,\,\delta x < \delta x^*, \end{aligned}$$
$$\begin{aligned}&\sigma \propto {\delta _x}^{\zeta ^+}\,\,\,\text {for}\,\,\,\delta x > \delta x^*. \end{aligned}$$
We derived \(\zeta ^-\sim 0.68\pm 0.05\) and \(\zeta ^+=0.4\pm 0.05\) for the simulation, which compare well to the exponents that were measured experimentally, respectively, \(\zeta ^-=0.60\pm 0.05\) and \(\zeta ^+=0.35\pm 0.05\) and which are also close to the values we found for \(H_d\) and \(H_p\). The cut-off scale between the two regimes is also similar in both the experimental and numerical cases: \(\delta x^*\sim 80\,\upmu\)m, comparable to the disorder correlation length \(l_c\), and to the length scales \(x^*\), below which the local propagation velocities are correlated.
(a) Advancement fluctuation \(\sigma\) along the crack fronts, as per Eq. (15), for the simulation (plain points) and an experimental data set from Santucci et al.27 (see their Fig. 4). Different self-affine behaviours are observed above and below the \(\delta x^*\) cut-off, with comparable Hurst exponents \(\zeta\). The dashed lines mark the slopes fitted on the simulation data for the two cases. The experimental points are from an experiment different from those of Tallakstad et al.28 to which the model was calibrated. (b) Power spectra of the simulated crack advancement, averaged over 10,000 consecutive fronts. It is shown both before (raw) and after (binned) binning the fronts to the experimental camera pixel size. The difference between these two plots shows an influence of the observation scale on the small-scale study of the crack morphology. The plain line is the approximation8 from Eq. (18), which is valid between \(l_c\) and \(R_c\), where the morphology is dominated by the material quenched disorder. Note that the scaling regime for scales above \(l_c\) was already studied by Cochard et al.8, while the model match to the experiment below this cut-off scale, shown in (a), is a new result, as already discussed in the "Local velocity distribution and fracture energy standard deviation" section.
For scales above this correlation length, Cochard et al.8 showed, by analytically analysing the same model as we here study, that the front morphology is dominated by the material quenched disorder with a Hurst coefficient approximating to \(\zeta ^+=0.5\). At even larger scales, above \(R_c\sim \pi l_c \alpha ^2 \overline{G_c}/(k_B T_0)\), they also showed8 that the roughness of the simulated cracks ceases to be governed by the quenched disorder but is rather dominated by the thermal (annealed) noise, with \(\sigma\) decaying logarithmically and with a Hurst coefficient tending to \(\zeta ^\infty =0\). With our set of parameters, \(R_c\) computes to \(6\,\) mm, which is close to, yet bigger than, the total analysed length of the front. The value \(\zeta ^+\sim 0.4\), that we have here inverted, arises then likely from the transition between the two regimes, \(\zeta ^+=0.5\) and \(\zeta ^\infty =0\), as already mentioned for the experimental case, in Ref.27. In addition to a theoretical Hurst exponent \(\zeta ^+=0.5\), Cochard et al.8 computed an analytical approximation for the fronts morphology power spectrum \(P_a(\Gamma )\), for the length scales \(\Gamma\) for which the effect of the quenched disorder prevails:
$$\begin{aligned} P_a(\Gamma )\sim \left( \frac{\delta G_c}{\overline{G}}\right) ^2\frac{\Gamma R_c}{4\pi ^2}. \end{aligned}$$
We show, in Fig. 11b, how this approximation also fits the power spectra of our modelled front.
We studied an interfacial fracture propagation model, based only on statistical and subcritical physics in the sense of an Arrhenius law (Eq. (1)) and on the elastic redistribution of stress along crack fronts (Eq. (2)). Following the work of Cochard et al.8, we here showed that it allows a good representation of the intermittent dynamics of fracture in disordered media, as it approximately mimics the scaling laws dictating the propagation of experimental fronts, such as their growth exponent, their local velocity distribution and space and time correlations, the size of their avalanches and their self-affine characteristics.
To run our simulations, we had to assume a given distribution for the toughness of the rupturing interface, as this quantity is not directly measurable in the laboratory. We proposed \(G_c\) to be normally distributed with a unique correlation length and, of course, this can only be a rough approximation of the actual fracture energy obtained by Tallakstad et al.28 by sintering two sand-blasted plexiglass plates. From this approximation, could arise discrepancies between our simulations and the experiments. We have indeed shown how some of the observed exponents were strongly dependent on the definition of the material disorder. We also have assumed a perfectly elastic crack front, when the local dynamics of creeping PMMA could be visco-elastic in part, particularly below the typical length scale \(r\,\sim \,GE/\sigma _y^2 \sim 30\,\upmu\)m for plasticity around crack tips (e.g.,1) in this material, where \(\sigma _y \sim 100\) MPa is the tensile yield stress of the polymer and \(E \sim 3\) GPa its Young modulus53.
These points being stated, the vast majority of the statistical quantities that we have here studied show a good match to those from the experimental observations, so that both the considered physical model and the interface definition are likely to be relevant. A further validation of this thermally activated model could derive from the comparison of its predictions with interfacial experiments at various background temperatures \(T_0\). However, such experimental data is, to our knowledge, not yet available. Of course, some of our considered parameters (e.g., \(G_c\), \(V_0\) or \(\alpha\)) may, in practice, be temperature dependent so that a straight transposition of the model to different background temperatures could prove to be too simple. Creep experiments in bulk PMMA at various room temperatures can however be found in the literature54, where only the mean front velocity versus the mean mechanical load are measured. In this case54, it is reported that the creep dynamics is compatible with an Arrhenius-like process. By submitting many different materials to a constant load, at various temperatures, their lifetime was also shown33,36 to follow an Arrhenius law, with an energy barrier that decreases with the applied stress. These materials include metals, alloys, non-metallic crystals and polymers (and PMMA in particular).
It should be noted that, as stated in our introduction, other models have been considered to numerically reproduce the interfacial PMMA experiments, notably, a non-subcritical threshold based fluctuating line model by Tanguy et al.29, Bonamy et al.4 or Laurson et al.5,20 and a fiber bundle approach by Schmittbuhl et al.6, Gjerden et al.30 or Stormo et al.31. The present manuscript does not challenge these other models per se, but rather offers an alternative explanation to the intermittent propagation of rough cracks. The former model, the fluctuating line model4,5,20,29, considers a similar redistribution of energy release rate G as proposed in Eq. (2), but with a dynamics that is thresholded rather than following a subcritical growth law. The fronts either move forward by one pixel5 if \(G>G_c\), or with a velocity proportional4 to (\(G-G_c\)). It is completely pinned otherwise (\(V=0\) for \(G<G_c\)). While reproducing several statistical features of the experiments, this non-subcritical line propagation model does not simulate the mean propagation of cracks in various loading regimes (as done by Cochard et al.8) or the distribution in local velocity55, and, in particular, the power law tail of this distribution (i.e., Fig. 5). By contrast, the latter model6,30,31, the fiber bundle one, can reproduce this particular power law tail. It is not a line model: the interface is sampled with parallel elastic fibers breaking at a given force threshold. This threshold is less in the vicinity of the crack than away from it (it is modelled with a linear gradient), explaining why the rupture is concentrated around a defined front, and it holds a random component in order to model the quenched disorder of the interface. An advantage of the fiber bundle model is to be able to describe a coalescence of damage in front of the crack56 rather than solely describing a unique front. This could likely also be achieved in a subcritical framework, but would require to authorise damage in a full 2D plane, or require a full 3D modelling (i.e., also authorise out-of-plane damage), rather than only the modelling of a 2D front. In practice, thermal activation and damage coalescence may occur simultaneously. The observation of actual damage nucleation, in the experiments that we reproduce, has however never been obvious. Instead, the experimental fronts look rather continuous. Coalescence could yet still be at play at length scales smaller than the observation resolution. This being stated, an advantage of our model is to only rely on the experimental observations, on stress redistribution and on statistical physics. Another clear advantage of the Arrhenius based model, when compared to the other ones, is to hold a subcritical description that is physically well understood and that is a good descriptor of creep in many materials1,36. For the record, we show in Table 2 a comparison between the different exponents predicted by the three models, that all successfully reproduce some experimental observables.
Table 2 Comparison of various exponents and cut-off scales derived experimentally27,28 (Expt.) and numerically with the present, Arrhenius Based, fluctuating Line model (ABL), the Fiber Bundle model30 (FB) and the Non-Subcritical fluctuating Line model4,5 (NSL).
Note that, if linearizing Eq. (1) with a Taylor expansion around \(\overline{G_c}-\overline{G}\), that is, for propagation velocities close to the mean crack speed \(\overline{V}=V_0\exp (-\alpha ^2[\overline{G_c}-\overline{G}]/[k_B T_0])\), one obtains
$$\begin{aligned} V \sim V_\text {cst}-\frac{\alpha ^2\overline{V}}{k_B T_0}(G_c-G), \end{aligned}$$
where \(V_\text {cst}\) is a constant equal to \(\overline{V}(1+\alpha ^2[\overline{G_c}-\overline{G}]/[k_B T_0])\). This simplified form for our subcritical model is mathematically similar to that of the overcritical model (in the sense that a non zero velocity is only obtained for \(G>G_c\)) of Bonamy et al.4, where \(V=\text {max}[\mu (G_c-G),0]\) and where the coefficient of proportionality \(\mu\) was named the 'effective mobility of the crack front'. Equation (19) may give some insight in the physical meaning of \(\mu\) in this alternative model4. While the above similitude in mathematical forms may explain the obtention of some similar exponents in the dynamics of the two models (see Table 2), Eq. (19) is only a crude approximation of our highly non-linear Arrhenius formalism, which, as discussed below, allows a more exhaustive description of the experimental intermittent creep dynamics. In our simulations, the exponential Arrhenius probability term, describing the crack velocity, ranges over more than three orders of magnitude while (\(G_c-G\)) ranges over less than two decades.
Continuing with the comparison of our model with pre-existing ones, we had, in our case, to calibrate the disorder to the experimental data, in particular to accurately reproduce the \(\eta\) exponent, that is, to reproduce the fat tail of the crack velocity distribution. Paradoxically, this exponent, which is not accounted for by the other line model, has been found to be rather constant across different experiments and experimental set-ups. It could indicate that, in practice, the disorder obtained experimentally from the blasting and sintering of PMMA plates has always been relatively similar. Such qualitative statement is of course difficult to verify, because there exists no direct way of measuring the fracture energy of the experimental sintered samples. From Fig. 5b, one can yet notice that the calibration of the disorder amplitude does not need to be particularly accurate to obtain a qualitative fit to the experimental velocity distribution. The spread of the \(\eta\) exponent, for large disorders, is not that important in our model for the range of considered \(\delta G_c\), which can also be seen in Fig. 5 of Cochard et al.8. Gjerden et al.57 suggested that the nucleation of damages, predicted by their fiber bundle model, led to a new - percolation - universality class for the propagation of cracks, explaining in particular the robustness of the exponent \(\eta\). Their studies are however also numerical and cover a finite range of disorders, and an extra analytical proof would be needed to show that a system of infinite size would lead exactly to the same exponent, for any disorder distribution shape and amplitude.
Despite the variety in models reproducing the rough dynamics of creep, the present work provides additional indications that a thermodynamics framework in the sense of a thermally activated subcritical crack growth is well suited for the description of creeping cracks. Such a framework has long been considered (e.g.,32,33,34,36,58,59), and, additionally to the scaling laws that we have here presented, the proposed model was proven to fit many other observable features of the physics of rupture8,37,39,40. It accurately recreates the mean advancement of cracks under various loading conditions8,37, including when a front creeps in a spontaneous (not forced) relaxation regime, which cannot be achieved with the other (non subcritical) models, predicting an immobile front. When coupled with heat dissipation at the fracture tip, our description also accounts for the brittleness of matter40 and for its brittle-ductile transition39.
Indeed, for zero dimensional (scalar) crack fronts, it was shown40 that the thermal fluctuation at the crack tip, expressed as a deviation of the temperature from \(T_0\) in Eq. (1), can explain the transition between creep and abrupt rupture, that is, the transition to a propagation velocity close to a mechanical wave speed \(V_0\), five orders of magnitude higher than the maximal creep velocity V that was here modelled. It was also shown, similarly to many phase transition problems, that such a thermal transition could be favoured by material disorder39. Thus, a direct continuation of the present work could be to introduce such a heat dissipation for interfacial cracks in order to study how brittle avalanches nucleate at given positions (typically positions with weaker \(G_c\)) to then expand laterally to become bulk threatening events.
Lawn, B. Fracture of Brittle Solids. Cambridge Solid State Science Series 2nd edn. (Cambridge University Press, 1993).
Gerard, D. & Koss, D. Porosity and crack initiation during low cycle fatigue. Mater. Sci. Eng. A 129, 77–85. https://doi.org/10.1016/0921-5093(90)90346-5 (1990).
Gao, H. & Rice, J. R. A first-order perturbation analysis of crack trapping by arrays of obstacles. J. Appl. Mech. 56, 828–836. https://doi.org/10.1115/1.3176178 (1989).
Article ADS MATH Google Scholar
Bonamy, D., Santucci, S. & Ponson, L. Crackling dynamics in material failure as the signature of a self-organized dynamic phase transition. Phys. Rev. Lett. 101, 045501. https://doi.org/10.1103/PhysRevLett.101.045501 (2008).
Laurson, L., Santucci, S. & Zapperi, S. Avalanches and clusters in planar crack front propagation. Phys. Rev. E 81, 046116. https://doi.org/10.1103/PhysRevE.81.046116 (2010).
Article ADS CAS Google Scholar
Schmittbuhl, J., Hansen, A. & Batrouni, G. G. Roughness of interfacial crack fronts: Stress-weighted percolation in the damage zone. Phys. Rev. Lett. 90, 045505. https://doi.org/10.1103/PhysRevLett.90.045505 (2003).
Danku, Z., Kun, F. & Herrmann, H. J. Fractal frontiers of bursts and cracks in a fiber bundle model of creep rupture. Phys. Rev. E 92, 062402. https://doi.org/10.1103/PhysRevE.92.062402 (2015).
Cochard, A., Lengliné, O., Måløy, K. J. & Toussaint, R. Thermally activated crack fronts propagating in pinning disorder: Simultaneous brittle/creep behavior depending on scale. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci.https://doi.org/10.1098/rsta.2017.0399 (2018).
Article MATH Google Scholar
Wiese, K. J. Theory and experiments for disordered elastic manifolds, depinning, avalanches, and sandpiles (2021). Preprint, arXiv:2102.01215.
Griffith, A. The phenomena of rupture and flow in solids. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 221, 163–198. https://doi.org/10.1098/rsta.1921.0006 (1921).
Irwin, G. R. Analysis of stresses and strains near the end of a crack traversing a plate. J. Appl. Mech. 24, 361–364 (1957).
Barés, J., Hattali, M. L., Dalmas, D. & Bonamy, D. Fluctuations of global energy release and crackling in nominally brittle heterogeneous fracture. Phys. Rev. Lett. 113, 264301. https://doi.org/10.1103/physrevlett.113.264301 (2014).
Turquet, A. L. et al. Source localization of microseismic emissions during pneumatic fracturing. Geophys. Res. Lett. 46, 3726–3733. https://doi.org/10.1029/2019GL082198 (2019).
Vu, C.-C. & Weiss, J. Asymmetric damage avalanche shape in quasibrittle materials and subavalanche (aftershock) clusters. Phys. Rev. Lett. 125, 105502. https://doi.org/10.1103/PhysRevLett.125.105502 (2020).
Santucci, S., Planet, R., Måløy, K. J. & Ortín, J. Avalanches of imbibition fronts: Towards critical pinning. Europhys. Lett. 94, 46005. https://doi.org/10.1209/0295-5075/94/46005 (2011).
Vives, E. et al. Distributions of avalanches in martensitic transformations. Phys. Rev. Lett. 72, 1694–1697. https://doi.org/10.1103/PhysRevLett.72.1694 (1994).
Dimiduk, D. M., Woodward, C., LeSar, R. & Uchic, M. D. Scale-free intermittent flow in crystal plasticity. Science 312, 1188–1190. https://doi.org/10.1126/science.1123889 (2006).
Durin, G. & Zapperi, S. Scaling exponents for Barkhausen avalanches in polycrystalline and amorphous ferromagnets. Phys. Rev. Lett. 84, 4705–4708. https://doi.org/10.1103/PhysRevLett.84.4705 (2000).
Sethna, J. P., Dahmen, K. A. & Myers, C. R. Crackling noise. Nature 12, 12. https://doi.org/10.1038/35065675 (2001).
Laurson, L. et al. Evolution of the average avalanche shape with the universality class. Nat. Commun. 4, 242–250. https://doi.org/10.1038/ncomms3927 (2013).
Jolivet, R. et al. The burst-like behavior of aseismic slip on a rough fault: The creeping section of the Haiyuan Fault, China. Bull. Seismol. Soc. Am. 105, 480–488. https://doi.org/10.1785/0120140237 (2014).
Rousset, B. et al. An aseismic slip transient on the north anatolian fault. Geophys. Res. Lett. 43, 3254–3262. https://doi.org/10.1002/2016GL068250 (2016).
Grob, M. et al. Quake catalogs from an optical monitoring of an interfacial crack propagation. Pure Appl. Geophys. 166, 777–799. https://doi.org/10.1007/s00024-004-0496-z (2009).
Lengliné, O. et al. Downscaling of fracture energy during brittle creep experiments. J. Geophys. Res. Solid Earthhttps://doi.org/10.1029/2010JB008059 (2011).
Lengliné, O. et al. Interplay of seismic and aseismic deformations during earthquake swarms: An experimental approach. Earth Planet. Sci. Lett. 331–332, 215–223. https://doi.org/10.1016/j.epsl.2012.03.022 (2012).
Måløy, K. J., Santucci, S., Schmittbuhl, J. & Toussaint, R. Local waiting time fluctuations along a randomly pinned crack front. Phys. Rev. Lett. 96, 045501. https://doi.org/10.1103/PhysRevLett.96.045501 (2006).
Santucci, S. et al. Fracture roughness scaling: A case study on planar cracks. Europhys. Lett. 92, 44001. https://doi.org/10.1209/0295-5075/92/44001 (2010).
Tallakstad, K. T., Toussaint, R., Santucci, S., Schmittbuhl, J. & Måløy, K. J. Local dynamics of a randomly pinned crack front during creep and forced propagation: An experimental study. Phys. Rev. E 83, 046108. https://doi.org/10.1103/PhysRevE.83.046108 (2011).
Tanguy, A., Gounelle, M. & Roux, S. From individual to collective pinning: Effect of long-range elastic interactions. Phys. Rev. E 58, 1577–1590. https://doi.org/10.1103/PhysRevE.58.1577 (1998).
Gjerden, K. S., Stormo, A. & Hansen, A. Local dynamics of a randomly pinned crack front: A numerical study. Front. Phys. 2, 66. https://doi.org/10.3389/fphy.2014.00066 (2014).
Stormo, A., Lengliné, O., Schmittbuhl, J. & Hansen, A. Soft-clamp fiber bundle model and interfacial crack propagation: Comparison using a non-linear imposed displacement. Front. Phys. 4, 18. https://doi.org/10.3389/fphy.2016.00018 (2016).
Brenner, S. S. Mechanical behavior of sapphire whiskers at elevated temperatures. J. Appl. Phys. 33, 33–39. https://doi.org/10.1063/1.1728523 (1962).
Zhurkov, S. N. Kinetic concept of the strength of solids. Int. J. Fract. 26, 295–307. https://doi.org/10.1007/BF00962961 (1984).
Santucci, S., Vanel, L. & Ciliberto, S. Subcritical statistics in rupture of fibrous materials: Experiments and model. Phys. Rev. Lett. 93, 095505. https://doi.org/10.1103/PhysRevLett.93.095505 (2004).
Santucci, S., Cortet, P.-P., Deschanel, S., Vanel, L. & Ciliberto, S. Subcritical crack growth in fibrous materials. Europhys. Lett. 74, 595–601. https://doi.org/10.1209/epl/i2005-10575-2 (2006).
Vanel, L., Ciliberto, S., Cortet, P.-P. & Santucci, S. Time-dependent rupture and slow crack growth: Elastic and viscoplastic dynamics. J. Phys. D Appl. Phys. 42, 214007. https://doi.org/10.1088/0022-3727/42/21/214007 (2009).
Lengliné, O. et al. Average crack-front velocity during subcritical fracture propagation in a heterogeneous medium. Phys. Rev. E 84, 036104. https://doi.org/10.1103/PhysRevE.84.036104 (2011).
Tallakstad, K., Toussaint, R., Santucci, S. & Måløy, K. Non-gaussian nature of fracture and the survival of fat-tail exponents. Phys. Rev. Lett. 110, 145501. https://doi.org/10.1103/PhysRevLett.110.145501 (2013).
Vincent-Dospital, T., Toussaint, R., Cochard, A., Måløy, K. J. & Flekkøy, E. G. Thermal weakening of cracks and brittle-ductile transition of matter: A phase model. Phys. Rev. Mater.https://doi.org/10.1103/PhysRevMaterials.4.023604 (2020).
Vincent-Dospital, T. et al. How heat controls fracture: The thermodynamics of creeping and avalanching cracks. Soft Matter 12, 12. https://doi.org/10.1039/d0sm010 (2020).
Hammes, G. G. Principles of Chemical Kinetics (Academic Press, 1978).
Freund, L. B. Crack propagation in an elastic solid subjected to general loading. J. Mech. Phys. Solids 20, 129–152. https://doi.org/10.1016/0022-5096(72)90006-3 (1972).
Santucci, S., Måløy, K. J., Toussaint, R. & Schmittbuhl, J. Self-affine scaling during interfacial crack front propagation. NATO, ASI, Geilo. In Dynamics of Complex Interconnected Systems (eds Skjeltorp, A. T. & Belushkin, A. V.) (Springer, 2006).
Hattali, M., Barés, J., Ponson, L. & Bonamy, D. Low velocity surface fracture patterns in brittle material: A newly evidenced mechanical instability. In THERMEC 2011, Materials Science Forum Vol. 706, 920–924 (Trans Tech Publications Ltd, 2012) https://doi.org/10.4028/www.scientific.net/MSF.706-709.920.
Perrin, G., Rice, J. R. & Zheng, G. Self-healing slip pulse on a frictional surface. J. Mech. Phys. Solids 43, 1461–1495. https://doi.org/10.1016/0022-5096(95)00036-I (1995).
Article ADS MathSciNet MATH Google Scholar
Dormand, J. R. & Prince, P. J. A family of embedded Runge–Kutta formulae. J. Comput. Appl. Math. 6, 19–26. https://doi.org/10.1016/0771-050X(80)90013-3 (1980).
Article MathSciNet MATH Google Scholar
Hairer, E., Nørsett, S. P. & Wanner, G. Solving Ordinary Differential Equations I, Nonstiff Problems (Springer, 1993).
MATH Google Scholar
Cochard, A., Lengliné, O., Vincent-Dospital, T., Toussaint, R., Santucci, S. & Måløy, K. J. CreepyFractures, Software Heritage, swh:1:dir:03ab8789332dc44bb7699b36c1bb6abfc061dbb4 (2021).
Zerwer, A., Polak, M. A. & Santamarina, J. C. Wave propagation in thin plexiglas plates: Implications for Rayleigh waves. NDT E Int. 33, 33–41. https://doi.org/10.1016/S0963-8695(99)00010-9 (2000).
Barabási, A.-L. & Stanley, H. E. Fractal Concepts in Surface Growth (Cambridge University Press, 1995).
Book Google Scholar
Soriano, J. et al. Anomalous roughening of viscous fluid fronts in spontaneous imbibition. Phys. Rev. Lett. 95, 104501. https://doi.org/10.1103/PhysRevLett.95.104501 (2005).
Måløy, K. J., Toussaint, R. & Schmittbuhl, J. Dynamics and structure of interfacial crack front. In 11th International Conference on Fracture 2005, ICF11 7 (2005).
Technical information, Altuglas sheets. Tech. Rep., Arkema (2017).
Marshall, G. P., Coutts, L. H. & Williams, J. G. Temperature effects in the fracture of PMMA. J. Mater. Sci. 9, 1409–1419. https://doi.org/10.1007/BF00552926 (1974).
Santucci, S. et al. Avalanches and extreme value statistics in interfacial crackling dynamics. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 377, 20170394. https://doi.org/10.1098/rsta.2017.0394 (2019).
Bouchaud, E., Bouchaud, J., Fisher, D., Ramanathan, S. & Rice, J. Can crack front waves explain the roughness of cracks?. J. Mech. Phys. Solids 50, 1703–1725. https://doi.org/10.1016/S0022-5096(01)00137-5 (2002).
Gjerden, K. S., Stormo, A. & Hansen, A. Universality classes in constrained crack growth. Phys. Rev. Lett. 111, 135502. https://doi.org/10.1103/PhysRevLett.111.135502 (2013).
Scorretti, R., Ciliberto, S. & Guarino, A. Disorder enhances the effects of thermal noise in the fiber bundle model. Europhys. Lett. 55, 626–632. https://doi.org/10.1209/epl/i2001-00462-x (2001).
Roux, S. Thermally activated breakdown in the fiber-bundle model. Phys. Rev. E 62, 6164–6169. https://doi.org/10.1103/PhysRevE.62.6164 (2000).
Toussaint, R. et al. How cracks are hot and cool: A burning issue for paper. Soft Matter 12, 5563–5571. https://doi.org/10.1039/C6SM00615A (2016).
Vincent-Dospital, T., Toussaint, R., Cochard, A., Flekkøy, E. G. & Måløy, K. J. Thermal dissipation as both the strength and weakness of matter. A material failure prediction by monitoring creep. Soft Matter 12, 12. https://doi.org/10.1039/D0SM02089C (2021).
Dugdale, D. Yielding of steel sheets containing slits. J. Mech. Phys. Solids 8, 100–104. https://doi.org/10.1016/0022-5096(60)90013-2 (1960).
Henkee, C. S. & Kramer, E. J. Crazing and shear deformation in crosslinked polystyrene. J. Polym. Sci. Polym. Phys. Ed. 22, 721–737. https://doi.org/10.1002/pol.1984.180220414 (1984).
Anderson, T. L. Fracture Mechanics: Fundamentals and Applications (Taylor and Francis, 2005).
The authors declare no competing interests in the publishing of this work. They acknowledge the support of the Universities of Strasbourg and Oslo, of the CNRS INSU ALEAS program and of the IRP France-Norway D-FFRACT. We thank the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644. We are also grateful for the support of the Lavrentyev Institute of Hydrodynamics, through Grant No 14.W03.31.0002 of the Russian Government.
ITES UMR 7063, Université de Strasbourg, 67084, Strasbourg, France
Tom Vincent-Dospital, Alain Cochard & Renaud Toussaint
SFF Porelab, The Njord Centre, Department of physics, University of Oslo, Oslo, Norway
Tom Vincent-Dospital, Knut Jørgen Måløy & Renaud Toussaint
ENS de Lyon, Université Claude Bernard, CNRS, Laboratoire de Physique, Université de Lyon, Lyon, France
Stéphane Santucci
Lavrentyev Institute of Hydrodynamics, Siberian Branch of the Russian Academy of Sciences, Novosibirsk, Russia
Tom Vincent-Dospital
Alain Cochard
Knut Jørgen Måløy
Renaud Toussaint
T.V.-D. developed and analyzed the model and the simulations, and redacted the first versions of the manuscript. A.C. set the basis for the numerical implementation of the model and the principles of the resolution algorithm. R.T. proposed the physical basis of the model and its mathematical formulation. K.J.M. and S.S. contributed in the interpretation of the model in fracture mechanics applications, and advised on the the experimental results. All authors participated to the redaction of the manuscript and agreed with the submitted version.
Correspondence to Tom Vincent-Dospital, Alain Cochard or Renaud Toussaint.
The authors declare no competing interests.
Supplementary Information.
Vincent-Dospital, T., Cochard, A., Santucci, S. et al. Thermally activated intermittent dynamics of creeping crack fronts along disordered interfaces. Sci Rep 11, 20418 (2021). https://doi.org/10.1038/s41598-021-98556-x
DOI: https://doi.org/10.1038/s41598-021-98556-x