Brittle Crack Roughness in Three-Dimensional Beam Lattices

The roughness exponent is reported in numerical simulations with a three-dimensional elastic beam lattice. Two different types of disorder have been used to generate the breaking thresholds, i.e., distributions with a tail towards either strong or weak beams. Beyond the weak disorder regime a universal exponent of 0.59(1) is obtained. This is within the range 0.4-0.6 reported experimentally for small scale quasi-static fracture, as would be expected for media with a characteristic length scale.

The roughness exponent is reported in numerical simulations with a three dimensional elastic beam lattice. Two different types of disorder have been used to generate the breaking thresholds, i.e., distributions with a tail towards either strong or weak beams. Beyond the weak disorder regime a universal exponent of ζ = 0.59 ± 0.01 is obtained. This is within the range 0.4 ζ 0.6 reported experimentally for small scale quasi-static fracture, as would be expected for media with a characteristic length scale. The interest in crack morphology can be traced back to Mandelbrot et al., who introduced a mathematical framework for describing rough surfaces in terms of fractal geometry [1]. More specifically, crack surfaces have been shown to be self-affine in the sense that they satisfy certain scaling laws, where a characteristic exponent governs the asymptotic behaviour [2]. Since the quantitative properties of crack surfaces have been made readily accessible through the development of modern imaging techniques, predictions of computer models [3] can be compared with the results of practical experiments to verify or discard theoretical assumptions.
Much of the theoretical work done so far has been based on the random fuse model [4], i.e., a regular lattice of conducting elements which irreversibly burn out once their individual thresholds have been exceeded. Breakdown is driven by a potential difference between two opposing boundaries and the analogy of Kirchhoff's equations with linear elasticity is the reason why this model is referred to as a scalar model of fracture. The result obtained for the roughness exponent in simulations with the fuse model is ζ = 0.74(2) in two dimensions [5]. Recently, a somewhat different result, ζ = 0.84(3), has been reported [6]. Experimental studies in two dimensions have been carried out by Poirer et al. [7], who considered a stacking of collapsible cylinders (drinking straws) to obtain ζ = 0.73(7), by Kertesz et al. [8] for tear lines in wet paper, obtaining ζ ≈ 0.73, and by Engøy et al. [9] for crack lines in thin wood plates, yielding a roughness exponent of ζ = 0.68 (4). Results also differ according to the dimensionality, i.e., in three dimensions the fuse model result is ζ = 0.62(5) [10]. This result, however, is very different from the universal value suggested by Bouchaud et al. [11], i.e., ζ = 0.8, a value which has been confirmed in experimental studies [12]. Hence the apparent agreement of the fuse model with experiment in two dimensions is a coincidence.
An important point to be considered with the fuse model, however, is the fact that it models electrical breakdown rather than brittle elastic fracture. A different model which incorporates elasticity is the Bornmodel [13]. Here the material is modeled as a network of elastic springs, each spring being free to rotate at its ends. In two dimensions ζ = 0.64 (5) is obtained with this model [14]. As with the fuse model, the exponent drops when going from two to three dimensions. The result, ζ = 0.5 [15], is also significantly different from the universal value of ζ = 0.8. Recent findings, however, indicate that there should exist two different regimes for the scaling behaviour of cracks. As was first shown by Milman et al. [16], a much lower exponent ζ ≈ 0.5 is found on small length scales. Daguier et al. identified two self-affine fracture regimes, where ζ ∼ 0.5 is associated with small length scales and slow crack growth and ζ = 0.8 involves dynamic effects and is associated with larger length scales [17].
A problem with the Born model is that it is not rotationally invariant [18]. A different model which does not suffer from this drawback is the elastic beam model [19]. Here the beams are rigidly connected at the nodes so as to preserve the angle between any two neighbouring beams. Rotations thus induce flexing and twisting deformations, while linear displacements induce transverse shear and axial forces. The exponent obtained with the beam model is ζ = 0.86(3) in two dimensions [20], and shows that the scalar and vectorial descriptions of fracture need not necessarily produce similar results. Since the beam model provides a more realistic description of brittle elastic fracture than does either the fuse model or the Born model it is of great interest to see what the result is in three dimensions.
The lattice used in our calculations is a regular cube of size L × L × L, where each node is connected to its nearest neighbours by linearly elastic beams of unit length. Forces acting on the nodes have been deduced from the effect a concentrated end-load has on a beam with no end-restraints [21]. A coordinate system is placed on each node, and the enumeration of the connecting beams fol- lows an anti-clockwise scheme within the XY -plane, i.e., beginning with the beam which lies along the positive Xaxis and ending with that which extends upwards along the positive Z-axis, see Fig. 1.
At each stage of the breaking process, the updated displacements for each node is obtained from which is solved iteratively via relaxation using the conjugate gradient method. This minimizes the elastic energy to obtain those displacements for which the sum of forces and moments on each node vanish, i.e. the mechanical equilibrium.
Six terms contribute to each of the force components in Eq. (1). With the neighbouring nodes fixed, a translational displacement δx = x j − x i of node i induces axial forces in beams 1 and 3, and transverse forces in beams 2, 4, 5 and 6. With the expressions for force and moment defined in Refs. [20] and [22], this gives with similar expressions for Y i and Z i . Likewise, a rotational displacement δu = u j − u i about the X-axis causes a torsional moment in beams 1 and 3, and flexure in beams 2, 4, 5 and 6. Hence, now with similar expressions for V i and W i . To express the thirty-six force components in Eq. (1) more compactly, and are defined for notational convenience, to keep track of the signs. The Kronecker delta, moreover, has been used to construct and operators which include or exclude s and t from the sum over neighbouring beams. Eq. (2) can then be stated on the form where, exchanging indices (1,3) and (2,4), and interchanging v with u and x with y, Y i is obtained from X i . In the same way, Z i is obtained from X i by an exchange of the indices (1,3) and (5,6), by interchanging w with u and x with z, and by replacing r j with −s j . Next, Eq. (3) can be stated on a similar form, i.e., where V i , for angular displacements about the Y -axis, is obtained from U i by the same substitutions as in X i → Y i . Likewise, the expression for W i , relevant to angular displacements about the Z-axis, follows from U i by use of the X i → Z i substitutions, now however, without changing r j to −s j . Instead, the change from z to x is made so that −δx replaces δz. Breaking thresholds are generating by raising a random number r on the interval [0, 1] to a power D. The magnitude of D then corresponds to the strength of disorder, with the distribution having a tail towards weak beams for D > 0 and towards strong beams for D < 0. This is the most simple way of incorporating scale invariance, i.e., those mechanisms responsible for the multifractal behaviour of disordered systems, into the distribution of thresholds [23]. For each sample two casts are generated for the thresholds, one for the flexural strength and one for the axial strength.
Fracture is initiated by imposing a uniform displacement on all nodes defining the top surface of the cube. The first beam to break is the one with the weakest axial strength, whereupon the location of subsequent breaks depends on a complex interplay between quenched disorder and the constantly evolving non-uniform stress-field. Each time a beam is removed from the lattice, Eq. (1) is used to obtain the new equilibrium displacements. This is continued until a surface appears which divides the sample in two separate parts, see Fig. 2.
Results obtained for D = 1, 2 and 4 are shown in Figs. 3, 4 and 5, respectively, where the roughness W has been calculated for systems ranging in size from L = 3 to L = 40. Disregarding finite-size effects for small L, the relationship between L and W is clearly self-affine. Hence, defining the roughness as the root-mean-square of the variance perpendicular to the (average) fracture  plane, i.e., we have W ∼ L ζ , with the results of W x and W y being statistically the same. The exponent in Fig. 3 is ζ = 0.77(1), and does not belong to the self-affine regime of fracture. This is because D = 1 in a three dimensional lattice generates too little disorder. The qualitative difference between the fracture surfaces obtained for D = 1 and higher values of D is shown in Fig. 2. A universal exponent emerges as the disorder is increased, however, with the exponents obtained for D = 2 and D = 4 being ζ = 0.62(2) and ζ = 0.592(5), respectively.
The result also seems to hold for disorders with a tail towards strong beams. For D = −4 the exponent is ζ = 0.65(2), see Fig. 6. Although this is slightly different from the best value obtained for D > 0, which we take to be ζ = 0.59(1), the small discrepancy is probably due to the fact that D = −4 lies just within the transient regime towards strong disorders for distributions of type D < 0. In two dimensions the D < 0 stress-strain curve displays a cross-over towards stable crack-growth when |D| ∼ 2. This cross-over value is probably higher in the cube lattice, due to the stronger constraints imposed on the crack path in three dimensions.
To summarize, the roughness exponent ζ = 0.59(1) is obtained numerically for a cubic lattice of elastic beams. This is lower than the experimental value relevant to large scales, ζ = 0.8, but within the range reported for small scales, 0.4 < ζ < 0.6. This lends further evidence to the existence of a different scaling regime for slow crack growth involving characteristic, or "small", length scales [25]. The agreement of our result with the exponents reported in this regime stems from the discretization in terms of a beam lattice, i.e., a system which involves a characteristic intermediate length scale, as in the case of Cosserat elasticity.