Discrete element modeling of brittle crack roughness in three dimensions

Crack morphology obtained in the fracture of materials with a disordered micro-structure is studied using numerical simulations. Physical properties are embedded on a regular three dimensional lattice as discrete stochastic elements which conform to the laws of linear elasticity. In this model, also known as the beam lattice, these elements are analogous to beams in that relative displacements between neighbouring nodes induce axial, bending and shearing forces, as in a real elastic solid. The stochastic nature enters via the introduction of random breaking thresholds on the individual elements. Using this model, the exponent characterizing the scaling with system size of the crack roughness perpendicular to the fracture plane is reported. Two different types of disorder have been used to generate the thresholds, i.e., distributions with a tail towards strong elements or with a tail towards weak elements. At weak disorders the self-affine regime seems to lie beyond the system sizes presently included. At stronger disorders a self-affine regime appears, for which we obtain exponents consistent with 0.6 for both types of disorder. The latter result is in fair agreement with the experimental value reported for large length scales, 0.50.


INTRODUCTION
The modern academic 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]. 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 behavior [2]. This means that a scaling of lengths with a factor λ within the fracture plane implies that lengths perpendicular to this plane scale with a factor λ ζ . Understanding fundamental aspects of crack roughness is also important in the engineering fields. Such is the case in the oil industry where flow properties within fractured rock masses is of great interest. Quantitative properties of crack surfaces are accessible through modern imaging techniques, and predictions of computer models [3,4] 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 [5], i.e., a regular lattice of conducting elements which irreversibly burn out once their individual thresholds have been exceeded. Breakdown is driven by a voltage 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. In two dimensions results reported for the roughness exponent with the fuse model are ζ = 0.74 (2) [6,7], and more recently ζ = 0.83 (4) [8].  (2) [10,11], and Nukala et al. obtained ζ = 0.52(3) [12].
How do these results compare with experimental findings? Initially it was thought that the universal roughness exponent should be ζ 0.8 [2]. Later it was realized that there should instead exist two different regimes for the scaling behavior of cracks [13,14]. It was proposed that a smaller exponent ζ 0.5 should be associated with short length scales and slow crack growth [15], while a larger exponent ζ 0.8 involving dynamic effects of fracture should be associated with longer length scales. The picture that has emerged more recently however is that, although there are indeed two scaling regimes, it is in fact the larger exponent ζ 0.8 which is relevant to short length scales. This is because of the screening of elastic interactions which takes place within the intermediate region between the crack tip and the surrounding medium, i.e., the fracture process zone. At larger length scales one observes ζ 0.45 [16,17].
Roughly speaking, then, the three dimensional fuse model results are in agreement with what is obtained experimentally at large length scales. An important point to be considered with the fuse model, however, is the fact that it models electrical breakdown rather than fracture in a fully elastic medium. A different lattice model which incorporates elasticity is the Bornmodel [18]. 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 [19]. As with the fuse model, this exponent is seen to drop somewhat when going from two to three dimensions. The result, ζ = 0.5 [20], is consistent with the range of exponents obtained with the random fuse model and agrees rather well with the scaling expected at length scales above those relevant to the fracture process zone.
A problem with the Born model is that it is not rotationally invariant [21]. A different model which does not suffer from this drawback is the elastic beam model [22,23]. Here the beams are rigidly connected at the nodes so as to preserve the angle between any two neighboring beams. Rotations thus induce flexing and twisting deformations, while linear displacements induce transverse shear and axial forces. The exponent obtained with the beam model in two dimensions is ζ = 0.86(3) [24], and shows that the scalar and vectorial descriptions of fracture need not necessarily produce similar results, although the more recent two dimensional random fuse result of ζ = 0.83(4) quoted by Zapperi et al. [8] is consistent with this. 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. An analogous study to this end has been made by Nukala et al. [25], obtaining both the global roughness exponent ζ = 0.48(3) and the local roughness exponent ζ loc , which they find to be equal to the global roughness exponent for the beam lattice, ζ = ζ loc . The local exponent in such cases is obtained by considering the scaling properties within "windows" of the crack surface for length scales l << L, where L is the external size of the L × L × L lattice. Presently we only regard the global roughness exponent.

MODEL
Our model is a deformable lattice in the form of a regular cube with size L × L × L, where each node is connected to its nearest neighbors by linearly elastic beams. Forces acting on the nodes have been deduced from the effect a concentrated end-load has on a beam with no end-restraints [26,27]. A coordinate system is placed on each node, and the enumeration of the connecting beams follows an anti-clockwise scheme within the XY-plane, i.e., beginning with the beam which lies along the positive X-axis and ending with that which extends upwards along the positive Z-axis, see Figure 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 [28,29]. 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. In Equation (1), x i , y i , and z i are the coordinate displacements of node i relative to its starting position before fracturing is initiated. Likewise, u i , v i , and w i are the angular displacements around the X-, Y-, and Z-axes, respectively (see Figure 1). Presently we use the same expressions for force and moment as those used in Herrmann et al. [23] and Skjetne et al. [24]. We thus have three constants defined as where E and G are Young's modulus and the shear modulus, respectively, A is the area of the cross section of the beam, is the beam length and I is the moment of inertia about the centroidal axis of the beam. Our choice of these parameters mirrors that of Herrmann et al. [23] and Skjetne et al. [24], i.e., the beam length is set to = 1 and we use α = 1, β = 30/7, and γ = 60/7. Additionally, we define the quantity where J is the moment of inertia for torsion, in analogy with Equation (2) since in three dimensions we require torque in addition to bending, axial and transverse forces. In this case we use = 1. In the following it is also convenient to define and similarly for the other five displacement coordinates. Six terms contribute to each of the force components in Equation (1). For instance, if we imagined the neighboring nodes to be fixed, a translation x i of the central node i would induce axial forces in beams 1 and 3 and transverse forces in beams 2, 4, 5, and 6. If we take into account the displacements of the neighboring nodes as well, the axial force on node i from beam 1 is while the transverse force on node i from beam 2, along the Xaxis, is given by In each case j refers to the neighbor depicted in Figure 1. Consequently, is how the force on node i along the X-axis depends on the displacements and rotations of its six nearest neighboring nodes. Similar expressions are deduced for Y i or Z i by considering translations along the Y-axis or the Z-axis, respectively. Likewise, an angular displacement u i about the X-axis with the neighboring nodes fixed would create torque in beams 1 and 3, and bending in beams 2, 4, 5, and 6. More generally, the torque in node i from beam 1 is while the bending moment from beam 2 is For the angular force on node i about the X-axis and its dependence on the displacements of the six neighboring nodes, we have now with similar expressions for V i and W i . To express the 36 force components in Equation (1) more compactly, and s j = −1 j r j (12) are quantities which we define for notational convenience, to keep track of the signs and contributions from neighboring beams.
The j in each case refers to the neighboring beams as shown in Figure 1. The Kronecker delta, moreover, has been used to construct i.e., an operator which includes s and t in the sum over all nearest neighbors, and which instead excludes s and t from the sum over all nearest neighbors.
For the six components making up the force on node i along the X-axis, i.e., Equation (7), we can now state this on a compact form as and Y i as In the same way, Z i becomes Next, Equation (10) for angular displacements about the X-axis is written out in full as and V i , for angular displacements about the Y-axis, becomes Lastly, for angular displacements about the Z-axis, we get www.frontiersin.org November 2014 | Volume 2 | Article 68 | 3 To include the effects of material disorder we generate a random number r on the unit interval [0, 1] and let this represent the cumulative threshold distribution. We then assign thresholds according to t F = r D , using two different types of distribution. In the case of D > 0 the distribution is a power law with a maximum threshold of 1 which has a tail extending toward zero. In the other case, D < 0, the power law is limited downwards by a minimum threshold of 1, except now it has a tail extending up toward infinity. The respective cumulative distribution functions are then given by for D > 0, where 0 ≤ t F ≤ 1, and by for D < 0, where 1 ≤ t F < ∞. The case of D = 0 corresponds to all thresholds being the same (t F = 1), i.e., we have a homogeneous medium without structural disorder. An increase in the magnitude of the exponent |D| causes the coefficient of variation with respect to any two random numbers r and r on the interval [0, 1] to increase. For the distribution types 0 ≤ t F ≤ 1 and 1 ≤ t F < ∞, the coefficient of variation for any two random numbers within the interval are reciprocal but otherwise the same. Therefore, large values of |D| correspond to strong disorders and small values to weak disorders. This prescription is the most simple way of incorporating scale invariance into the distribution of thresholds. As system size diverges, only power law tails toward zero or infinity should be important [30][31][32]. The use of D as a parameter is then very convenient and enables the asymptotic behavior of the fracture process to be fully explored as a function of the disorder.
For each sample two thresholds are generated for each beam on the lattice, one for the amount of bending that particular beam can sustain before breaking (in pure bending mode), and one for the amount of axial stress it can sustain before breaking (in pure axial loading). Each beam is assumed to be linearly elastic up to the breaking threshold. The two thresholds are combined in an interaction formula for combined loading, i.e., where t A and t M are the axial and bending thresholds, respectively, and where M is the largest of the bending moments at the two beam ends i and j. Essentially this is the same failure criterion as that used in Herrmann et al. [23], except that |M| here generally represents the largest of the values obtained by comparing bending within two planes. Hence, in the case of beams 1 and 3, where xy denotes bending within the XY-plane, etc. In these calculations we have not considered breaking due to torque. The engineering literature suggests a range of interaction formulae relevant to combined loading, but we have only used the one outlined above. 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 nonuniform stress-field. Each time a beam is removed from the lattice, Equation (1) is used to obtain new equilibrium displacements before Equation (23) is used to identify the next beam to fail. This is continued until a surface appears which divides the sample in two separate parts, see Figure 2. In this process one beam at a time is removed from the lattice and the redistribution of stress is assumed to occur much faster than the breaking of individual beams.

RESULTS AND DISCUSSION
In our calculations we have used open rather than periodic boundary conditions. This introduces edge effects, especially for small systems, but we assume it does not affect the roughness calculations significantly for our largest systems. For each trace along the surface we subtract the drift in the vertical direction, and overhangs which might appear are removed. The roughness is obtained as the root-mean-square of the variance perpendicular to the (average) fracture plane, i.e., where z x (i) is the height in the vertical direction of the first intact node encountered when moving down toward the lower remaining part of the cube (those parts shown in Figure 2). Tracing parallel to the X-axis we obtain W x , and parallel to the Y-axis we obtain W y , results for W x and W y are statistically the same. We have evaluated four different disorders, having made three sets of calculations for type D > 0 disorder, i.e., D = 1, D = 2, and D = 4, and one set for type D < 0 disorder, i.e., D = −4. In each case a large number of samples were calculated for each size. The system sizes we show in the plots are from L = 3 up to L = 40 or L = 50, and for D = 1 we include L = 63. In a log-log plot the self-affine relationship between the outof-plane crack roughness W and the system size L shows up as a straight line, and we determine the roughness exponent from the slope of this line. Apart from the smallest systems, for which edge effects are significant, a self-affine relationship W ∼ L ζ emerges well within the range of system sizes included in the present study if we regard D = 2, D = 4, and D = −4. Presently system sizes for which we have a sufficiently large number of samples lie in the range L = 3 to about L = 50. Results for D = 2 and D = 4 are shown in Figures 3, 4, respectively. The exponent obtained for D = 2 is ζ = 0.62(2), based on a straight-line fit to systems L = 10, 13, 16, 19, 20, 22, 25, 32, and 40. For D = 4 we obtain ζ = 0.590(4), based on a straight-line fit to L = 10, 13, 16, 20, 25, 32, 40, and 50. At weak disorders the change in the fracture process from disorder dominated fracturing, where small cracks appear randomly throughout the material, to stress dominated fracturing, where a large crack emerges and breaks the sample, occurs at an early stage. Hence the computer time needed to break such each sample is less than that required to break a sample at strong disorder. Despite this, a very much larger number of samples is needed to obtain solid statistics at weak disorders. A comparable number of samples was generated for D = 2 and D = 4, but the statistics are stronger in the latter case.
In the case of D = 1 it is perhaps more questionable whether a self-affine regime has appeared within the range of system sizes studied, as can be seen from Figure 5. In Figure 5 a line with slope ζ = 0.78 has been fitted to the values for L = 16, 20, 25, and 32, corresponding to an apparently linear regime. This is seen to cross over into a regime with a smaller roughness exponent. At present, however, we cannot tell what this exponent is, based on the system sizes used currently. For the number of samples calculated, the average values for L = 50 and L = 63 are subject to a certain amount of statistical fluctuation but these are not large enough to place them on the fitted slope. We assume that a self-affine regime manifests itself for larger system sizes.
For weak disorders the fracture surface will appear as mainly flat plateaus, slightly stepped up or down in relation to each other, due to the surface height variations being smaller than the lattice course-graining. This effect is especially accentuated for small system sizes, but improves as system size increases. Because of this, the length scale at which the cross-over occurs decreases with increasing disorder and for this reason it appears that larger system sizes are required in order to reliably estimate the roughness exponent for weak disorders. So although D = 1 represents a comparatively strong disorder in two dimensions [24], where the resulting crack displays a very rough surface, this is not the case in three dimensions. In three dimensions the added constraint on the crack surface from the extra dimension causes  the two dimensional traces along the three dimensional surface to be more smooth than what is the case for a trace along a purely two dimensional crack. Consequently, as can also be seen by examining the fracture surfaces shown in Figure 2, D = 1 in the threshold distribution corresponds to weak disorder in three dimensions. Hence, the inclusion of larger system sizes, with sufficiently many samples, is required to determine the exponent for D = 1 or weaker disorders.
We next consider fracture in the D < 0 regime. At weak disorder one has a more or less uniformly distributed strength throughout the elastic medium, interrupted here and there by the odd strong element. Fracture in this case is stress dominated from the onset-as the crack advances the next beams encountered along the path are most likely comparable in strength to those already broken and since the stress is highest at the crack tip these are the ones that will break next. Consequently, crack growth is localized and unstable from the beginning, resulting in a fracture surface that is quite smooth. As D increases, however, a region of stable crack growth emerges, since now the number of strong beams in the distribution of thresholds becomes significant. Fracture is now disorder dominated from the beginning and a rougher interface results as the randomly distributed micro-cracks gradually merge into a macroscopic fracture surface. Finally, at some point, fracture becomes unstable and the sample breaks completely. In studies of stress and strain in two dimensions, using a very large range of disorders [33], this turnover, from systems which are unstable from the onset to systems within which there is initially a complex coupled interaction between stress and disorder, is seen to occur for disorders approximately within the interval −2 < D < −3.
Presently, we include a single disorder for distributions with a tail toward strong elements, i.e., D = −4. The result obtained is ζ = 0.65(2), based on a straight-line fit to the average roughness values for L = 10, 13, 16, 20, 25, 32, and 40. This is shown in Figure 6, and the exponent is thus quite consistent with the results obtained for D > 0 distributions. With a given D expected to produce less roughness in three dimensions than in two dimensions, D = −4 here is probably close to this transient regime, where a controlled regime of micro-cracking appears. This could possibly account for the small difference obtained, i.e., ζ = 0.65 as opposed to ζ = 0.59-62 at high D > 0 disorders.

CONCLUSIONS
To summarize, a roughness exponent of approximately ζ 0.6 is obtained numerically for a discrete element model in the form of a cube lattice where forces are derived in analogy with elastic beams. Compared to other reported values in numerical works, this value is somewhat higher. One cannot really compare with results obtained in scalar models of fracture, such as the random fuse model, although the result is still in fair agreement with the various exponents obtained, i.e., 0.4 ζ 0.6, as reported in Batrouni and Hansen [9], Räisänen et al. [10], Räisänen et al. [11], Nukala et al. [12]. The Born model, which improves on the fuse model in that it includes elastic behavior, has yielded exponents in agreement with experiment, i.e., ζ 0.5 [20]. This model, as mentioned in the introduction is, however, not rotationally invariant.
The most relevant result with which to compare our current results are those of Nukala et al. [25], i.e., ζ 0.5. Our results are in fair agreement with this. A difference between the two studies is that in the former case periodic boundary conditions were used and also a slightly different fracture criterion. In two dimensions [24] we chose open boundary conditions to avoid having to address the problem of fitting a sine curve to the fracture surface. This because topologically the periodic lattice in two dimensions is a cylinder and the intersection of this with a plane is a sine curve. It may be difficult to judge whether a crack surface which appears to be a sine curve really is a sine curve or if this shape is something which appears as truly a part of the surface roughness. We chose to use open boundary conditions in three dimensions for similar reasons, and also generally to emulate a physical cube as it would appear in a fractured laboratory sample. To our knowledge, the significance of using different interaction formulae for combined loading is a subject which has not yet been investigated in studies using stochastic models. This is a subject which should be addressed in future studies.
Future studies should also address the issue of lattice topology. In two dimensions our calculations gave ζ = 0.86(3) [24] which is quite different from the two dimensional beam lattice result quoted in a study using a triangular lattice topology, i.e., ζ = 0.64 (2) [34]. The triangular topology in two dimensions is actually required for a Cosserat medium (such as the beam lattice models) to be able to reproduce the macroscopic properties of a real elastic solid. In three dimensions, the cubic lattice topology is inadequate in this respect.