# Soft-Clamp Fiber Bundle Model and Interfacial Crack Propagation: Comparison Using a Non-linear Imposed Displacement

Arne Stormo, Olivier Lengliné, Jean Schmittbuhl, Alex Hansen

► To cite this version:

Arne Stormo, Olivier Lengliné, Jean Schmittbuhl, Alex Hansen. Soft-Clamp Fiber Bundle Model and Interfacial Crack Propagation: Comparison Using a Non-linear Imposed Displacement. Frontiers in Physics, 2016, 4, 10.3389/fphy.2016.00018. [https://doi.org/10.3389/fphy.2016.00018](https://doi.org/10.3389/fphy.2016.00018). [arXiv:2012.00018](https://arxiv.org/abs/2012.00018)

HAL Id: hal-03501003

https://hal.science/hal-03501003

Submitted on 23 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Soft-Clamp Fiber Bundle Model and Interfacial Crack Propagation: Comparison Using a Non-linear Imposed Displacement

Arne Stormo¹, Olivier Lengliné¹, Jean Schmittbuhl¹* and Alex Hansen²

¹École et Observatoire des Sciences de la Terre, Institut de Physique du Globe de Strasbourg, Centre National de la Recherche Scientifique, Université de Strasbourg, Strasbourg, France, ²Department of Physics, Norwegian University of Science and Technology, Trondheim, Norway

We compare experimental observations of a slow interfacial crack propagation along an heterogeneous interface to numerical simulations using a soft-clamped fiber bundle model. The model consists of a planar set of brittle fibers between a deformable elastic half-space and a rigid plate with a square root shape that imposes a non-linear displacement around the process zone. The non-linear square-root rigid shape combined with the long range elastic interactions is shown to provide more realistic displacement and stress fields around the crack tip in the process zone and thereby significantly improving the predictions of the model. Experiments and model are shown to share a similar self-affine roughening of the crack front both at small and large scales and a similar distribution of the local crack front velocity. Numerical predictions of the Family-Viscek scaling for both regimes are discussed together with the local velocity distribution of the fracture front.

**Keywords:** interfacial crack, asperity pinning, self-affine roughness, local crack velocity, Family-Vicsek scaling, fiber bundle model

OPEN ACCESS

**Edited by:**
Daniel Bonamy,
Commissariat à l'énergie atomique et
aux énergies alternatives, France

**Reviewed by:**
Purusattam Ray,
The Institute of Mathematical
Sciences, India
Jonathan Bares,
Duke University, USA
Laurent Ponson,
Centre National de la Recherche
Scientifique, France

*Correspondence:
Jean Schmittbuhl
jean.schmittbuhl@unistra.fr

**Specialty section:**
This article was submitted to
Interdisciplinary Physics,
a section of the journal
Frontiers in Physics

Received: 04 November 2015

Accepted: 19 April 2016

Published: 11 May 2016

**Citation:**

Stormo A, Lengliné O, Schmittbuhl J
and Hansen A (2016) Soft-Clamp
Fiber Bundle Model and Interfacial
Crack Propagation: Comparison
Using a Non-linear Imposed
Displacement. Front. Phys. 4:18.
doi: 10.3389/fphy.2016.00018

## 1. INTRODUCTION

Crack propagation in heterogeneous media is a rich problem which involves the interplay of various physical processes. The problem has been intensively investigated theoretically, numerically, and experimentally, but a unifying model capturing all the experimental features has not been entirely achieved [1, 2] despite its broad range of implications in engineering and Earth sciences problems [3, 4].

During most regular fracture experiments, the fracture surface can only be observed post mortem [5]. Since Mandelbrot's discovery [6] of fracture surfaces in metal that display fractal properties, there have been many attempts to understand the scaling relations of fracture [5]. In particular, the self-affine nature of fracture surfaces have come under much scrutiny. A self-affine surface $h(x)$ has the following scaling relation [7]

$$h(x) \propto \lambda^{\zeta} h(\lambda x), \quad (1)$$

where $\lambda$ is an arbitrary scaling factor and $\zeta$ is the roughness exponent.

The slow propagation of a crack front where long range elastic interactions are dominant, is of crucial importance to fill the gap between experiments and models. In an attempt to study the in-situ properties of this process, Schmittbuhl and Måløy [8–10] developed an experiment where
they forced the fracture to advance along a weak heterogeneous plane embedded in a transparent medium (Plexiglas). In such an experimental configuration, the propagation of the fracture front can be monitored in great detail using high resolution optical devices from which different properties characterizing the fracture dynamics can be inferred. For instance, the distribution of the local fracture velocities [11] or the multi-scale morphology of the fracture front [12] can notably be measured and are found to be robust parameters when tested against several experimental conditions.

Several theoretical and numerical studies have been devoted to quasi-static models [1, 13]. Such models give rise to an intermittent local activity characterized by a depinning transition and can be viewed as a critical phenomenon [2, 14, 15]. However, these models do not provide a complete modeling of the system as they fail to reproduce all the experimental observations. Notably the front morphology in these models does not display a cross-over length scale with two different roughness exponents above and below it. Such a cross-over was observed experimentally in 2010 by Santucci et al. [12]. They found that the in-plane fractures also had a scale-dependent roughness exponent. At small scales they obtained $\zeta^{-} = 0.60$, while at large scales they showed that $\zeta^{+} = 0.35$. Attempts to model the behavior of the in-plane fracture process had a long time resulted in only the large scale roughness exponent [16, 17], or the small scale roughness exponent in the context of a stress-weighted percolation [18]. Moreover, the scaling relations for fractures have now been established as being both scale and direction dependent [19–22].

Recently, Gjerden et al. [23, 24] presented for the first time roughness exponents compatible with both the large scale and small scale values found in the experiments. Their model was a development of the soft clamp fiber bundle model [25, 26] first introduced by Batrouni et al. [27]. Gjerden et al. implemented this model with a linear spatial gradient in the threshold distribution thereby creating a crack front whose dynamics and morphology could be compared to the experimental observations [24].

Even though the characteristics of the crack front that was created in the soft clamp fiber bundle model through the introduction of a linear gradient matched the experimental observations well, it cannot be ignored that the crack front was generated in a very different way in the experiment. It is therefore important to implement the same boundary conditions in the soft clamp model as was used in the experiments to test whether this change does lead to different results. This is the primary aim of this paper.

We also push further the comparison between theory and experiments, by focusing essentially on (i) the roughness exponent $\zeta$ that characterizes the geometry of the front at some given time (ii) the dynamic exponent $\kappa$ that characterizes the time evolution of the crack front roughness and (iii) the statistics of the velocity fluctuations that show a scaling behavior with an exponent $\eta$ that is still poorly understood. We note in particular that $\kappa$ was not measured directly by Gjerden et al. [24], but rather inferred from assuming that the Family-Vicsek scaling assumption holds.

In the following, we first review in Section 2 the main aspects of the experimental setup. In Section 3, we give details on the numerical model based on the soft clamp fiber bundle approach but modified to include both long range elastic interactions owing to the elastically deformable clamped block and a non-linear imposed displacement around the process zone. In Section 4, we compare the experimental observations and numerical results. The predictions of the transient and steady-state scaling of the front morphology is discussed in the light of the experimental evidences. Finally the comparison of the local crack front velocity distribution is presented.

## 2. EXPERIMENTAL SETUP

We compare our numerical results to experimental results obtained with the following setup. It consists of two transparent welded plexiglas (PMMA - poly methyl methacrylate) plates [e.g., 8, 11, 28, 29]. Disorder in the strength of the interface between the two plates is introduced by sandblasting using glass beads of variable diameters (180–300 µm), the opposite surfaces of the two plates that are to be welded together. Welding of the two plates is then achieved by imposing a normal load on the assembled plates while heating the system at 190°C. This thermal annealing produces a cohesive interface weaker than the bulk, along which the sample will break apart under normal opening loading. When the system is loaded from one side, a propagating crack—a fracture front—evolves that moves away from the loaded side. It is this crack we focus upon.

The dimensions are typically 20 cm long, 3 cm wide and 5 mm thick for the bottom plexiglas plate and 1 cm thick for the upper plate. The upper plate is attached to a stiff aluminum frame while a load is applied over the top side of the bottom plate in a direction normal to the plate interface. The vertical displacement imposed on the bottom plate induces stable mode I propagation of a planar fracture along the prescribed weak interface (mode II and III components related to the finite size of the sample and friction along the loading rod are neglected). A sketch of the experiment is given in **Figure 1**.

The fracture front creates an optical contrast between the broken and unbroken sections of the interface [28]. This fracture front is then tracked by an optical digital camera (a Nikon D700 with a resolution of 4256 × 2832 pixels at up to 5 images/s for slow monotoring or an Optronis Camrecord 600 with a resolution of 1000 × 1000 pixels at 1000 images/s for fast tracking) in order to extract its position through time [29–31]. The front propagates along the *y* axis where the origin is defined at the loading point and is positive in the direction of crack propagation (see **Figures 1, 2**). The *x* axis is perpendicular to *y* and defines the coordinate of a point along the front. The large scale bending of the plate at the sample size is well reproduced by the elastic beam theory [29]. At a smaller scale i.e., at the scale of the process zone, the elastic beam theory might however not be any more valid in describing the shape of the plate and typical crack solutions eventually provide more a realistic description as discussed below.
<figure><img src="image_2.png" /><figcaption>**FIGURE 1 | Sketch of the experimental setup.** Two PMMA plates of thickness $h_1$ and $h_2$ are welded and a mode I fracture propagates at the interface between the two plates. Fracture propagation is caused by the opening, $M$, between the two plates. The distance between the loading point and the crack front is noted $L_a$. The crack front propagates along the y direction. The gray rectangle marks the area around the crack front (see **Figure 2**) we will reproduce in our simulations.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**FIGURE 2 | Picture of the crack front area obtained during an experiment.** The crack is propagating from top to bottom, the gray part refers to already broken part and the dark area to the unbroken part. The optical contrast between this two parts defines the fracture front. The dimension of the picture in the x-direction is 4 mm.</figcaption></figure>

## 3. MODEL

The soft clamp fiber bundle model was first introduced by Batrouni et al. [27] in 2002. We use this model in order to reproduce the phenomena occurring at the process zone scale around the fracture front by introducing the appropriate boundary conditions. We are indeed interested in describing the behavior of the fracture front on a short range from its edge. We thus consider our model as a square of dimension $L_f \times L_f$ around the fracture front. This area of the system is illustrated in **Figure 1** by the dark square covering a part of the fracture front. This area is divided into $N \times N$ sub-squares, each of them containing a single fiber. A discretization size is then $dl = L_f/N$. The two facing Plexiglas plates are represented as infinite facing half-spaces, and act as clamps attached to the $N^2$ fibers. The half spaces are supposed to have an elastic behavior with a Young's modulus $E$ and a Poisson ratio of $\nu$. An equivalent simpler configuration can be obtained by assuming that one of the half spaces is infinitely stiff. Thus, $E$ will only describe the elastic response of the other half-space. The system is periodic along the x-axis.

The fibers are linearly elastic up to a strength threshold where they break irreversibly and can no longer support any load. On the contrary to Gjerden et al. [24], the thresholds are uniformly spatially distributed between $f_{min}^t$ and $f_{max}^t$, where we have set $f_{min}^t = 0$ and $f_{max}^t = 1$. The force, $f_i$, experienced by each individual fiber is given by the Hooke's law:

$$f_i = -k(u_i - D_i), \quad (2)$$

where $i$ is the index of the fiber, $k$ is the spring constant, equal for all the fibers, $D_i$ is the local position of the rigid plate (e.g., constant for all fibers if the plate is plane and horizontal), $u_i$ is the local displacement of the facing soft elastic half-space created by the strain around the fiber by of this elastic half-space. When $D_i = u_i$, i.e., the imposed displacement is balancing the strain in the elastic half-space, the force in the fiber is nil. The spatially dependent deformation is calculated by Love's law [32]:

$$u_i = \sum_{j \neq i}^{N^2} G_j f_j \quad (3)$$

where $j$ runs over all the fibers in the array (except the self-induced contribution). The Green function $G_j$ is given by:

$$G_j = \frac{1 - \nu^2}{\pi E d l^2} \int_{-dl/2}^{+dl/2} \int_{-dl/2}^{+dl/2} \frac{dx' dy'}{|\vec{r}_i(0, 0) - \vec{r}_j(x', y')|}, \quad (4)$$

Integrals in Equation (4) run over each $dl^2$ square associated with the fiber $j$. The $x'$ and $y'$ coordinates have their origin in the middle of the discretization area. The vector $\vec{r}_i(x', y')$ represents the position of the fiber $i$. $E$ and $\nu$ are respectively the Young and Poisson coefficients.

During loading in the experiments, the large scale mode I opening, induced by the imposed displacement of the bottom plate, is well reproduced by an elastic beam model [30]. On the contrary, at the process zone scale around the crack tip, the opening is strongly influenced by the fracture mechanics [3]. The opening of the crack $\delta(y)$ is then expected to have a square root shape behind the crack front [33]:

$$\delta(y_c - y) = \frac{8(1 - \nu^2)}{\sqrt{2\pi}} \left(\frac{K}{E}\right) \sqrt{y_c - y} \quad (5)$$

where $y_c$ is the average position of the crack tip and $K$ is the stress intensity factor.

Since our investigations focus on the small scales around the fracture front, we assume that elastic crack opening is dominant. Accordingly, we introduced a square root shape of the rigid plate (see **Figure 3**). The imposed displacement $D(y)$ to the fiber bundle is then proposed to be proportional to the square root of the distance $y_c - y$ to the fracture front [33].

We tested the approach using a non-linear shape of the rigid plate by considering a system with half of the fibers broken (see insert of **Figure 4**) and computing the resulting opening $\delta(y) = D(y) - u(y)$ in the domain of broken fibers, in the limit of infinite strength of the surviving fibers in the rest of the domain. We
<figure><img src="image_4.png" /><figcaption>**FIGURE 3** | Sketch of the loading of the fiber-bundle model. Individual fibers are represented by dark vertical lines and can be broken or unbroken. The upper plate is the soft elastic block with a Young modulus *E*. The shape of the rigid bottom plate follows the displacement $D \propto \sqrt{L_f - y}$, as given by Equation (5).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**FIGURE 4** | Insert: sketch of the half broken configuration considered in the figure in the case of a flat rigid plate (gray) facing a soft elastic block (white). Fibers in red are only set in the right half of the system. For this test, there are supposed to have an infinite strength, i.e., no possible failure. Main figure: Crack tip opening, $\delta(y)$, as a function of the distance to the first intact fiber, $y_c$. Two cases are considered for the imposed shape of the rigid plate: a flat shape in red and a crack-like shape in blue. A square-root crack fit is displayed as a black line.</figcaption></figure>

compared the configuration of a flat shape $D(y) = D_0$ and a square root shape $D(y) = D_0\sqrt{L_f - y}$ of the rigid plate (see **Figure 4**). The red line shows the opening field $\delta(y)$ between the rigid and the elastic clamps as a function of $y$ for a system of size $N = 256$. The linear crack tip fracture front is at $y_c = 129$.

When we fit the opening $\delta(y)$ with a square root shape similar to Equation (5), we see the difference between both configurations. In the case of a flat plate, the relation only holds in the immediate vicinity of the crack tip. As we increase the

<figure><img src="image_6.png" /><figcaption>**FIGURE 5** | Profile of the force $f(y)$ applied on the fibers ahead of the crack tip for two shapes of the rigid plate (flat in red or square root in blue) in a configuration similar to Figure 4.</figcaption></figure>

distance back to the front, $\delta(y)$ soon deviates from the theoretical crack tip opening. This deviation mostly results from the periodic property of the flat shape. With a square root shape illustrated by the blue line in **Figure 4**, the resulting crack tip opening is nicely fitted by Equation (5) and holds for a much larger portion of the system, over all the process zone of the front.

The force experienced by the fibers ahead of the crack tip for the different loading scheme is shown in **Figure 5**. We observe, for both the flat and square-root rigid plate shape that the force on the fibers, depends on the distance from the crack tip as an inverse square root profile consistently with elastic crack theory [3]. As we can see from the figure, the imposed square-root shape improves significantly the system compliance except at very large distance ahead of the crack front preventing failure of the fibers in this domain.

We see from both **Figures 4, 5** that the loading and the response of the model become non-physical close to the boundaries of the system. In order to keep the process zone around the fracture tip away from the boundaries, we used the “conveyor belt” technique [34]. Since only the surviving fibers matter in the force field calculations and fibers fail irreversibly, we remove the first broken (i.e., left side in **Figure 3**) row from our calculations and add a fully intact one at the last (i.e., right side) row. $n_c$ is the number of times a new intact row is introduced during the computation.

We prevented the last $(1/2)N \log(N)$ rows of fibers from breaking in order to prevent wrap-around of the failed fibers due to the periodic boundary conditions. For $L = 256$, this meant that four rows were unbreakable.

During the computation of the crack front advance, the system is updated iteratively and quasistatically by removing the fibers one by one following an event-driven algorithm [16]. At each step we compute a new prefactor $D_0$ that will correspond to a force field where only one of the fibers is experiencing a force equal to its strength threshold. This fiber is actually the closest one to failure and thus requires the least increase of $D_0$ to break. The
fiber is then removed and the procedure is continued until there is no intact fiber left.

When the distances are measured in units of the discretization size $dl$, the relative influence of the elastic interactions with respect to the global loading has been found to be described by a single parameter named the reduced Young modulus [35]

$$e = \frac{Edl}{N}. \qquad (6)$$

where $E$ is the Young modulus, $N$ the number of fibers per system length $L_f$. Indeed, the contribution of the Green function in the force computation of each fiber can be turned on or off, depending on the value $e$. In the stiff regime, i.e., large $e$, the fracture process is dominated by a diffusive damage controlled by the strength heterogeneities in the system. When $e$ decreases, the system goes through a cross-over regime where the influence of the elastic forces rises. When $e$ is sufficiently low, the competition between the elastic redistribution and the strength heterogeneities are balanced and the system enters a so-called soft regime. This regime is dominated by localized damage. As we can see from Equation (6), decreasing $E$, while increasing $N$ in a way that keeps $e$ constant is not affecting the influence of the Green function. We will thus associate the soft systems with large length scales. Similarly, we will associate the stiff systems with short length scales.

<figure><img src="image_7.png" /><figcaption>FIGURE 6 | Example of process zone during a simulation for a stiff system (e = 62.5). The dark points refer to broken fibers while yellow points are intact fibers. The system size is N = 32.</figcaption></figure>

# 4. RESULTS

## 4.1. Front Morphology: Self-Affine Scaling and Cross-Over Length Scale

Previous experimental studies reported that an interfacial crack front in the present configuration is self-affine with a roughness exponent $\zeta \simeq 0.6$ [e.g., 8, 28]. More recent data extracted from numerous experiments and at various scales show that actually two distinct regimes emerge depending on the scale of investigation: at small scales the scaling regime is characterized by a roughness exponent $\zeta^{-} \simeq 0.60$ while at large scale the exponent is lower and is found around $\zeta^{+} \simeq 0.35$ [12]. We aim at a direct comparison of the scaling properties of the simulated crack front and the experimental observations. For this, we identify the fracture front as the interface between the cluster of broken fibers spacially connected (i.e., in the sense of a link between first neighbors) to one side of the system and the cluster of surviving fibers connected to the other side. When this interface is identified, it is approximated as a function of the x-coordinate for each time step $t$: $h(x, t)$ by using a solid-on-solid (SOS) algorithm. The front line is indeed deduced by taking the first y-value hit when searching along the y-axis from above or below. An example of this front is given by the red line in Figure 6 for stiff systems, and in Figure 7 for soft systems.

We investigated the scaling properties of the crack front morphology in our simulations in a similar way to the experiments. In order to find the roughness exponent (cf Equation 1) of the fracture fronts produced by our model, we analyzed them using the average wavelet coefficient method

which is shown to be very efficient for data over a limited range of scales [36]. The average wavelet coefficient $W(a)$ is defined as:

$$W(a) = \langle W(a, b) \rangle_b \qquad (7)$$

<figure><img src="image_8.png" /><figcaption>FIGURE 7 | Same as Figure 6 for a soft system (e = 7.6 · 10⁻⁶).</figcaption></figure>
where a is the scale factor of the wavelet transform and b is the translation factor. It scales for a self-affine front as

$$W(a) \propto a^{\zeta+1/2}. \qquad (8)$$

In Figure 8, we see the average wavelet coefficient (using Daubechies wavelets of order 4) for the stiff regime for different system sizes. We find the roughness exponent of the small systems by making a fit over the three smallest coefficients. We exclude the very smallest, as the discretization size begins to influence the results. This gives us a roughness exponent of $\zeta_-$ = 0.53 ± 0.02.

It is clear from Figure 8 that the estimated roughness exponent is size dependent. Nevertheless, as the system size increased we observe a convergence of the roughness exponent toward the experimental value represented by the dashed line.

We also performed simulations in the soft regime, which is equivalent to simulating the system on large length scales. We extracted front positions and applied the average wavelet coefficient method to the front morphology. We obtained the roughness exponent by fitting the four largest coefficients for $L = 128$. This gives: $\zeta_+ = 0.35 \pm 0.01$. In Figure 9, we see the wavelet coefficient scaling for this long range regime.

As opposed to the short scale regime, the estimated value of $\zeta_+$ is nearly insensitive to the system size. We also observe that the experimental roughness exponent observed at large length scale [12] is similar to the one we obtained from our simulations.

## 4.2. Family-Vicsek Scaling

As the fracture front progresses through the system, the memory of the earlier positions fades progressively. We explore this process by creating a differential front $h_m$ with respect to the fracture front at time $t_0$:

$$h_m(x, t) = h(x, t) - h(x, t_0). \qquad (9)$$

<figure><img src="image_9.png" /></figure>

The front width is defined as the root mean square (rms):

$$\sigma(L, t) = \sqrt{\sum_{x=0}^{N} (h_m(x, t, t_0) - \langle h_m(x, t, t_0) \rangle)_x^2}. \qquad (10)$$

We chose to use the Family-Vicsek (FV) scaling [37] framework to describe the transient evolution of the differential front:

$$\sigma(t, L) \propto L^{\zeta} f\left(\frac{t}{L^{\kappa}}\right) \qquad (11)$$

where

$$f(x) = \begin{cases} x^{\zeta/\kappa} & x \to 0, \\ \text{constant} & x \to \infty. \end{cases} \qquad (12)$$

and $\kappa$ is the dynamic exponent.

Since there is no inherent time in the model equations, we have to introduce a posteriori time definition in order to be able to compare our results to the properties deduced from the experimental setup. We propose to use a simple definition: the average position of the front $h_m$ as a proxy of the time $t$. Using this time definition, Figures 10 and 11 show the best FV fit for both stiff and soft regimes respectively. It corresponds to the evolution of the rms $\sigma(L, \langle h_m \rangle)$ of the front of a system of width $L$ after an average propagation $\langle h_m \rangle$. Curves for different system sizes $L$ have been superimposed at best, adjusting the roughness exponent $\zeta$, the dynamic exponent $\kappa$ and fulfilling the related slope condition: $\beta = \zeta/\kappa$.

In both stiff and soft systems, we find the roughness exponent $\zeta$ in very good agreement with the results from the average wavelet coefficient method (see Figures 8, 9). The best collapses of the roughness growth curves $\sigma(L, \langle h_m \rangle)$ give: $\kappa = 0.85$ and $\beta = 0.71$ for the stiff regime and $\kappa = 0.6$ and $\beta = 0.7$ for the soft regime.

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 10 | Family-Vicsek scaling for stiff systems i.e., small scale systems (e = 62.5). $\zeta_{-} = 0.6$, $\kappa_{-} = 0.85$ and $\beta_{-} = 0.71$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11 | Family-Vicsek scaling for soft systems i.e., large scale systems (e = 7.6 · 10⁻⁶). $\zeta_{+} = 0.35$, $\kappa_{+} = 0.6$ and $\beta_{+} = 0.58$.</figcaption></figure>

### 4.3. Distribution of Local Velocity

A direct way for studying the dynamics of the fracture propagation is to measure the distribution of the local velocity along the crack front. To obtain this we use the waiting time matrix technique developed by Måløy et al. [11] on experimental data. The waiting time matrix W gives at each sub-squared (i.e., pixel), the time spent by the crack front at this location. This waiting time is measured in our simulations, by identifying the fracture front at each time step (i.e., each broken fiber) and increasing each element of W covered by the crack front by one. After doing this for the entire fracture propagation, we obtain the final estimate of W. From W, we can create the velocity matrix V. Each element of V is indeed the inverse element of

<figure><img src="image_13.png" /><figcaption>FIGURE 12 | Comparison of the local normalized velocity maps obtained for an experiment at two scales (top): at large scale (top right), and a detail at small scale (top left); and two simulations: in the soft regime (bottom right) and in the stiff regime (bottom left) (bottom). White dots correspond to areas where the local velocity was not measurable, i.e., too fast for the local waiting time method.</figcaption></figure>

W times the discretization length $dl^1$. This creates a map of the local velocities of the fracture front. Examples of velocity maps are given in Figure 12 for both regimes together with a map at two scales from experimental measurements [31]. It provides a direct visual observation of the pinning of the crack front in our simulation. Comparisons with experimental data show a good agreement and confirm that the stiff regime describes the details of the propagation at small scales and the soft regime at large scales.

Experimentally the distribution of the local velocities v have been found to follow the relation:

$$P(v/\langle v \rangle) \propto (v/\langle v \rangle)^{-\eta}, \quad \text{for } v/\langle v \rangle > 1, \quad (13)$$

where $\eta = 2.55 \pm 0.15$ and is found to be a robust feature over many experiments [11, 31, 38].

We see that both regimes (small and large e) produce a power-law decay of the velocity distribution above the average crack velocity consistently with experimental results. The best power law fits for the two regimes, stiff and soft regimes, give $\eta_{-} = 2.9$ and $\eta_{+} = 2.1$ respectively.

### 5. DISCUSSION

To better capture the physical meaning of the controlling parameter e, it is of interest to rewrite e as: $e = dl^2/L_f^2 \cdot L_f \cdot E$ which reads as: $e = L_f \cdot E/N^2$ where $N^2$ is the number of fibers in the system. If one has in mind that $L_f \cdot E$ is the effective stiffness of the elastic medium and $N^2 \cdot 1$ is the effective stiffness of the set of fibers in parallel (when they are all present and considering a unit stiffness of each fiber), one realizes that e is the ratio of the stiffness of the bulk to the stiffness of the interface: $e = L_f \cdot E/(N^2 k_f)$ with $k_f$ the stiffness of each fiber. The controlling parameter is then the dimensionless ratio of the two stiffnesses: that of the bulk to that of the interface. When the bulk is very stiff compared to the interface (large e), everything is dominated by

<footnote n="1">¹When W = 0, we replaced the inverse by 9999.9999 in order to ensure that they are indeed treated as large velocities by the program.</footnote>
the elastic bulk and when the interface is very stiff compared to the bulk, one has a so called soft system.

Using a similar numerical approach, Gjerden et al. [23] obtained advances of the fracture front by introducing a linear gradient in the breaking thresholds with a flat rigid plate. The scaling behavior they found was consistent with a gradient percolation modeling at small scales (large $e$): $\zeta_{-} = 2/3$ and with the contact line model at large scales (small $e$): $\zeta_{+} = 0.39$. In the present problem where there is no gradient in the breaking thresholds but at non-linear displacement loading (a square root profile), we obtained slightly different results. Indeed, the non-linearity may change the scaling exponents in the same way as the introduction of non-linear gradients in the percolation problem [39]. We will pursue this issue in a future publication.

The roughness exponent at small scales with the present model is $\zeta_{-} = 0.60$ both from a direct measurement along the fronts using the average wavelet coefficient method and the Family-Vicsek scaling from finite size scaling. This value is very close to the experimental estimate: $\zeta_{-} = 0.60 \pm 0.05$. At large scales, also both numerical and experimental values match well: $\zeta_{+} = 0.35$. Interestingly a recent study [22] reports on a crossover between two scaling regimes of the fracture roughness but invoking a very different mechanism: the porosity of the medium. Below the pore scales, it is argued that the roughness exponent is isotropic with a value close to $\zeta_{-} = 0.8$ but at large scale, the roughness exponent is anisotropic with a roughness exponent in the range [0.27 – 0.35]. Actually, it would be of interest to see how the porosity change of the medium is affecting the stiffness of the medium or of the interface and therefore the $e$ modulus.

In Gjerden et al. [24], the dynamical exponent is $\kappa_{-} = 1.47$ at small scales and $\kappa_{+} = 0.75$ at large scale. Here we obtained $\kappa_{-} = 0.85$ and $\kappa_{+} = 0.6$ which are rather different. Estimates of $\kappa$ are however not obtained in the same way. In Gjerden et al. [24], they are deduced from a collapse of the front width for different gradients of the strength but not from the finite size scaling. This might explain the discrepancy together with the difference between a linear and a non-linear force perturbation around the fracture tip. $\kappa_{-}$ has been measured in an earlier version of the model [18]: $\kappa_{-} = 0.67$ and $\kappa_{+}$ with a different approach in Duemmer and Krauth [40]: $\kappa_{+} = 0.70$.

Comparisons with experimental values of the dynamical exponent $\kappa$ is rather speculative since there is no experimental measure of $\kappa_{-}$ ($e \gg 1$) and $\kappa_{+}$ ($e \ll 1$) but of an effective $\kappa$ for an intermediate $e$ modulus ($e \approx 1$). An example of an experimental measure is $\kappa = 1.2$ reported by Måløy and Schmittbuhl [41]. An on-going work is devoted to the analysis of the role of the $e$ modulus on the Family-Vicsek scaling and the possible onset of anomalous scaling [42].

The local velocity distribution has also been measured by Gjerden et al. [24] using a linear gradient in the threshold distribution to create the crack front. They followed the same procedure for computing the waiting time matrix as we have done here, and found that, at large scale i.e., for soft systems, the velocity distribution is well fitted by an equation of the form of Equation (13) with an exponent $\eta_{+} = 2.53$ close to the experimental value.

<figure><img src="image_14.png" /><figcaption>**FIGURE 13** | Local velocity distributions from experimental data and stiff and soft models. The best power law fits of the velocity distributions provide estimates of the power law exponents: $\eta_{-} = 2.9$ (stiff system in red) and $\eta_{+} = 2.1$ (soft system in blue). For the stiff regime, only velocities in the pinning zones are considered.</figcaption></figure>

In **Figure 13** we can see the experimental data for the velocity distribution [31], as well as the distribution for the stiff and soft regimes in our model based on a non-linear loading of the system.

# 6. CONCLUSIONS

By introducing a physical loading into the fiber bundle model using a rigid square root shaped plate facing a fully elastic half-space, we obtain very similar results to the experimental values. The square root shape is shown to mimic the elastic opening behind the fracture tip. The scaling of fracture front is shown to exhibit two self-affine regimes using the average wavelet coefficient method (AWC): at small scale (i.e., stiff regime), the roughness exponent is $\zeta_{-} = 0.60$ and at large scale (i.e., soft regime), $\zeta_{+} = 0.35$ very similarly to the experimental measures reported by Santucci et al. [12] ($\zeta_{-} = 0.60$ and $\zeta_{+} = 0.35$). The cross-over between both regimes is controlled by the reduced $e$ modulus which is shown to be related to the ratio of the stiffness of the bulk and the stiffness of the interface.

Transient roughening is studied in the Family-Vicsek scaling framework. We measured the evolution of the root mean square of the fracture front roughness from a reference front as a function of the average position of the front. We obtained good collapses of the roughness growth for different system sizes using a roughness exponent $\zeta$ consistent with the direct measurements and a dynamical exponent $\kappa$. For the small scale regime, i.e., the stiff regime, we measured: $\kappa_{-} = 0.85$. For the large scale regime, i.e., the soft regime, we estimated: $\kappa_{+} = 0.60$. Both
values are consistent with the independently measured $\beta = \zeta/\kappa$:
$\beta_{-} = 0.71$ for the stiff regime and $\beta_{+} = 0.58$ for the soft
regime. Our roughness exponents are close to that of Gjerden
et al. [23]. The Gjerden analysis for the hard case was based
on the assumption that the problem was in the universality
class of ordinary percolation. The dynamical exponents $\kappa$ are on
the contrary rather different which is expected to be related to
the differences in loading: linear gradient of the fiber strength
vs. non-linear imposed displacement. The dynamical exponent
seems accordingly more sensitive to the difference in the models.
Comparisons with experiments of the dynamical exponent is
difficult since there is no measure of the dynamical exponents
for the extreme case $e \ll 1$ and $e \gg 1$ but only for intermediate
cases. An interesting forthcoming study will deal with the effect
of the $e$ modulus on the dynamical scaling of the fracture
front.

We also compared the local velocity distribution of our model
to the distribution found in experiments. Both are consistent with
a power law distribution: $P(v/\langle v \rangle) \propto (v/\langle v \rangle)^{-\eta}$ with $\eta_{exp} = 2.55$

for the experiment, $\eta_{-} = 2.9$ for the stiff regime and $\eta_{+} = 2.1$
for the soft regime.

## AUTHOR CONTRIBUTIONS

AS did the programming, produced data, some analysis and a
first draft. OL did experimental work, produced data and some
analysis. He also updated the text. JS did data analysis and
updated the manuscript. AH did data analysis and updated the
manuscript.

## ACKNOWLEDGMENTS

We thank K. S. Gjerden and K. J. Måløy for insightful discussions and A. Steyer for technical support for the experimental work. We furthermore acknowledge support from ANR SUPNAF, Labex G-EAU-THERMIE PROFONDE (grant no. ANR-11-LABX-0050) and the Norwegian Research Council through grant no. 1999770.

## REFERENCES

1. Alava MJ, Nukala PK, Zapperi S. Statistical models of fracture. *Adv Phys.* (2006) **55**:349–476. doi: 10.1080/00018730300741518

2. Bonamy D, Bouchaud E. Failure of heterogeneous materials: a dynamic phase transition? *Phys Rep.* (2011) **498**:1–44. doi: 10.1016/j.physrep.2010.07.006

3. Anderson TL, Anderson T. *Fracture Mechanics: Fundamentals and Applications*. Boca Raton, FL: CRC Press (2005).

4. Scholz CH. *The Mechanics of Earthquakes and Faulting*. Cambridge: Cambridge University Press (2002).

5. Bouchaud E. Scaling properties of cracks. *J Phys.* (1997) **9**:4319. doi: 10.1088/0953-8984/9/21/002

6. Mandelbrot BB, Passoja DE, Paullay AJ. Fractal character of fracture surfaces of metals. *Nature* (1984) **308**:721–2. doi: 10.1038/308721a0

7. Barabási AL. *Fractal Concepts in Surface Growth*. Cambridge: Cambridge University Press (1995).

8. Schmittbuhl J, Måløy KJ. Direct observation of a self-affine crack propagation. *Phys Rev Lett.* (1997) **78**:3888. doi: 10.1103/PhysRevLett.78.3888

9. Delaplace A, Schmittbuhl J, Måløy KJ. High resolution description of a crack front in an heterogeneous plexiglas block. *Phys Rev E* (1999) **60**:1337–43. doi: 10.1103/PhysRevE.60.1337

10. Måløy KJ, Schmittbuhl J. Dynamical events during slow crack propagation. *Phys Rev Lett.* (2001) **87**:105502. doi: 10.1103/PhysRevLett.87.105502

11. Måløy KJ, Santucci S, Schmittbuhl J, Toussaint R. Local waiting time fluctuations along a randomly pinned crack front. *Phys Rev Lett.* (2006) **96**:045501. doi: 10.1103/PhysRevLett.96.045501

12. Santucci S, Grob M, Toussaint R, Schmittbuhl J, Hansen A, Måløy KJ. Fracture roughness scaling: a case study on planar cracks. *Europhys Lett.* (2010) **92**:44001. doi: 10.1209/0295-5075/92/44001

13. Sahimi M. Non-linear and non-local transport processes in heterogeneous media: from long-range correlated percolation to fracture and materials breakdown. *Phys Rep.* (1998) **306**:213–395. doi: 10.1016/S0370-1573(98)00024-6

14. Chen K, Bak P, Obukhov S. Self-organized criticality in a crack-propagation model of earthquakes. *Phys Rev A* (1991) **43**:625. doi: 10.1103/PhysRevA.43.625

15. Dickman R, Muñoz MA, Vespignani A, Zapperi S. Paths to self-organized criticality. *Brazil J Phys.* (2000) **30**:27–41. doi: 10.1590/S0103-97332000000100004

16. Schmittbuhl J, Roux S, Vilotte JP, Måløy KJ. Interfacial crack pinning: effect of nonlocal interactions. *Phys Rev Lett.* (1995) **74**:1787. doi: 10.1103/PhysRevLett.74.1787

17. Rosso A, Krauth W. Roughness at the depinning threshold for a long-range elastic string. *Phys Rev E* (2002) **65**:025101. doi: 10.1103/PhysRevE.65.025101

18. Schmittbuhl J, Hansen A, Batrouni GG. Roughness of interfacial crack fronts: stress-weighted percolation in the damage zone. *Phys. Rev. Lett.* (2003) **90**:045505. doi: 10.1103/PhysRevLett.90.045505

19. Ponson L, Bonamy D, Bouchaud E. Two-dimensional scaling properties of experimental fracture surfaces. *Phys Rev Lett.* (2006) **96**:035506. doi: 10.1103/PhysRevLett.96.035506

20. Ponson L, Auradou H, Vié P, Hulin JP. Low self-affine exponents of fractured glass ceramics surfaces. *Phys Rev Lett.* (2006) **97**:125501. doi: 10.1103/PhysRevLett.97.125501

21. Bonamy D, Ponson L, Prades S, Bouchaud E, Guillot C. Scaling exponents for fracture surfaces in homogeneous glass and glassy ceramics. *Phys Rev Lett.* (2006) **97**:135504. doi: 10.1103/PhysRevLett.97.135504

22. Cambonie T, Bares J, Hattali M, Bonamy D, Lazarus V, Auradou H. Effect of the porosity on the fracture surface roughness of sintered materials: from anisotropic to isotropic self-affine scaling. *Phys Rev E* (2015) **91**:012406. doi: 10.1103/PhysRevE.91.012406

23. Gjerden KS, Stormo A, Hansen A. Universality classes in constrained crack growth. *Phys Rev Lett.* (2013) **111**:135502. doi: 10.1103/PhysRevLett.111.135502

24. Gjerden KS, Stormo A, Hansen A. Local dynamics of a randomly pinned crack front: a numerical study. *Front Phys.* (2014) **2**:66. doi: 10.3389/fphy.2014.00066

25. Pradhan S, Hansen A, Chakrabarti BK. Failure processes in elastci fiber bundles. *Rev Mod Phys.* (2010) **82**:499. doi: 10.1103/RevModPhys.82.499

26. Hansen A, Hemmer PC, Pradhan S. *The Fiber Bundle Model*. Weinheim: Wiley-VCH (2015).

27. Batrouni GG, Hansen A, Schmittbuhl J. Heterogeneous interfacial failure between two elastic blocks. *Phys. Rev E* (2002) **65**:036126. doi: 10.1103/PhysRevE.65.036126

28. Delaplace A, Schmittbuhl J, Måløy KJ. High resolution description of a crack front in a heterogeneous Plexiglas block. *Phys Rev E* (1999) **60**:1337. doi: 10.1103/PhysRevE.60.1337

29. Lengliné O, Toussaint R, Schmittbuhl J, Elkhoury JE, Ampuero JP, Tallakstad KT, et al. Average crack-front velocity during subcritical fracture propagation in a heterogeneous medium. *Phys Rev E* (2011) **84**:036104. doi: 10.1103/PhysRevE.84.036104
30. Lengliné O, Schmittbuhl J, Elkhoury JE, Ampuero JP, Toussaint R, Måløy KJ. Down-scaling of the fracture energy during brittle creep experiments. (2011) *J Geophys Res.* **116**:B08215. doi: 10.1029/2010JB008059

31. Lengliné O, Elkhoury JE, Daniel G, Schmittbuhl J, Toussaint R, Ampuero JP, et al. Interplay of seismic and aseismic deformations during earthquake swarms: an experimental approach. *Earth Planet Sci Lett.* (2012) **331**-332:215-23. doi: 10.1016/j.epsl.2012.03.022

32. Love AEH. The stress produced in a semi-infinite solid by pressure on part of the boundary. *Philos Trans R Soc London Ser A* (1929) **228**:377-420.

33. Rice JR. First-order variation in elastic fields due to variation in location of a planar crack front. *J Appl Mech*. (1985) **52**:571-9. doi: 10.1115/1.3169103

34. Gjerden KS, Stormo A, Hansen A. A model for stable interfacial crack growth. In: *Journal of Physics: Conference Series*. Vol. 402. Bristol, UK: IOP Publishing (2012). p. 012039.

35. Stormo A, Gjerden KJ, Hansen A. Onset of localization in heterogeneous interfacial failure. *Phys Rev E* (2012) **86**:025101. doi: 10.1103/PhysRevE.86.025101

36. Simonsen I, Hansen A, Nes OM. Determination of the Hurst exponent by use of wavelet transforms. *Phys Rev E* (1998) **58**:2779. doi: 10.1103/PhysRevE.58.2779

37. Family F, Vicsek T. Scaling of the active zone in the Eden process on percolation networks and the ballistic deposition model. *J Phys A Math Gen*. (1985) **18**:L75. doi: 10.1088/0305-4470/18/2/005

38. Tallakstad KT, Toussaint R, Santucci S, Schmittbuhl J, Måløy KJ. Local dynamics of a randomly pinned crack front during creep and forced propagation: an experimental study. *Phys Rev E* (2011) **83**:046108. doi: 10.1103/PhysRevE.83.046108

39. Hansen A, Schmittbuhl J. Origin of the universal roughness exponent of brittle fracture surfaces: stress-weighted percolation on the damage zone. *Phys Rev Lett.* (2003) **90**:045504. doi: 10.1103/PhysRevLett.90.045504

40. Duemmer O, Krauth W. Depinning exponents of the driven long-range elastic string. *J Stat Mech Theory Exp.* (2007) **2007**:P01019. doi: 10.1088/1742-5468/2007/01/P01019

41. Måløy KJ, Schmittbuhl J. Dynamical event during slow crack propagation. *Phys Rev Lett.* (2001) **87**:105502. doi: 10.1103/PhysRevLett.87.105502

42. López JM, Schmittbuhl J. Anomalous scaling of fracture surfaces. *Phys Rev E* (1998) **57**:6405.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer JB and handling Editor declared a current collaboration and the handling Editor states that the process nevertheless met the standards of a fair and objective review.

Copyright © 2016 Stormo, Lengliné, Schmittbuhl and Hansen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.