Local dynamics of a randomly pinned crack front: A numerical study

We investigate numerically the dynamics of crack propagation along a weak plane using a model consisting of fibers connecting a soft and a hard clamp. This bottom-up model has previously been shown to contain the competition of two crack propagation mechanisms: coalescence of damage with the front on small scales and pinned elastic line motion on large scales. We investigate the dynamical scaling properties of the model, both on small and large scale. The model results compare favorable with experimental results.


I. INTRODUCTION
The motion of a fracture front through a disordered medium is for obvious reasons of great practical importance. Nevertheless, it is a very complex problem which have kept materials scientists and mechanical engineers occupied for almost two centuries [1]. In an attempt to simplify the problem, Schmittbuhl et al. [2] proposed to study the motion of a fracture front moving along a weak plane, thereby supressing the out of plane motion of the front. In a seminal experimental study, Schmittbuhl and Måløy [3] followed this idea up by sintering two sandblasted plexiglass plates together and then plying them apart from one edge. The sintering made the otherwise due to the sandblasting, opaque plates transparent. Where the plates were broken apart again, the opaqueness returned. Hence, it was possible to identify and follow the motion of the front since the plates would be transparent in front of the crack front and opaque behind it. There have been a large number of studies on this system following this initial work. Most of these studies have been theoretical or numerical in character and they have mostly dealt with the roughness of the crack front, see e.g. [4][5][6][7][8]. See also the recent reviews by Bonamy and Bonamy and Bouchaud [9,10].
The roughness is only one aspect of the motion of the crack front in this system. A first study of the dynamics of the front was published by Måløy and Schmittbuhl in 2001 [11]. The velocity distribution was measured and found to decay slower than an exponential. This was followed by a study by Måløy et al. [12] introducing the waiting time matrix technique making it possible to measure pixel by pixel the time the front sits still locally. This revealed an avalanche structure where portions of the front would sweep an area S with a velocity v before halting again. Both the distribution of areas and velocities were power laws. In two later papers, Tallakstad et al. [13] and Lengliné et al. [14], have continued and refined this work.
The fluctuating line model, based on the idea that the crack front moves as a pinned elastic line [2,15], dominates the theoretical descriptions of the constrained crack problem. It is a top-down approach where the motion of the crack front is derived from continuum elastic theory. Its earliest use was to calculate the roughness exponent controlling the scaling properties of the crack front. If h(x, t) is the position of the crack front at time t and position x along its base line, then one finds the heightheight correlation function (h(x + ∆x) − h(x)) 2 scaling as |∆x| ζ where ζ is the rougness exponent [2]. The most presise measurement of the roughness exponent withing the fluctuating line model to date is ζ = 0.388 ± 0.002 [6]. However, the experimental measurements of the ζ gave systematically a much larger value, namely around 0.6 -Schmittbuhl and Måloy [3] found ζ = 0.55 ± 0.05 and Delaplace et al. [5] ζ = 0.63 ± 0.03.
Santucci et al. [8] has reanalyzed data from a number of earlier studies, including [5], finding that the crack front has two scaling regimes: one small-scale regime described by a roughness exponent ζ − = 0.60 ± 0.05 and a large-scale regime described by a roughness exponent ζ + = 0.35 ± 0.05. Hence, there is a large scale roughness regime which is consistent with the fluctuating line model. However, the small scale regime is too rough. Laurson et al. [16] has linked this regime to correlations in the pinning strength below the Larkin lenght within the fluctuating line model. This approach assumes that the physics behind both scaling regimes is describable within the same top-down model.
A very different approach to explain the roughness of the crack front has been proposed by Schmittbuhl et al. [7]. It is a bottom-up approach based on a stressweighted gradient percolation process. This idea in turn has its origin in the proposal by Bouchaud et al. [17] that the crack front does not advance not only due to a competition between effective elastic forces and pinning forces at the front, but also by coalescence of damage in front of the crack with the advancing crack itself. The coalescence By refining the model proposed by Batrouni et al. [18,19], consisting of fibers clamped between a hard and a soft block and with a gradient in the breaking thresholds, we have identified two scaling regimes for the roughness of the advancing crack front [20]. On large scales we recover the roughness seen in the reanalysis of experimental data by Santucci et al. [8], ζ = 0.39 ± 0.04, consistent with the fluctuating line model, whereas on small scales with identify a gradient percolation process, leading to ζ = 2/3 [21], which is also consistent with Santucci et al. reanalysis.
We study in this paper the crackling dynamics in the fiber bundle model of Gjerden et al. [19,20]. We find quantitative consistency with the experimental measurements. Hence, we demonstrate that this simple model contains the observed experimental features of constrained crack growth. It contains the essential physics of the problem. Since we have control over the crossover length scale between the crack advancing due to crack coalescence and due to pinning of the crack front as an elastic line, we analyse our results in light of this.
In Section II, we describe the model in detail. Section III continues the analysis from Ref. [20] of the roughness of the crack front. In [20], we measured the roughness locally using the average wavelet coefficient method [24,25]. Here we base our analysis on the variance of the front, again finding two regimes: a small-scale regime consistent with ordinary gradient percolation and a largescale regime consistent with the fluctuating line model. In Section IV we examine the velocity distribution during the motion of the crack front. Our data are consistent with the findings of Måløy et al. [12] and Tallakstad et al. [13]. Section V is devoted to an analysis of the geometry of the avalanches that occur during the motion of the crack front. The last section before the conclusion, Section VI contains a tying up of loose ends through a discussion of the results.
The main aim of this paper has been to demonstrate that the simple model we use is capable of reproducing the experimental data. Hence, this simple fiber bundle model seems to contain the essential physics of the constrained crack problem.

II. MODEL
The model we base our calculations on is a refinement of the fiber bundle model [26] used by Schmittbuhl et al. [7] and introduced by Batrouni et al. the year before [18]. L × L elastic fibers are placed in a square lattice between two clamps. One of the clamps is infinely stiff whereas the other has a finite Young modulus E and a Poisson ratio ν. All fibers are equally long and have the same elastic constant k. We measure the position of the stiff clamp with respect its position when all fibers carry zero force, D. The force carried by the fiber at position (i, j), where i and j are coordinates in a cartesian coordinate system oriented along the edges of the system, is then where u (i,j) is the elongation of the fiber at (i, j) and k is the spring constant of the fibers. It is assumed to be the same for all fibers. The fibers redistribute the forces they carry through the response of the clamp with finite elasticity. The redistribution of forces is accomplished by using the Green function connecting the force f (m,n) acting on the clamp from fiber (m, n) with the deformation u (i,j) at fiber (i, j), [27] where a is the distance between neighboring fibers, ν the Poisson ratio and E the elastic constant. The integration in Eq. 2b is performed over the voronoi cells around each fiber. r (i,j) (0, 0) − r (m,n) (0, 0) is the distance between fibers (i, j) and (m, n). The equation set, Eqs. (2) and (2b), is solved using a Fourier accelerated conjugate gradient method [28,29].
The Green function, Eq. (2b), is proportional to (Ea) −1 . The elastic constant of the fibers, k, must be proportional to a 2 . The linear size of the system is aL. Hence, by changing the linear size of the system without changing the discretization a, we change L → λL but leave (Ea) and k unchanged. If we on the other hand change the discretization without changing the linear size of the system, we simultaneously set L → λL, (Ea) → λ(Ea) and k → k/λ 2 . We define the scaled Young modulus e = (Ea)/L. Hence, changing e without changing k is equivalent to changing L -and hence the linear size of the system -while keeping the elastic properties of the system constant [30].
The fibers are broken by using the quasistatic approach [31]. That is, we assign to each fiber (i, j) a threshold value t (i,j) . They are then broken one by one by each time identifying max (i,j) (f (i,j) /t (i,j) ) for D = 1. This ratio is then used to read off the value D at which the next fiber breaks.
In the constrained crack growth experiments of Schmittbuhl and Måløy [3], the two sintered plexiglass plates were plied apart from one edge. In the numerical modeling of Schmittbuhl et al. [7], an asymmetric loading was accomplished by introducing a linear gradient in D. Rather than implementing an asymmetric loading, we use a gradient in the threshold distribution, t (i,j) = gj + r (i,j) , where g is the gradient and r (i,j) is a random number drawn from a flat distribution on the unit interval. In the limit of large Young modulus E, i.e., both clams are stiff, the Green function vanishes and the system becomes equivalent to the gradient percolation problem [32,33].
In order to follow the crack front as the breakdown process develops, we implement the "conveyor belt" technique [19,34] where a new upper row of intact fibers is added and a lower row of broken fibers removed from the system at regular intervals. This makes it possible to follow the advancing crack front indefinitely. The implementation has been described in detail in [19].

III. ROUGHNESS AND DYNAMICAL SCALING
We demonstrate in this Section the existence of two scaling regimes for the roughness of the crack front by providing evidence beyond that which was presented in Gjerden et al. [20].
We identify the crack front by first eliminating all islands of surviving fibers behind it and all islands of failed fibers in front of it. We measure "time" n in terms of the number of failed fibers. After an initial period, the system settles into a steady state. We then record the position of the crack front j = j(i, n = 0) after having set n = 0. We then define the position at later times n > 0 relative to this initial position, This is the same definition as was used by Schmittbuhl and Måløy [3]. The front as it has now been defined will contain overhangs. That is, there may be multiple values of h i (n) for the same i and n values. We only keep the largest h i (n), i.e. we implement the Solid-on-Solid (SOS) front.
We define the average position of the front as and the front width as We start by examining the front in the range of values of the scaled elastic constant e for which the front has the character of a gradient percolation hull, e = 3.1 [20]. Fig.  1 shows one-parameter scaling of the front width w versus the average height of the front h with the gradient g andas the scaling parameter. There is data collapse for a range of system sizes and gradients. Based on measuring the fractal dimension of the front and using the average wavelet coefficient method to determine the roughness exponent, we concluded in [20] that the front roughness is in the ordinary gradient percolation universality class. This implies that the two scaling exponents in Fig. 1 are α = β = ν(1 + ν) = 4/7, where ν = 4/3, the percolation value [21]. We now assume that the front width obeys Family-Vicsek scaling [35], i.e., where For h ≪ L κ , and assuming that we are dealing with ordinary percolation, we may relate h to the number of fibers below and including the front, A, by the fractal dimension D A = 91/48 [36], h DA ∼ A. Invoking Mandelbrot's area-perimeter relation, we may also relate the length of the front, l to the number of fibers below and including the front, [33,37]. Eliminating A between these two expressions leads to l ∼ h DH DA/2 . We relate the length of the front l to its width w by the relation l ∼ w DH (L/w) [20]. Eliminating l between these two relations gives where 288/637 ≈ 0.45. The straight line in Fig. 1 has this slope. We may now calculate the dynamical exponent κ in Eq. (6). We find κ = 637/432 ≈ 1.47. Fig. 2 shows the corresponding one-parameter scaling as in Fig. 1, but now for systems with a much lower e = 7.8 × 10 −4 , placing us in the realm of the fluctuating line model [20]. The two scaling exponents that produce data collapse are α = β = 0.4. By a least squares fit, we determine w ∼ h 0.52±0.05 . This value is very close to the one seen for large e values, and consistent with the measurement of Tallakstad et al. [13], who found 0.55 and Schmittbuhl et al. [38] who found 0.52, both groups using sintered PMMA plates. Assuming Family-Vicsek scaling, we determine the dynamical exponent to be κ = 0.75±0.07. This is consistent with the value reported by Duemmer and Krauth for the fluctuating line model, κ = 0.70 ± 0.5.
The dynamical exponent reported by Måloy and Schmittbuhl [11] for the PMMA system was κ = 1.2. We find κ = 1.47 on small scales and κ = 0.75 on large scales.
It was observed by Santucci et al. [8] that the largescale roughness regime shows a gaussian distribution of the height fluctuations, whereas in the small-scale it broadens to a non-gaussian distribution. We investigate this in our model. We define the height fluctuations as ∆h(1) = |h i+1 (n) − h i (n)|. By changing the elastic constant, E in our model while keeping the step size constant, the result is equivalent to changing the step size in the experiment. We show in Fig. 3 the distribution of fluctuations for different E. We see that for soft systems, the distribution is gaussian. As the system stiffens, the distribution broadens as in the experiment.

IV. LOCAL VELOCITY DISTRIBUTION
We examine the distribution of local velocities using the waiting time matrix (WTM) method [12]. The WTM method consists of letting the fracture move across a matrix with all elements equal to zero in discrete time steps. Every time step, the matrix elements that contain the front are incremented by one. When the front has passed, the matrix contains the time the front spent in each position. Hence, matrix elements contaning low values correspond to rapid front movements and large values correspond to slow, pinned movement. From the WTM, one can then calculate the spatiotemporal map of velocities by mapping the coordinates of the front at each step in the simulation onto the local velocities V (x, y) calculated from the inverse of the waiting time, thus obtaining v(x, t). From this map we then calculate the global average velocity v and the distribution of local velocities P (v). It has been found experimentally that this distribution follows with the exponent η = 2.55 ± 0.15 [12,13]. The results from our simulations are shown in Fig. 4 for scaled elastic constant e = 1.6 × 10 −4 and 7.8 × 10 −4 . The slope in the figure is obtained through a linear fit and has the value 2.53. Based on this figure, we follow Tallakstad et al. [13] and define a pinning regime where v < v and a depinning regime where v > v .

A. Space and time correlations
From the local velocities, we calculate the normalized correlation functions in space, G(∆x), and time, G(∆t) (10b) σ x/t and x/t denote standard deviation and average over x or t. The results for G(∆x) is presented in figure  5. The results are in qualitative agreement with what has been reported experimentally [13]. Due to the small system sizes available numerically, the cutoff we observe is extremely sharp. Still, the data appears to be consistent with a power law distribution with an exponential cutoff G(∆x) = A(∆x/L) −τx exp(−∆x/0.1L) with τ x = 0.4. This is in the same range as the values reported by Tallakstad et al. [13], where τ x = 0.53 ± 0.12.
In figure 6, we show the results of the analysis of the temporal velocity distribution. Again, we attempt to fit the data to a power law with exponential cutoff, G(∆t) = B(∆t v ) −τt exp(−∆t v /0.25), with exponent τ t = 0.43. The scaling resemble closely that observed by Tallakstad et al. [13].
We note that data collapse is obtained in Fig. 5 when the correlation function is scaled by g 1/2 and ∆x by L, and in Fig. 6, the correlation function is scaled by gL 1/2 .

V. CLUSTER ANALYSIS
In the following we analyze the geometrical structure of the areas swept by the crack front during avalanches.
Following Tallakstad et al. [13], we construct thresholded velocity maps V C in such a way that in the depinning regime and in the pinning regime. C is a constant in the range of 2−12 for depinning and 1.5−6 for pinning for our simulations. For larger system sizes, larger C are possible. Examples of these matrices are shown in figure 7, where the structural difference between depinning and pinning clusters can be seen. In the spirit of de Saint-Exupéry we may characterize the pinning clusters as snakelike, whereas the depinning clusters are elephant-in-snakelike [40]. The following analysis uses only L = 64 and L = 128, and are based on 250 and 100 samples, respectively. g = 0.025 for the smaller systems, and g = 0.0125 for the larger systems. For L = 64, e = 1.6 × 10 −4 , and for L = 128, e = 7.8×10 −4 . The systems will for these e-values behave virtually identical [20].

A. Size distribution of clusters
We denote the number of bonds in a cluster the size or area of a cluster S. Both experiments [12,13] and numerical simulations using the fluctuating line model [22] on in-plane fracture, show a power law distribution P (S) ∼ S −γ for S with exponent close to γ = 1.56 [13] or γ = 1.65 [22].
We show in Fig. 8 the probability distribution of S both for pinning and depinning for scaled elastic constant e = 1.6 × 10 −4 . Our data are consistent with previous measurements, but due to the limited size of the numerical simulations currently available, the exponent cannot be accurately determined. One prominent feature of the clusters is that due to the small system size, the presence of large clusters for low values of C suppresses the number of smaller clusters, and for high values of C, the number of large clusters is reduced by the filtering process in C. This is a pure finite size effect, which we assume will disappear for larger system sizes. The experiments, which would correspond to a very large simulation size, shows no such dependence on C. A direct comparison between the experimental and our data is given in Figure 9.

B. Aspect ratio of clusters
Next, we extract the linear extensions as length and height from the cluster sizes. This is done by enclosing each cluster in a bounding box of size l x by l y . Although the appearance differs greatly between pinning and depinning clusters, the method for determining l x and l y are the same. As seen qualitatively in Fig. 7, pinning clusters are much longer than they are wide, and in a worst-case scenario they could approximate the diagonal in the bounding box, which would lead to a gross overestimation of l y . However, such behavior was not seen for L = 64, 128 and we judge the approximation of l y to be equally sound in both regimes. Figure 10 shows the size distributions of l x and l y for scaled elastic constant e = 1.6 × 10 −4 . As in Fig. 8, there is a pronounced dependence on C. The data indicate power law distributions P (l x ) ∼ l −βx The longest data series (unfilled markers) are experimental data obtained with permission from the authors [13]. The shorter data series (filled markers) are numerical data from Fig. 8. The experimental data are unchanged from [13], so S ′ = S * γ and S * are as given in [13]. Due to lack of direct physical parameter comparison, the numerical data are scaled to match the experimental data to verify similar behavior.
1.93 and β y = 2.36 [13]. Both pinning and depinning data show behaviour consistent with the experimental data.
Turning to examine the relation between l y and l x , we plot these against each other in Fig. 11. This data is more well-behaved than the probability distributions, so scaling relations are more pronounced. The data are consistent with a relation of type as studied both by Tallakstad et al. [13] and Bonamy et al. [22] in respectively the experimental system and with the fluctuating line model. We see a clear dependency of H on parameters such as e and C. Considering the depinning case first for e = 7.8 × 10 −4 , we divide the data between C high ∈ [6 − 12], a group containing data from only the highest local velocities, and C low ∈ [3 − 5], containing lower local velocities. There is a clear visual difference in behavior between the two groups: The group with the higher C values, have the higher slope. We measure the corresponding exponents to be H high d = 0.9 and H low d = 0.65. In the pinning regime we similarly divide between low and high local velocity effects, except that the range of available C is much smaller, so we find C high = 2, 3 and C low = 1. This yields the two exponents H high p = 0.65 and H low p = 0.4. It has been suggested that H is another measure of ζ, the roughness exponent [12,16,22]. In order to test this, we have added data both for pinning and depinning of a stiff system (e = 7.8) in Fig. 11. If the conjecture is right in this case, we would expect a H exponent equal to 2/3. We have added such slopes, and indeed this seems fulfilled. We may then speculate whether the exponents that we measure for the softer system (e = 7.8 × 10 −4 ), indeed are roughness exponents -at least for some groups of C values.
Another way to examine these scaling relations is to check the relation between l x , l y and S. We expect to see a power law relationship of type The data are plotted in Fig. 12, in which data for the highest local velocities are not included, C = C low for both pinning and depinning. In the depinning regime, we find the exponents α d x = 0.6 and α d y = 0.4. Eliminating S, we get a measure of H through resulting in H d = 0.65. In the pinning regime we find α p x = 0.75 and α p y = 0.3, yielding H p = 0.4. These results are to be compared to the experimental measurements [13] giving α x = 0.62 ± 0.04 in pinning and the depinning regime and α y = 0.41 ± 0.06 or 0.34 ± 0.05 in the depinning regime depending on how the bounding box was defined. This leads to the value H = α y /α x = 0.66 or 0.55 depending on the bounding box used. Bonamy et al. [22] found H = 0.65 ± 0.05.

VI. DISCUSSION
In our model, e and g are the most important input parameters. g effectively controls the loading conditions. A high value for g amounts to a large angle of contact between the materials, or slow loading conditions, and vice versa. As long as 1/g is in a range where it allows the front to develop but not extend out of the system, we see no change in any scaling exponents we have investigated. This behavior is supported by the experimental results where no dependence on the loading conditions was found and the average propagation velocity of the front spanned more than four and a half orders of magnitude [13].
The parameter e controls the material behavior. A small value for e can amount to model a system of low modulus of elasticity or a long length scale, and vice versa. Based on the behavior alluded to in our model, we may propose the following scenario. For low values of e, the fracture process is due to stress concentration leading to damage forming on the fracture front. This regime is classified by a roughness exponent of ζ large = 0.39 and an underlying burst process of same exponent H = ζ large . As e increases, the burst process changes behavior to H = ζ small = 0.67, while the front is still characterized by ζ large . As e is increased further, the effect of the burst exponent begins to show on the scale of the front and damage begins to form ahead of the front. At this point, the roughness exponent changes from 0.39 and grows into 0.67. At higher e, both the burst process and the roughness of the front is now characterized by the larger exponent H = ζ small = 0.67.

VII. CONCLUSION
We have developed a minimal model which we have previously shown to contain two different fracture growth processes which reproduce results from experiments with fracture fronts of roughness ζ small = 0.67 on small length scales and ζ large = 0.39 on larger length scales [20]. In this paper we have shown that this quasistatic model also produces data consistent with the experimental analysis of the dynamics of the fracture front. Most notably we recover the existence and behavior of a depinning and a pinning regime. In the pinning regime we recover the same 0.39-scaling exponent of the front in the bursts. In the depinning regime we can see traces of this behavior, but more importantly, we recover the 0.67-scaling exponent of the front seen on shorter length scales. From this, we conjecture that the collective behavior of the underlying bond breaking process is first and most visible in the damage bursts. The global front roughness changes according to the change in burst process when one process suppresses the other, i.e., on a different scale.