Mathematical models of supersonic and intersonic crack propagation in linear elastodynamics

This paper presents mathematical models of supersonic and intersonic crack propagation exhibiting Mach type of shock wave patterns that closely resemble the growing body of experimental and computational evidence reported in recent years. The models are developed in the form of weak discontinuous solutions of the equations of motion for isotropic linear elasticity in two dimensions. Instead of the classical second order elastodynamics equations in terms of the displacement field, equivalent first order equations in terms of the evolution of velocity and displacement gradient fields are used together with their associated jump conditions across solution discontinuities. The paper postulates supersonic and intersonic steady-state crack propagation solutions consisting of regions of constant deformation and velocity separated by pressure and shear shock waves converging at the crack tip and obtains the necessary requirements for their existence. It shows that such mathematical solutions exist for significant ranges of material properties both in plane stress and plane strain. Both mode I and mode II fracture configurations are considered. In line with the linear elasticity theory used, the solutions obtained satisfy exact energy conservation, which implies that strain energy in the unfractured material is converted in its entirety into kinetic energy as the crack propagates. This neglects dissipation phenomena both in the material and in the creation of the new crack surface. This leads to the conclusion that fast crack propagation beyond the classical limit of the Rayleigh wave speed is a phenomenon dominated by the transfer of strain energy into kinetic energy rather than by the transfer into surface energy, which is the basis of Griffiths theory.


Introduction
Dynamic fracture mechanics is a well establish field of research with significant relevance in many areas of engineering (Baker 1960;Craggs 1960;Freund 1972;Itou 1978;Ravi-Chandar and Knauss 1984a, b;Freund 1990;Xu and Needleman 1994;Washabaugh and Knauss 1994;Gao 1996;Sharon et al. 1996;Willis and Movchan 1997;Morrissey and Rice 1998;Hauch et al. 1999;Broberg 1999;Aranson et al. 2000;Cramer et al. 2000;Belytschko et al. 2003;Cox et al. 2005;Gross and Seelig 2011). It has been generally accepted that in an elastic medium cracks can only grow and propagate at speeds below the speed of Rayleigh waves (Fre-und 1990;Xu and Needleman 1994;Washabaugh and Knauss 1994;Gao 1996;Sharon et al. 1996; Willis and Movchan 1997;Morrissey and Rice 1998;Hauch et al. 1999;Broberg 1999). Beyond this speed linear elastic solutions can be shown to lead to a vanishing energy release rate (at least in Mode I fracture) which is incompatible with the energy required to create the new crack surface. However, recent results both experimental and computational have shown crack propagation at significantly higher speeds both in the intersonic range and even beyond the speed of sound (Rosakis et al. 1999;Needleman 1999;Abraham and Gao 2000;Huang and Gao 2001;Gao et al. 2001;Rosakis 2002;Abraham et al. 2002;Guo et al. 2003;Buehler et al. 2003;Hao et al. 2004). Moreover, some of the published results exhibit features akin to shock waves of supersonic flow where solution discontinuities are seen to emerge from the crack tip as it travels along the material (Petersan et al. 2004;Willmott and Field 2006;Radi and Loret 2008;Bizzarri et al. 2010;Schubnel et al. 2011;Barras et al. 2018;Yue et al. 2019;Mai et al. 2020). No analytical models have been reported that predict such behaviour in linear elasticity. This paper will provide mathematical models with solutions that closely resemble some of these results. These are not inconsistent with existing analytical theories in so far as they are based on full transfer from strain into kinetic energy and therefore neglect the energy required to create the new surface.
Classical mathematical analysis of dynamic crack propagation has been based on the solution of the governing second order equation for the displacement field as a function of time and position (Freund 1972;Burridge 1973;Freund 1990;Andrews 1976;Freund 1979;Burgers and Freund 1980;Broberg 1989Broberg , 1994Andrews 1994;Willis 2003, 2008). Freund (1972) obtained a fundamental solution for Mode I type fracture in an infinite linear elastic domain by considering the solution of the second order equation for a half-plane in conjunction with cohesive type tractions on the fracture plane, leading to rupture velocities limited by the Rayleigh wave speed. Burridge (1973), using similar methods in the case of Mode II type fracture, showed that in the case of friction type forces (lacking cohesion), it is possible to obtain rupture velocities of the value of the dilatational wave velocity. Andrews (1976), resorting to a finite difference type method, reconciled Freund's and Burridge's solutions by considering a slip-weakening type model which combined friction and cohesion. This was further confirmed by Freund in Freund (1979), Burgers and Freund (1980) and Broberg in Broberg (1989) and Broberg (1994), the latter for the special case of a Mode II rupture velocity equal to √ 2 the shear wave speed. Andrews in Andrews (1994) extended his previous work to the case of mixed-mode shear cracks, albeit this time employing analytic solutions and the boundary element method, obtaining similar conclusions, that is, mixed-mode ruptures can propagate faster than the Rayleigh speed. In Obrezanova and Willis (2003), Obrezanova and Willis (2008), the authors study the stability of intersonic shear cracks and provide stress fields in the vicinity of the crack, first in plane strain and subsequently in three dimensions.
Numerical methods and solutions enumerated above are appropriate for smooth solutions but side-step possible solutions where the first derivatives of the displacement field are discontinuous across moving surfaces. It is now well-accepted that this discontinuous type solutions exist and they exhibit complex Mach cone type profiles for velocities and strains/stresses. Indeed, numerous references have reported these type of solutions through either experiments or atomistic simulations (Abraham and Gao 2000;Gao et al. 2001;Rosakis 2002;Abraham et al. 2002;Buehler et al. 2003;Hao et al. 2004;Petersan et al. 2004;Schubnel et al. 2011;Yue et al. 2019;Abraham 2001). These are known as weak solutions and satisfy certain jump conditions across the moving discontinuities, known as shocks (Eringen and Suhubi 1975;Gurtin et al. 2010;Bonet et al. 2021). Motivated by improving numerical simulations, significant research has been recently devoted to re-expressing the fundamental equations of solid dynamics as a set of conservation laws for the linear momentum and the deformation (Lee et al. 2013(Lee et al. , 2014Aguirre et al. 2014;Bonet et al. 2015;Aguirre et al. 2015;Haider et al. 2017). These conservation laws lead to a system of first order differential equations for the velocity field and the deformation together with associated jump conditions across moving shocks. These equations are well known in the classical literature but typically only used as a step towards a final second order equation in terms of displacements. Using them directly allows for weak solutions to be obtained that could otherwise be by-passed. First order conservation laws with moving discontinuities are commonly used in the solution of fluid dynamic problems (LeVeque 1992;Toro 1999;Trangenstein 2009) where solu-tion patterns incorporating Mach shock waves are frequently encountered.
The aim of the paper is to provide mathematical solutions of steady-state crack propagation in the context of linear elastic elastodynamics, where the governing equations are expressed as a set of first order linear hyperbolic conservation laws with accompanying jump conditions. The resulting solutions are two-dimensional extensions of well known so-called "local" solutions in one dimension (Bažant and Belytschko 1985). It is not the premise of the paper to claim that these solutions can provide an exact description of the physical reality as they are constrained by the conservative nature of linear elasticity; whereas, clearly, fracture is an inherently dissipative process. Energy losses through the creation of new crack surfaces and inelastic dissipative effects near the tip are ignored by the purely elastic solutions proposed. Instead, the aim of the paper is to show that Mach type of shock wave patterns as well as intersonic and supersonic crack propagation can be explained in linear elasticity under the assumption that surface energy is negligible in relation to the available strain energy in the material. This represents the opposite end of the energy balance spectrum to the classical Griffiths theory, where all the available strain energy is assumed to be used to create the new surfaces. It is clear that real physical phenomena will be found somewhere in between these two extremes. The contribution of this paper is therefore to suggest that fast crack propagation leading to shock wave formation is governed by the transfer of strain energy into kinetic energy rather than into surface energy. An additional intended contribution of the paper is to provide the outcome of the proposed solutions as initial and boundary conditions for computational models where dissipative effects can then be incorporated with relative simplicity. This will be the subject of future research.
The paper is organised as follows. As a way of motivation, Section 2 provides a mathematical description of the simple 1-D fracture of a bar in terms of first order differential equations describing the evolution of the velocity and strain with associated jump conditions in linear elasticity. This leads to a simple so-called "local" solution in which strain energy is transferred into kinetic energy in its entirety and the physical implications of this transfer are discussed. Section 3 extends the first order linear elastodynamics equations to 2-D in terms of velocity vector and displacement gradient.
The speed of pressure (or dilatational) and shear (or distortional) shock waves are derived from the jump conditions and linear elasticity constitutive equations. Section 4 postulates a two-shock solution for supersonic crack propagation in mode I fracture and determines the constitutive conditions for this solution to exist. As in the 1-D model exact energy preservation is shown to hold through the transfer of elastic energy into kinetic energy. Section 5 examines supersonic mode II crack propagation and obtains very similar solutions to mode I, albeit with more severe restrictions in relation to the range of material properties for which the model is valid. Section 6 presents an intersonic solution where the pressure shock travels ahead of the crack tip and the shear shock. Section 7 provides some concluding remarks and recommendations for further work.

Motivation: the local one-dimensional fracture solution
In order to motivate the more complex two-dimensional solutions developed in the following sections, consider first the simplest possible fracture problem consisting of a one-dimensional bar under uniaxial tension as shown in Fig. 1. For simplicity, it will be assumed that the bar is initially at rest, linearly elastic and under the action of a constant external traction T as shown in the Figure. At time t = 0 fracture takes place at the centre of the bar and the resulting motion is governed by a set of first order hyperbolic equations given as, where ρ denotes the density, v(x, t) the velocity field, σ (ε) = Eε is the stress, ε(x, t) the strain field and E the Young's modulus of elasticity. Note that the above set of equations represent a set of first order linear hyperbolic equations from which the standard second order wave equation for the displacement u(x, t) can be easily derived by noting that v = ∂u/∂t and ε = ∂u/∂ x. However, in order to establish weak solutions, that is solutions that contain discontinuities, it is preferable to consider the first order equations directly together with where the notation a = a + − a − has been used to indicate the jump of a variable across a shock or discontinuity travelling at speed c along the bar. Noting that σ = E ε and substituting (2b) into (2a) gives the speed of propagation of a discontinuity as, Note that this corresponds to the speed of the sound wave in the bar which can be obtained considering smooth solutions of classical wave propagation theory.
Consider now the material fracturing at t = x = 0 as shown in Fig. 1. This creates an initial condition where the stress at the fracture point vanishes and two release waves travel in opposite directions along each half of the bar as shown in the Fig. 1. The above jump conditions can be applied between the values of the variables ahead and behind the shock in either the left or right halves of the broken bar to give, and therefore, an instantaneous release velocity at either side of the fracture point is generated at time t = 0 with a value given by, v = σ ρc .
This represents a complete transfer of elastic energy at the point of fracture into kinetic energy as shown by the simple identity, Essentially, fracture is modelled in the above equations as a sudden removal of the equal and opposing internal traction forces acting on the section where fracture takes place. The above model, often referred to as the local solution (Bažant and Belytschko 1985), neglects the amount of energy required to generate the new surfaces at the point of fracture and therefore energy is exactly transferred from elastic to kinetic energy. From the point of view of hyperbolic differential equations, this results in a Riemann problem in each half of the broken bar with a simple analytical solution (LeVeque 1992;Toro 1999;Trangenstein 2009). The boundary fluxes at the point of fracture, namely stress and velocity, must be compatible with the jump conditions. In linear elasticity this compatibility leads to exact energy balance since the shock wave travels at the speed of sound. It is the end velocity given in equation (5) that would be responsible for a slender wire curling in a whipping like fashion as it snaps under tension or a rubber band snapping back after fracture. Note that from a mathematical point of view the problem is well posed in elastodynamics without the need to introduce dissipative or length effects, provided that the first order conservation equations and the jump conditions are used (Bažant and Belytschko 1985).
The energy balance in the solution described above represents the opposite end of the spectrum to the classical Griffiths theory where all the available strain energy is transferred to surface energy. Arguably, all brittle fracture process would occupy a point somewhere between these two extremes, where some of the strain energy is transferred into kinetic energy and some into surface energy and some dissipated as heat. Of course, if the aim is to assess the possible onset of crack growth in a metal sheet, it is prudent to consider the energy balance as stated in classical Griffiths theory, leading to a quasi-static fracture problem. Alternatively, modelling the curling motion of a snapping cable (Tjavaras et al. 1998) should be approached from the other end. For instance, a high strength steel cable would snap at a strain energy level in the order of magnitude of 1 MJ/m 3 , whereas the order of magnitude of the surface energy is 1J/m 2 . Hence, the energy required to generate the new surfaces is contained in a sliver of cable of thickness in the region of 1μm, which is negligible in comparison to the length of the cable.
From a mathematical point of view, dissipative effects to account for the surface energy can be introduced in the above hyperbolic equations by including conductive fluxes proportional to the gradient of the variables. This can be done through a dependency of stress on the gradient of the velocity, leading to a visco-elastic model (Schreyer and Chen 1986;Simo and Ju 1987;Sluys and de Borst 1992;Etse and Willam 1999;Bažant and Jirásek 2002), or a dependency on the gradient of ε which results in a classical weakly non-local model (Bažant and Jirásek 2002;Belytschko and Lasry 1989;De Borst and Muehlhaus 1992;Fleck and Hutchinson 1993;De Borst 2001). From a computational point of view, many alternative methodologies have been derived, ranging from smeared crack theory (Pijaudier-Cabot and Bažant 1987;Jirásek and Zimmermann 1998), strong discontinuity approaches (Simo et al. 1993;Chen and Sulsky 1995), cohesive zone models (Barenblatt 1962;Falk et al. 2001;Camacho and Ortiz 1996), extended finite element methods (Moës and Belytschko 2002;Dolbow et al. 2001), configurational force methods (Miehe and Gürses 2007), phase field models (Miehe et al. 2010;Hesch and Weinberg 2014;Hesch et al. 2017) and others.
The sections that follow will essentially extend the above local solution to 2-dimensional mode I and mode II fracture configurations.

Linear elastodynamics in 2-D
The aim of this section is to review the equations of linear elasticity in the form of first order hyperbolic equations with associated jump conditions, so that conceptually similar solutions to the above one-dimensional model can be derived in the case of 2-dimensional isotropic elasticity.

First order hyperbolic governing equations
The general dynamic equilibrium in linear elasticity is established by the momentum conservation equation which in the absence of external body forces is given by Eringen and Suhubi (1975), Gurtin et al. (2010), Bonet et al. (2021), where div indicates the divergence operator (the contraction of the gradient by the last index) and σ denotes the Cauchy stress tensor, which in linear elasticity can be related to the displacement gradient G = ∇u via a symmetric fourth order elasticity tensor C as, In this expression λ, μ denote the Lamé coefficients, I the second order identity tensor, I the fourth order identity tensor andĨ the fourth order tensor that maps second order tensors to their transpose, that isĨ : G = G T . Equation (8) is more commonly written in terms of the symmetric strain tensor ε = (G+G T )/2 rather than the displacement gradient G but the outcome will be the same when equation (8) b is used as the constitutive relationship.
Note that plane stress can be obtained by simply replacing λ in equation (8) b by its equivalent or effective valueλ defined as, No distinction will be made in what follows between plane strain and plane stress on the understanding that λ is replaced byλ for the latter.
A first order hyperbolic equation for the evolution of G similar to equation (1b) is given by, This equation simply states that the rate of change of G is given by the gradient of the velocities. Equation (10) is the small strain counterpart of the evolution equation for the deformation gradient described in references Lee et al. 2013Lee et al. , 2014Aguirre et al. 2014;Bonet et al. 2015;Aguirre et al. 2015;Haider et al. 2017).
Equations (7) and (10) represent a set of first order linear hyperbolic equations for which analytical solutions can be obtained under appropriate initial and boundary conditions. However, the type of solutions sought are weak solutions containing moving discontinuities, which require the formulation of the jump conditions associated to the conservation laws (7) and (10). These are derived in the following sections.

Jump conditions and wave speeds
Consider a shock in the solution, that is a discontinuity surface with unit normal n travelling inside the solid with normal speed c as shown in Fig. 2. Variables ahead of the shock are denoted with a "+" sign, whereas the value of variables behind the shock is identified with the superscript "-".
In the presence of discontinuous solutions, the conservation equations (7) and (10) have associated jump conditions given, respectively, by Eringen and Suhubi (1975), Gurtin et al. (2010), Bonet et al. (2021), Bonet et al. (2015), Using the constitutive relationship between the stress tensor and the displacement gradient given by equations (8) and substituting equation (11b) into equation (11a) gives, after simple algebra, an eigenvalue equation in terms of the acoustic tensor A as, There are two solutions to this eigenproblem, corresponding to pressure shock waves (also known as dilatation or longitudinal) and shear waves (Eringen and Suhubi 1975;Gurtin et al. 2010;Bonet et al. 2021Bonet et al. , 2015. In the first instance the jump in velocity is proportional to the normal vector and the shock travels at the speed of sound as, v = v n n; The second solution is obtained by taking velocity jumps parallel to the surface tangent ν leading to shear shocks travelling at the speed of shear wave as, The ratio between the pressure and shear shock speeds will be a useful material parameter in the derivations that follow and can be variously expressed in terms of the Lamé coefficients or Poisson's ratio ν as, Plane stress: Note that in the incompressible limit this ratio tends to infinity in the case of plane strain but in plane stress the speed of the pressure wave cannot exceed twice the speed of the shear wave. In fact, in plane stress the pressure wave should be more correctly described as normal or longitudinal wave as it involves changes in the thickness and therefore transversal out of plane distortion.

Mode I crack propagation
This section describes a mathematical solution for a propagating crack under mode I fracture consisting of a pair of shock waves, a pressure and a shear wave, propagating from the crack tip. Prior to the arrival of the crack the material is in a state of uniform stress given by a vertical traction σ = σ 22 e 2 ⊗ e 2 as shown in Fig. 3.
Inspired by the one-dimensional model described in Sect. 2, a solution is postulated whereby two shock waves, a shear and a pressure wave, separate three regions of constant stress, strain and velocity as shown in Fig. 3, where due to symmetry only the upper half of the continuum is depicted. Notice that a similar Mach cone type solution was experimentally reported in Fig. 2 in Reference (Petersan et al. 2004) and in Figs. 4 and 5a in Reference (Yue et al. 2019). The regions I, II and III are formally defined as: • Region I: unfractured body, satisfying: • Region II: transition region with stress, strain and velocity to be determined by the jump conditions; • Region III: unloaded region with stress tensor satisfying no traction on the crack surface: These three regions are separated by straight pressure and shear shocks originating from the crack tip and travelling at angles α p and α s , respectively. At the crack tip, the intersection of the shear and pres-sure waves with the crack imply that the speed of the crack is related to the pressure and shear wave speeds as, Note that the existence of a solution as described above entails supersonic crack propagation as equation (18) a implies that v c ≥ c p ≥ c s . Intersonic solutions can be found if the pressure shock wave is allowed to travel ahead of the crack tip. These solutions will be explored in a later section.

Velocity and stress relationships
The constant state of stresses and velocities in the three regions described above implies that the differential equations (7) and (10) are trivially satisfied. The key to determine the validity of the postulated solution is the satisfaction of the jump conditions at each shock wave and appropriate boundary conditions at the crack surface. Consider first the enforcement of the momentum jump condition (11a) at the pressure shock wave as, Across the pressure shock the velocity jumps, and therefore the traction jumps, will be parallel to the nor- where the shock condition, equation (19), implies that the normal velocity and normal traction jumps, v n and t n respectively, are related as, Consider next the enforcement of the displacement gradient jump condition (11b) across the pressure shock as, Substituting for the jump in velocity from equations (20) a and (21), gives the displacement gradient in the intermediate region as, Multiplying the above equation by the elasticity tensor gives the stress tensor at the intermediate region as, The same derivations can now be repeated across the shear shock where the velocity and traction jumps are parallel to the tangential unit vector as, The jump condition now gives a relationship between the jump in tangential velocity and shear stress as, The displacement gradient jump condition across the shear shock combined with the above relationship between jumps in tangential velocity and shear stress gives, Finally, multiplying by the elasticity tensor and noting that μ = ρc 2 s yields an expression for the stress tensor in region III as, Combining equations across both shocks enables direct relationships to be established between the velocity and stresses in region III and those in region I. For instance, using equations (20) a , (21), (25) a and (26) and noting that the velocity in region I vanishes gives an expression for v (I I I ) as, Similarly, combining equations (24) and (28) gives the stress in region III as, Note that the displacement field u(x, y, t) in the proposed solution can be obtained by integrating in time the velocity field at any given point as the shock fronts go past. Since the velocity and gradient fields are constant in each region, the resulting displacement is inevitably linear in t, x and y. Therefore, the classical second order equation in u(x, y, t) is trivially satisfied. Essentially, weak solutions of the form proposed here require the use of the first order equations and their associated jump conditions and are missed by the second order classical wave equations.

General solution
In order to explore the existence of solutions in the form postulated in the sections above, the stress and velocity in region III need to satisfy appropriate boundary conditions on the crack surface. In particular, the traction on the crack surface must vanish, giving two scalar conditions as, e 1 · σ (I I I ) e 2 = 0; e 2 · σ (I I I ) e 2 = 0.
In order to derive a third condition, it is necessary to determine whether the crack propagates symmetrically or from a free surface as shown in Fig. 4. This section will consider only the symmetric case, in which case a third condition can be found by imposing that the velocity v (I I I ) in region III should be vertical, that is, v (I I I ) · e 1 = 0.
Substituting for the stress and velocity in region III from equations (29) and (30) into equations (31) and (32) gives a set of three equations with three unknowns, namely t n , τ and either the angle α p or α s since they are connected to the speed of crack propagation (33a) ν p,s = − cos α p,s e 1 + sin α p,s e 2 . (33b) Multiplying equation (29) by e 1 and employing the above trigonometric relationships together with the crack propagation condition (18) b gives, after simple algebra, a relationship between the jumps in shear stress and normal traction as, Enforcing condition (31) a and noting again that μ = ρc 2 s , gives after simple trigonometric algebra, c 2 s t n sin 2α p = c 2 p τ cos 2α s , which, when combined with condition (34) and taking into account the angle relationship provided by equation (18) b gives a further geometric relationship between the angles of the shear and pressure shocks as, Adding the square of this equation to the square of equation (18) This leads to a simple equation for the shear shock angle as, where only the positive solution has been considered as the negative value of the cosine leads to the geometrically obvious symmetric solution below the crack. The above equation for α s has a valid solution provided that c p ≤ 2c s or λ ≤ 2μ , that is the speed of the pressure wave is not higher than twice the speed of the shear wave. For the case of plane stress, where λ should be replaced byλ as given in equation (10), it is easy to show that this is the case for any value of the Poisson's ratio in the physically meaningful range. For the plane strain case, valid solutions will only exist for Poisson's ratios ν ≤ 1/3. Beyond this value the solution postulated is no longer possible in isotropic linear elasticity. Finally, combining equation (18) b with equation (38) it is easy to show that the angle of the pressure shock is, These expressions can now be used to determine the complete velocity and stress fields in the unloaded and transition regions. For instance, the velocity in region III, which by construction will only have a vertical component, can be derived using equations (42)  Note that this gives a velocity field identical to the one derived in the one-dimensional case in equation (5).
The state of stress in region III is by construction limited to a horizontal traction as, The above expression can also be re-written in terms of the Poisson's ratio for the plane strain and plane stress cases respectively as, Observe that by construction the above state of stress and velocity would be compatible with a crack propa-gating symmetrically in the left and right directions as shown in Fig. 4a. As the shear shock lines reach each other on the central symmetry axis reflections would take place into region II which would interact with the pressure shock in a complex manner as shown in Fig. 6. The analysis of this interaction and reflection is beyond the scope of the simple solution proposed here.
For completeness, it is useful to derive the velocity and stress fields in the transition region II. This is easily achieved using Eqs. (20) a , (21), (24) and (42) The solution described above represents a mode I crack propagating symmetrically at speeds equal or higher than the speed of sounds in a manner consistent with the equations of linear isotropic elasticity. In the plane strain case, solutions in the simple form postulated here do not exist beyond ν = 1/3 at which point the speed of propagation becomes infinite. The same effect takes place in plane stress at the incompressible limit ν = 1/2. For a Poisson's ratio different to zero, viable solutions where the horizontal stress in region III vanishes, which would correspond to cracks propagating from a free surface as shown in Fig. 4b have not yet been found with shocks emanating jointly from the crack tip, ie., in the supersonic range. However, inter- Fig. 6 Propagation of a symmetric crack showing the shock reflections as dashed lines at the symmetry axis sonic solutions where the pressure shock moves away from the crack tip exist and will be described in Section 6.

Particular solutions
It is interesting to consider particular cases of the general solution derived above. For instance the case ν = 0, or equivalently λ = 0, leads to an interesting particular solution where the angle of the shear shock is 45 • and the pressure shock travels at 90 • as shown in Fig. 7 (refer also to Fig. 5b). In this case the crack propagates at sonic speed (refer to Fig. 5a) and region III is fully unstressed. The pressure wave acting as a compressive shock wave first transforms the initial state of uniaxial vertical stress into a state of pure shear stress given as, The shear shock then eliminates this state of shear stress and leaves the material fully unstressed in region III as shown in Fig. 7. The absence of horizontal stress in region III makes this solution valid for both the symmetric crack propagation case as well as for a crack propagating from a free surface.
As a further example of a particular case obtained from the equations developed in the previous sections, consider the case where λ = μ and therefore c p /c s = √ 3. This corresponds to a Poisson's ration of 0.25 in the case of plane strain or 1/3 in plane stress. The angles of the shear and pressure shock waves are 30 • and 60 • respectively and the crack propagates at a speed of, The state of stresses in regions III and II is given for this particular case as, where matrix notation has been used for simplicity in equation (50) b .

Energy balance
In the one-dimensional case presented in Sect. 2, it is clear that strain energy before fracture is converted into kinetic energy fully respecting energy conservation. It is also possible to show that this is the case in the proposed 2-D model, given that it has been developed using the equations of linear elastodynamics which do not incorporate any dissipative effects. In particular, it will be shown in this section that the elastic energy in region I is transformed into a smaller amount of strain energy in region III associated with the horizontal stress that emerges due to Poisson's effect and the difference becomes kinetic energy in region III. Consider the fixed control area shown in Fig. 8. It is simple to show that the energy flux across the boundaries of this region vanish so that the total energy inside The energy flow in region I is trivially zero as the velocity vector vanishes. In region III the energy flow across (I I I ) also vanishes as the traction vector is horizontal whereas the velocity vector is vertical. Showing that the energy flow across the section of the control boundary in region II also vanishes requires a more detailed evaluation using equations (47) to give, Consequently, the energy in the control area must be conserved as the crack advances provided that there are no dissipative effects as the crack opens, which is the case in linear elastodynamics. It is easy to see that the area of region II inside the control volume remains constant as the crack progresses whereas the areas of regions I and III decrease and increase respectively in equal measure. Hence the energy per unit volume in regions I and III must coincide to ensure conservation of energy. In order to show that this is the case, note first that the strain energy is obtained in terms of the stress tensor and the 2-dimensional inverse elasticity tensor as, Using this equation at regions I and II gives after simple algebra, It is a simple exercise to show that the difference between these two energy expressions coincides with the kinetic energy per unit volume in region III as, Hence, in the proposed linear elastic solution strain energy is not dissipated or transformed into surface energy but exactly transformed into kinetic energy as the crack progresses. This represents the opposite end of the spectrum to the Griffiths energy balance where strain energy is transformed into surface energy.

Mode II crack propagation
The possibility of obtaining similar solutions to those derived above for the case of mode II fracture is discussed in this section. The unfractured stress field will be assumed to be a pure state of shear where the stress The release conditions in the fractured region III will now imply a horizontal velocity together with a vanishing traction vector at the crack surface. The resulting configuration is shown in Fig. 9 and leads to the following set of conditions, e 2 · σ (I I I ) e 2 = 0; e 1 · σ (I I I ) e 2 = 0; v (I I I ) · e 2 = 0.
Note that in this case there is no vertical symmetry so the solution below the crack will need to be determined explicitly and independently but will need to satisfy appropriate skew-symmetry conditions. Substituting for the velocity and stress fields into equations (57) a,c from equations (29) and (30) gives, after simple trigonometric algebra and use of (33), a set of equations for the normal and shear stress jumps and the shock angles as, t n cos α p + τ sin α p = 0; (58a) t n cos 2α s + τ sin 2α s = 0.
These equations lead to a relationship between the shock angles as, which in the 0 • to 360 • interval has the solutions, Notice that similar angular relationships can be observed in Figs. 2 and 4 in Reference (Abraham and Gao 2000), Figs. 5 and 8 in Reference (Abraham 2001), Fig. 1 in Reference , Fig. 26 in Reference (Rosakis 2002), Fig. 4 in Abraham et al. (2002), Fig. 4b in Reference (Buehler et al. 2003) and Fig. 20c in Reference (Hao et al. 2004), where a Mach cone type solution is obtained using an atomistic simulation. In their case, for an equivalent Poisson's ratio of ν = 0.3, angles close to α s = 30 • and α p = 60 • are shown, which match those depicted in Fig. 5b. A similar Mach cone type solution can also be found in Fig. 8 in Reference (Schubnel et al. 2011). Combining these angle relationships with the crack propagation condition (18) b gives the same angles as those obtained for mode I fracture plus a second solution which represents symmetric angles respect the horizontal axis as shown in Fig. 9. In particular, In order to complete the solution, it is necessary to determine the normal and shear stress jump in terms of the shear stress in region I. This can be done by using condition (57) b , which after some algebra leads to, Using these expressions, the values of the horizontal velocity and stress in the unloaded top and bottom regions emerge as, = 2τ (I ) 1 ν = τ for plane stress, respectively. In both these Eqs. (64) and (65), the plus sign corresponds to the solution above the crack and the minus to the solution below the crack. Note that these equations become singular as λ → 0 which corresponds to ν → 0 (refer to Eqs. (66), (67) and Fig. 10). In fact, well before reaching this value, the above state of stresses would imply a state of shear at 45 • greater than the original shear in region I so the above solution has a limited range of physical validity. For completeness, the velocity and stress in the transition region II can be found by using Eqs. (20) a , (21), (24) and (63)

Energy balance
A similar energy balance to that presented in Sect. 4.6 can be established for mode II fracture. In this case, the energy balance equation is slightly more cumbersome as the energy flow is non-zero in a part of the boundary of the control volume. The energy flow in region I is zero as the velocity vector vanishes. In region II, the energy flow across (I I ) is also zero as In region III, the energy flow is zero across the horizontal sides of (I I I ) as the traction vector across them vanishes (i.e. σ (I I ) e 2 = 0), reducing to As the area of region II inside the control volume remains constant as the crack progresses, the decrease of the area of region I is balanced by the increase of area of region III. Thus, referring to Fig. 11, the time rate of the energy in the control area must be equal to the energy flux as the crack advance, which can be expressed as It is straightforward to demonstrate that above energy balance equation is fulfilled for the stress and velocity states defined previously.

Intersonic mode I solutions
The two shock model described above can be used to derive intersonic solutions by allowing the pressure shock to move ahead of a shear shock which is still attached to the crack tip. Similar intersonic shock distribution has been reported in Fig. 21 Reference (Rosakis 2002) and Fig. 4a in Reference (Buehler et al. 2003). This implies that the relationship between α p and α s given in Eq. (18) a is no longer valid, leading to an additional unknown in the system of equations. The resulting configuration for mode I fracture is shown in Fig. 12 and generalizes the particular solution depicted in Fig. 7 to cases where the Poisson's ratio is greater than zero. It will be shown below that the angles of the pressure and shear shocks will be 90 • and 45 • respectively as shown in the figure.
The additional equation required to solve the system in the absence of (18) a can be obtained by noting that the symmetry condition along the horizontal axis implies that the shear stress in region II must vanish. With the help of Eq. (24), this gives the following condition for α p as, 0 = σ (I I) : (e 1 ⊗ e 2 ) = σ (I ) : (e 1 ⊗ e 2 ) − t n ρc 2 p λI + 2μ(n p ⊗ n p ) : (e 1 ⊗ e 2 ) Equation (72) implies that α p = 90 • as shown in Fig. 12. The angle of the shear shock can be derived by enforcing the same condition in region III to give,  Thus leading to α s = 45 • . The resulting crack propagation speed is therefore given by, which falls in the intersonic range except for the case where the Poisson's ratio vanishes in which case the crack propagates at the speed of sound. It is now possible to find solutions for both the symmetric propagation case, where v (I I I ) · e 1 = 0, and the free surface case where σ (I I I ) = 0 . It is easy to show that in the former case the shear and normal stress jumps are related by, whereas in the latter the jumps are related by, In both cases, the condition of vanishing vertical stress in region III, gives a relationship between τ , t n and σ (I ) 22 as, This equation combined with either Eq. (75) for the symmetric case or Eq. (76) for the case where the crack propagates from a free surface will fully determine the states of stress and velocity in regions II and III. Specif-ically, for the symmetric case, these are v (I I ) = t n ρc p e 1 ; σ (I I ) = − t n e 1 ⊗ e 1 + τ e 2 ⊗ e 2 ; (78a) v (I I I ) = t n ρc p e 2 ; and for the free surface case, these are v (I I ) = t n ρc p e 1 ; Regarding the energy balance, a similar analysis can be carried out to that presented in previous sections for supersonic modes I and II of fracture. In this case, it is easy to notice that the energy flow across the boundaries of the control volume is zero for the three regions I, II and III. However in this scenario, due to the intersonic nature of the fracture mode analysed, the area of region II is not preserved and the energy balance becomes which can be easily verified by substitution.

Conclusions
The paper has found novel local solutions to mode I and mode II fracture in isotropic linear elasticity that describe supersonic and intersonic crack propagation showing remarkable similarities to experimental and computational models previously reported in the literature Rosakis 2002;Abraham et al. 2002;Guo et al. 2003;Buehler et al. 2003;Hao et al. 2004;Petersan et al. 2004;Willmott and Field 2006;Radi and Loret 2008;Bizzarri et al. 2010;Schubnel et al. 2011;Barras et al. 2018;Yue et al. 2019;Mai et al. 2020). Despite these results, no theoretical predictions describing Mach cone solutions had been reported in the literature before. In order to derive the solutions reported here, the equations of linear isotropic elastodynamics has been formulated in terms of first order hyperbolic equations for velocity and displacement gradients with their associated jump conditions rather than using the more common second order equation for the displacement field. The solutions obtained describe the steady state propagation of a crack in terms of pressure and shear waves emanating from the crack tip which separate states of constant velocity, strain and stress. In line with the linear elasticity model used, energy is exactly preserved and simply transformed from elastic energy before fracture into a combination of a smaller amount of elastic energy and kinetic energy. This neglects the energy required to create the new surfaces which dominates fracture mechanics in the Griffith theory model. There are a number of physical limitations of the proposed mathematical solutions. First and foremost, the solutions developed are based on linear elasticity and therefore do not include any form of dissipative effects. This includes the absence of surface energy as the crack grows as explained above. In addition, the solutions consider only the steady state situation and provide stress fields that are uniform both ahead and behind the crack. It is obvious that just before the crack begins to grow a complex stress field with high concentrations around the tip will exist, as predicted by classical fracture models based on intensity factors. The transition between these two states is complex and beyond the scope of the analysis provided here. Finally, the state of stresses in the regions beyond the shocks can in some cases imply a state of stress more demanding than the original which would contradict the assumption of the crack propagating along the horizontal axis unless a prior weakness is assumed. In essence the paper has provided a set of very simple particular solutions that satisfy the linear elasticity equations of motion plus far field and crack surface boundary conditions under mode I and mode II fracture. It is likely that similar but more complex solutions with interacting or curved shock fronts will exist, or that through use of superposition techniques with homogenous solutions it will be possible to overcome some of the simplifying assumptions contained in the present work.
Regardless of the above physical shortcomings, the paper sheds light from an analytical point of view on the transfer of elastic to kinetic energy that takes place during brittle fracture. It provides a useful counterpoint to quasi-static fracture models where the strain energy is fully transformed into surface energy through the introduction of a length scale effect. It is likely that brittle fracture for most materials will include a component of kinetic energy transfer as well as a component of surface energy and dissipation. The balance between these different effects will be dictated by the problem configuration and physical considerations and obtained through experimental observations and computational models.
Given the growing number of experimental observations showing intersonic and supersonic crack growth with the presence of Mach type of shock waves very similar to those predicted by the theory above (refer to Abraham and Gao 2000;Gao et al. 2001;Rosakis 2002;Abraham et al. 2002;Buehler et al. 2003;Hao et al. 2004;Petersan et al. 2004;Schubnel et al. 2011;Yue et al. 2019;Abraham 2001), it is clear that high speed crack propagation in these regimes is governed by the transfer of strain energy into kinetic rather than into surface energy or dissipation.
The paper also provides a starting set of boundary and initial conditions for computational simulations where more complex material models and dissipative effects can be incorporated leading to more physically meaningful solutions. This will be explored in future work in the context of Smooth Particle Hydrodynamics (Lee et al. , 2016(Lee et al. , 2019Ghavamian et al. 2021) discretization methodologies where the introduction of dissipative and/or length based effects to account for surface energy, hyperelasticity and inelastic effects will be possible.