The transition from subsonic to supersonic cracks

1. Introduction

One of the conclusions frequently drawn from the continuum theory of fracture mechanics is that cracks in tension cannot travel faster than the speed at which sound travels over a flat surface, the Rayleigh wave speed [1–3]. The reason for this assertion is that the motion of a crack requires energy to break bonds, but energy flux seems to become nonsensical for supersonic cracks. In particular, the rate at which energy flows into a crack tip per time is

where v is the crack speed, c_l and c_t are longitudinal and transverse wave speeds, respectively, , , μ is a Lamé constant, and K_I, the mode I dynamic stress intensity factor, is the coefficient of a universal singularity that develops outside of cracks that run as they are pulled symmetrically in tension from above and below. The denominator of this expression vanishes when the crack speed v reaches the Rayleigh wave speed, and for slightly higher velocities it becomes negative. Once the crack speed exceeds the transverse wave speed, the expression becomes imaginary. An expression saying that cracks moving above the Rayleigh wave speed need negative energy seems physically impossible, and an expression requiring imaginary energy seems even worse. One resolution of these problems is to conclude that cracks travelling faster than the Rayleigh wave speed in tension are not physically allowed.

However, there is evidence that such cracks exist after all. The first indications came from measurements of earthquakes [4]. These suggested computer simulations [5,6], and laboratory experiments [7] for cracks that move by sliding faces past each other, showing that cracks in shear (mode II) can move faster than the transverse wave speed c_t (and hence the Rayleigh wave speed as well) and reach speeds close to the longitudinal speed c_l. Dynamic fracture theory was extended to include these ‘intersonic’ cracks [2,8,9].

Cracks in rubber under tension were found to have a wedge-like tip suggestive of supersonic motion [10]. Additional experimental work confirmed that the cracks do travel faster than the transverse wave speed [11,12]. The experiments led to theoretical descriptions for supersonic cracks in tension [13,14]. Both explicit numerical solutions for atomic equations of motion and the corresponding analytical solutions show the supersonic cracks do in fact exist [15]. This last reference contains an extended discussion of how crack speed depends upon loading for the models studied in this manuscript.

Supersonic cracks were also observed for cracks with hyperelastic constitutive laws that cause stiffening near the tip [16]. As shown, for example, in the current article supersonic cracks exist equally well in materials without hyperelastic constitutive laws. Thus, we are not sure that hyperelasticity is actually needed to promote supersonic cracks.

Furthermore, the original objection from continuum theory about the impossibility of supersonic cracks has not directly been addressed. Doing so is a goal in this article. What we will explain in detail is the specific way that the continuum argument breaks down. In particular, we will show that, as a crack approaches the Rayleigh wave speed, the stress intensity factor K_I, upon which the continuum theory depends, stops being defined. The deformations around the crack tip take a completely different form from the one they have in the subsonic regime.

Our main technical tool to provide this discussion is the atomic theory of dynamic cracks. Analytical solutions at the atomic level for moving cracks were first found by Slepyan. His original calculations applied to out-of-plane (mode III) cracks in an infinite square lattice [17]. Slepyan and co-workers then generalized the solutions to in-plane (mode I) cracks in an infinite square lattice [18]. Many additional solutions due to Slepyan are found in [19].

Additional solutions for cracks in a finite strip formed from a triangular lattice were found by Marder & Gross [20]. This paper also considers linear and nonlinear instabilities of the solutions. Many additional observations on the experimental implications of the solutions in discrete lattices are contained in [21].

In all these analytical solutions, the central quantity obtained from the analytical methods is a relationship between the external load applied to make a crack move and the speed v at which the crack travels. In fact, the speed v of the crack enters the theory as an input parameter, and the calculations give the load as an output. For subsonic cracks, one can equivalently say that one obtains K_I(v), the stress intensity factor as a function of crack speed. For supersonic cracks, the end result is the external system strain e_yy as a function of crack speed [14].

The theory provides more than just this relation. In principle, the theory makes it possible to find the behaviour of every atom as a function of time for a crack moving in steady state. This information is relatively easy to obtain for cracks moving in anti-plane shear [20]. However, for two-dimensional cracks where atoms move freely in the plane, the expressions are so lengthy that to our knowledge no one has completed the process of extracting analytically the motion of every atom.

Thus, in this paper, we have two aims. The first is to present for the first time analytical solutions for every atom in a two-dimensional lattice where a crack moves in plane. The algebra is extremely lengthy, so we just present an outline of the work. Then, as an application of these exact solutions, we turn to the transition between subsonic and supersonic cracks. We examine the precise way that the stress intensity factor vanishes across the transition point and show what the cracks look like once they pass the Rayleigh wave speed. Because the methods produce solutions very rapidly, we are also able to discuss how these phenomena depend upon system size.

The structure of the article is as follows: in §2, we describe briefly the means we have used to solve for atomic motions around steady-state in-plane lattices. In §3, we examine the structure of the solutions in the vicinity of the Rayleigh wave speed, and explain how the stress intensity factor disappears. In §4, we make final remarks and conclude.

2. Explicit in-plane solution

(a) Lattice system

We carry out our computations in a triangular lattice of atoms with 2(N+1) rows. The motion of each atom is described by the displacement u_m,n from its equilibrium position, with m the column and n the row. The index m takes integer values while n takes values of the form (k+1/2), with k an integer ranging from −(N+1) to N.

The atoms in the lattice interact linearly with their nearest neighbours. The interaction is a function of both parallel and perpendicular displacements, with two spring constants k_∥ and k_⊥. We must make the restriction k_⊥=0 along the crack line, otherwise the mathematical formalism fails. We can make up for this by setting for the two rows along the crack. The bonds between atoms snap when they are displaced a distance 2u_f, and the interaction becomes zero.

We note that, in the long-wavelength limit, the motions of atoms in our lattice are described by isotropic continuum elasticity. In particular, longitudinal and transverse wave speeds are given, for unit mass and lattice spacing a, by

and In turn, this correspondence implies that crack motion in this lattice contains all results of dynamic linear elastic fracture mechanics as a special case.

To study crack motion, one must impose loading on the system. We do this by displacing the top row of atoms in our strip by a fixed vertical amount U_N. The crack consists in a separation between the middle rows of the lattice, where n=±1/2. We study the system when the crack moves at constant speed, in a steady state. Steady crack motion at velocity v in the continuum means that elastic fields are functions of x−vt. For cracks in a lattice, one cannot employ this definition. However, there is a related symmetry, which is that in a lattice of spacing a

This means that what an atom at some height does now, its neighbour to the right will repeat exactly at a time a/v later. Employing this relation lets us eliminate the m index and relate the components of different columns in the lattice to each other by a simple time symmetry. In particular, since we will be able to drop the index m and set about finding u_n(t).

(b) Analytical steps

The lattice in which we solve for crack motion is shown in figure 1. The displacement, for example, in the first direction is Δ₁=u_m−1,n+1−u_m,n. To find the force on an atom, define e_∥j and e_⊥j to be unit vectors in the unstretched lattice. The unit vectors e_∥j point from atom 0 to atoms j=1…6 in figure 1. Each of the unit vectors e_⊥j is perpendicular to the corresponding e_∥j. That is,

The total force on any given atom, except those along the crack line where some bonds may be broken, is and the motion of atoms, of unit mass, is given by The addition of the Stokes dissipation through the second term on the right-hand side is necessary to break time reversal invariance and tell the crack whether it should advance or retreat as time runs forward.

• Download figure
• Open in new tab
• Download PowerPoint

Define x to be the horizontal direction (the direction along which the crack moves) and y to be the vertical direction (along the width of the strip). The forces on atoms above the crack line, n>1/2, take the explicit component form

and

The number g_n helps keep track of the location on the triangular lattice. It is defined by

Much of the algebra needed to solve this system appears in detail in [20]. There is a missing factor of 12 in the equations of motion (VI.21, p. 48) that appear in [20], although the results reported in the body of the paper are correct. The corrected forces along the crack line, n=1/2, are

and where U(t) is a linear combination of u^x(t) and u^y(t) multiplying the Heaviside θ(t) functions

Fourier transforming the equations of motion gives

and where

These equations can be combined and rearranged to give finally

Here the function Q(ω) can be calculated explicitly from the Fourier transformed equations of motion for the lattice, but its explicit expression is so lengthy that it would serve no purpose to type it out here; we compute it from MAXIMA or Mathematica scripts, and for the purposes of rapid numerical evaluation find it in Fortran. In addition, Q₀=Q(0), and α is a small positive constant introduced when rewriting the boundary condition as U_N e^−α|t|. The Fourier transform of this function is easier to deal with than δ(ω). We let α tend to zero at the end of the calculation.

The Wiener–Hopf technique [22] lets us write

where Q⁻ is analytic in the lower half-plane, and Q⁺ is analytic in the upper half-plane. This decomposition is given by the explicit formula:

We can split court (2.16) into two pieces, one which is analytic only in the upper half-plane, and one which is analytic only in the lower half-plane:

The left- and right-hand sides of the above equation are analytic in opposite sections of the complex plane; therefore, they must both equal a polynomial P(ω). We must have P(ω)=0, otherwise U⁻(t) and U⁺(t) will include unphysical δ(t) terms. This lets us write

From the two functions U⁺(ω) and U⁻(ω), one can obtain the linear combination of horizontal and vertical motions defined in equation (2.12) from Finally, one can obtain the horizontal and vertical motions of any atom in the system from and The functions and , like Q(ω), are very lengthy to write down here, but can be obtained from essentially straightforward algebra resulting from the equations of motion. Again, we obtain explicit expressions using symbolic algebra programs, and implement the results for rapid evaluation in Fortran.

(c) Addition of the Kelvin dissipation

The solution we have described so far included the Stokes dissipation , which was present in equation (2.6). However, the Stokes dissipation is not physically realistic. It acts as though atoms are embedded in some sort of ether that slows them down according to the magnitude of their velocity in some arbitrarily chosen reference frame. A better form of dissipation from the physical point of view is the Kelvin dissipation, which produces a force opposing atomic motion according to the relative motion of adjacent atoms. To preserve generality, we keep both the Stokes and Kelvin dissipation in the equation of motion, which becomes almost, but not quite,

The reason this expression is not quite right is that when bonds break the time derivative produces unphysical discontinuities on the right-hand side. A solution of this problem was laid out in [23]. The solution is to define

Then the equations of motion on the crack line are modified to read

and

The functions W^± are defined just like U^± in equation (2.15). From this point, the analysis proceeds as before, finding a new function Q(ω) that is not appreciably more complicated because of the presence of β. From this function, one determines W using the Wiener–Hopf technique as

and Once U(ω) is in hand, the motion of every atom can again be determined. There is also a compact result relating the velocity v of a crack (which goes into the computation of Q) and the boundary extension U_N. The boundary extension can be expressed in terms of a dimensionless load Δ, Then, defining ω₀=1/β, Δ can be computed from Q(ω) directly as

(d) Numerical calculations

All atomic positions in time can be evaluated from rapid numerical operations involving nothing but algebra and fast Fourier transforms. First we sample the function Q(ω) over an interval of width 2ω_max and use fast Fourier transforms to compute U⁻(ω) and u_n(ω) from equation (2.18). Before taking the inverse Fourier transform, any singular behaviour of the form 1/(iω) must be subtracted from . The coefficient is determined from the asymptotic behaviour of and is given by

The subtracted function is then added back analytically after the inverse Fourier transform. This completes the numerical calculation of u_n(t). We can check the accuracy of the components knowing the exact asymptotic behaviour along the crack line,

Also, in the limit, the load is given by . This can be compared to direct numerical integration of equation (2.31). We can choose an ω_max that minimizes the difference between these two Δ values. Figure 2 shows components of the atomic position calculated in this manner. We choose the boundary condition U_N in such a manner that at the moment bonds break they have increased in length by u_f=1. Thus, we say that all displacements in the plots are measured in units of u_f. Keeping the breaking length fixed means in turn that all the computations correspond to physical systems with the same fracture energy. As the number of rows 2(N+1) increases, the extension U_N needed to bring the system to the point of fracture increases as . For supersonic cracks, the extensions become large multiples of u_f.

• Download figure
• Open in new tab
• Download PowerPoint

It is worth pausing to ask how figure 2 might be produced were one not using the Wiener–Hopf technique. It would be possible. It would require integrating the equations of motion for a crack in a system 1002 rows high. In order for atomic motions to reach steady state, the crack would need to run for a distance around 10 times the height of the system [24]. Although this could be sped up with a cutting and pasting procedure, accurate results would require a system 3000 columns long. Thus, one would have to run to steady state a system with around 3 million atoms; this would require a supercomputer. By contrast, producing figure 2 required just a bit over a second using a single processor on a laptop.

3. Behaviour of solutions

(a) Supersonic transition as a function of velocity

Having developed a tool that allows us very quickly to find the time history of atoms in the vicinity of a running crack, we now turn to the question with which we began. We ask what happens near a crack tip as the motion of the crack moves from subsonic to supersonic motion. The important wave speed is the Rayleigh wave speed c. For central forces (k_⊥=0), it is given explicitly in our model by

which is the root of the denominator in equation (1.1), and the limiting speed for cracks in continuum fraction mechanics.

Supersonic solutions in the lattice model look very different from subsonic solutions. Figure 3 shows how the vertical displacement of an atom on the crack line varies as the crack speed increases through the Rayleigh wave speed from v=0.9c to v=1.05c. We can see two distinct behaviours for subsonic and supersonic cracks. In both cases, the atom is nearly motionless until the crack approaches. After the crack passes, in subsonic solutions, the vertical displacement approaches the boundary condition U_N with small-amplitude, high-frequency oscillations that continue for a time of the order of 1/β as shown in figure 3a. These are phonons carrying energy left over after the bonds along the crack line have snapped [25].

• Download figure
• Open in new tab
• Download PowerPoint

In supersonic solutions, the vertical displacement also approaches the boundary displacement, but with large-amplitude, low-frequency oscillations as seen in figure 3h. The phonons of the subsonic solutions are no longer present in supersonic solutions. Interestingly, the subsonic phonons and supersonic large oscillations appear simultaneously in the subsonic solutions just below the Rayleigh wave speed, as we see in figure 3c–e.

The continuum solutions for in-plane fracture are graphed in [2]. Using the same notation as equation (1.1), the vertical displacement is

where R=4αβ−(1+β²)², and x=vt is the horizontal position of the crack tip as it moves in steady state. This same behaviour of u^y is present in the subsonic lattice solutions shown in figures 2b and 3a. We could use this square root displacement to extract a stress intensity factor K_I from the discrete system, as a complement to methods involving energy balance used in the past. However, as the crack speed increases towards the Rayleigh wave speed, the square root profile becomes increasingly indistinct, and it is completely lost at the Rayleigh speed and above.

As an illustration of the difference between subsonic and supersonic solutions, in figure 4 we plot the asymptotic behaviour of subsonic and supersonic solutions for atomic motion near the crack tip. For the subsonic solutions, the displacement rises as as expected. For supersonic solutions, it rises instead as .

• Download figure
• Open in new tab
• Download PowerPoint

(b) Dependence upon system size

The subsonic lattice solutions have a vertical displacement that matches the continuum result equation (3.2), superposed with small-amplitude, high-frequency oscillations that result from periodic bond breaking, and carry off all the energy flux to the crack tip not absorbed by the bond-breaking process.

The supersonic solutions feature large-amplitude, low-frequency oscillations. There is a simple explanation for these oscillations, which take place at what we call the supersonic block frequency ω_b. This frequency is what we find if the lattice oscillates vertically as a block, with the top row held fixed. We can solve for this frequency using the equations of motion along the crack line equations (2.10) and (2.11). These oscillations occur after the crack has passed, so we can eliminate the θ(t) functions multiplying U⁻(t). Fourier transforming these gives equations (2.13) and (2.14) with U⁻(ω) replaced by zero. The components u_3/2(ω) can be written as linear combinations of the u_1/2(ω) components, which gives us a system to solve. The determinant of this system for u_1/2(ω) must vanish, which lets us solve for ω_b, the lowest normal mode frequency. The coefficients in equations (2.22) and (2.23) are so lengthy that this procedure is extremely cumbersome and not very informative.

It is more useful to consider the lattice as a continuous block and use the equations of linear elasticity [1]. The top of the block is fixed while the bottom is free to oscillate. If the relaxed length of the block is L, the lowest normal mode frequency can be shown to be

In our lattice model, the result becomes . This expression is approximate to the extent that our lattice system is not actually a continuum, and because the waves are not infinite in horizontal wavelength. Figure 5 shows how the vertical displacement along the crack line for a fixed supersonic speed varies with the system size. As the system grows from 402 rows high to 3202 rows high, the oscillation frequency decreases in proportion. A fit to the calculated results gives ω_b∼3.1742/N. This result agrees with the result in equation (3.3) within 5%.

• Download figure
• Open in new tab
• Download PowerPoint

We have seen how the vertical displacement along the crack line changes as the crack speed becomes supersonic for a fixed system size (figure 3). The subsonic solutions just below the Rayleigh wave speed have both small-amplitude, high-frequency oscillations as well as large-amplitude, low-frequency oscillations. The low-frequency oscillations as seen in figure 3h begin to appear in figure 3c,d alongside the high-frequency ones, which disappear completely at the Rayleigh speed (figure 3f) and above.

In figure 6, we examine the vertical motion of atoms along the crack line for a fixed subsonic speed as the system size varies. For the small system with 202 rows (N=100), the solution shares the character of both subsonic and supersonic solutions. Phonon oscillations are visible, but at the same time there is a long-wavelength oscillation with the frequency of the supersonic block frequency ω_b. Now we increase the system size and monitor the location of the first long-wavelength peak. As one can see in figure 6 it slides to the right and diminishes in amplitude. The location of the peak is given approximately by t_b∼0.7150N, just as the period of the large-amplitude oscillations goes as T∝N in figure 5. The height of the peak decreases roughly as and approaches the boundary condition U_N asymptotically from above as seen in figure 6d. Thus for any given subsonic crack speed below the Rayleigh wave speed, there exists a system size sufficiently large to produce behaviour like the continuum solution equation (3.2), and this makes it possible to find a stress intensity factor.

• Download figure
• Open in new tab
• Download PowerPoint

4. Conclusion

We have presented the full analytical solution for atomic motions accompanying steady-state in-plane crack motion in a brittle triangular lattice. This solution allows rapid numerical evaluation of atomic motion for very large systems, which is necessary to make comparisons with results from continuum fracture mechanics.

Equation (1.1) has been interpreted to mean that cracks cannot propagate faster than the Rayleigh wave speed. However, supersonic solutions do exist in lattice systems. Subsonic solutions are characterized by small-amplitude, high-frequency oscillations in the time dependence of the vertical displacement of an atom along the crack line. The overall envelope of the displacement agrees with the continuum result in equation (3.2). Supersonic solutions are characterized by large-amplitude, low-frequency oscillations, which can be attributed to the lattice oscillating as a block after the crack has passed. We call this frequency the supersonic block frequency ω_b and show that its value can be approximated using linear elasticity. We have also shown that subsonic solutions close to the Rayleigh wave speed exhibit small-amplitude, high-frequency and large-amplitude, low-frequency oscillations at the same time. This mixing of the character of subsonic and supersonic solutions can be eliminated by increasing the system size N sufficiently.

There are three natural objections to what we have presented. The first is that systems 2000 atoms high are computationally challenging, but they are a far cry from macroscopic systems 10¹² atoms high. Thus, we have not really reached the macroscopic limit. The second is that we are drawing general conclusions about crack behaviour, but we have solved only a particular atomic model with very special and particularly simple interatomic forces. The third is that since the experiments of Schardin [26] it has been known that cracks do not reach the Rayleigh wave speed, let alone exceed it, so discussions of supersonic cracks are not relevant.

In response to the first objection, we have focused on scaling our solutions in such a manner that the large N behaviour becomes apparent. For example, in figure 6, we assert that in a system with N=10¹² the plot of vertical displacement versus time for an atom just above the crack tip would look indistinguishable from figure 6d. In figure 5, we assert that, in a system with N=10¹², the period of the large oscillations would continue to increase in proportion to N, and in order for the oscillations to remain visible the dissipation would need to decrease accordingly.

In response to the second objection, we invoke the universality of fracture mechanics. Although our atomic force laws are much simpler than any forces that genuinely would arise from the quantum mechanics of interactions, they capture the essence of brittle solids, and include all of dynamic linear elastic fracture mechanics as a special case in the long-wavelength limit. Thus, we believe that conclusions we draw about the continuum (large N) limit of our theory should have general validity. The details of what happens on small scales should be indicative of what would happen in real brittle systems, but would ultimately depend upon realistic microscopic details. For examples of what more realistic systems might look like, see [24].

In response to the third objection, we note that cracks in brittle isotropic materials are usually limited by dynamic instabilities to terminal velocities below the Rayleigh wave speed [25]. However, in some cases, the instabilities are suppressed and supersonic cracks become possible. Rubber provides one instance [11], and materials with highly anisotropic fracture energy might provide others. Supersonic cracks are not inevitable just because instabilities have been suppressed [27]: sufficient energy must be supplied to drive cracks above the Rayleigh wave speed. But they do become possible.

The calculations presented here are just the first application of this method. There are some additional questions still pending. Because the stress intensity factor stops being defined as the crack passes the Rayleigh wave speed, the conventional continuum view of energy flux stops applying. This observation provides an unsatisfying answer to the question of why the crack speed is not limited by energy flux. Once the crack travels faster than the Rayleigh wave speed, it is clear that it must extract the energy needed for bond breaking from a local environment of finite size.

Thus, we intend to investigate the energy flux described by (1.1) more closely from an atomic point of view, and ask whether the seemingly pathological behaviour it predicts above the Rayleigh wave speed can be given any physical interpretation, or whether there exists some other quantity that describes energy flux and can be employed in both the subsonic and supersonic regimes.

Funding statement

Support from the National Science Foundation Condensed Matter and Materials Theory Program (DMR 1002428) is gratefully acknowledged.

Footnotes

One contribution of 19 to a theme issue ‘Fracturing across the multi-scales of diverse materials’.