Spacetime simulation of dynamic fracture with crack closure and frictional sliding

We combine the asynchronous spacetime discontinuous Galerkin (aSDG) method, an interfacial-damage fracture model, and a dynamic contact model to simulate dynamic fracture and crack closure in brittle materials. The contact model enforces specialized Riemann solutions for bonded, separation, slip and stick conditions while preserving elastodynamic characteristic structure across fracture interfaces. Powerful adaptive spacetime meshing tracks dynamic evolution of fracture-surface networks and captures moving solution features. We present numerical examples to demonstrate the model’s ability to reveal fine details of fracture response in problems that range from dynamic crack initiation, growth, closure, and arrest along a pre-defined planar path to fragmentation of rock by an explosively loaded wellbore with stochastic nucleation, free propagation, and coalescence of fracture surfaces.


Introduction
Despite decades of development, numerical simulation of dynamic fracture remains a challenging and open problem. The problem is inherently multiscale in space and time, requires stochastic modeling, and can involve complex networks of fractures that arise spontaneously and evolve rapidly over time. In this work we focus on dynamic fracture problems that involve dynamic contact, including sliding with friction, in which the contact is driven by crack closure or fracture in the presence of compressive confining stresses. A more complete discussion of available methods for dynamic fracture simulation can be found in [1]. Here we restrict ourselves to brief descriptions of some of the main families of methods used to address this challenging application.
Modeling the growth and dynamically evolving connectivity of fracture networks is one of the key problems in dynamic fracture simulation. In general, there are two approaches to this challenge. Explicit geometry models represent the fracture network directly, typically as the union of segments which may coincide with inter-element boundaries or use other segmented representations. Examples of implementations with explicit geometry models include cohesive zone models [2][3][4][5][6], extended and generalized finite element methods (XFEMs and GFEMs) [7][8][9][10][11][12][13][14][15], and interfacial damage models [1,[16][17][18][19][20]. These methods are well-suited to sharp-interface crack models, but they either admit undesirable constraints on crack paths, restrict their use to linear elastic fracture mechanics (LEFM) theory, or place extreme demands on adaptive meshing to support accurate crack tracking.
Implicit geometry models, such as variational [21], phase-field [22][23][24], and thick levelset [25] methods, have enjoyed recent popularity in computational fracture mechanics. These methods easily navigate changes in crack connectivity and, because the mesh is not required to track the fracture network, they do not require special adaptive meshing techniques. However, they do require expensive mesh refinement to attain suitably narrow crack widths and tend to generate blunt crack-tip profiles that can alter crack-tip fields and introduce conservation errors during crack extension.
In this work, we adopt the asynchronous spacetime discontinuous Galerkin (aSDG) method for elastodynamics [26][27][28] with extensions for dynamic fracture [1,29]. The aSDG method replaces the usual spatial discretization advanced by a temporal integration scheme with a discontinuous Galerkin (DG) finite element discretization of spacetime. An unstructured, asynchronous mesh covers the spacetime analysis domain with clusters of simplicial cells called patches, such that each patch has causal (space-like) boundaries. The causal property in combination with the DG discretization, wherein solution degrees of freedom are private to individual elements, establishes a partial ordering of patches whereby each new patch can be solved locally and without approximation using an implicit Galerkin projection.
We use the Tent Pitcher algorithm [30,31] with adaptive extensions [28,32] to generate causal spacetime meshes. The basic Tent Pitcher method advances the front, a space-like simplicial mesh, by incrementing the time coordinate of one front vertex at a time. It forms a new patch for each vertex advance by generating a small cluster of simplicial cells to cover the spacetime region between the old and new fronts. Adaptive extensions of Tent Pitcher [32] generate special patches that perform common remeshing operations (eg coarsening, edge flips, and mesh motion) within their interiors. Patch generation, patch solution, and adaptive meshing operations share a common granularity and are interleaved in the aSDG algorithm. Detailed descriptions of the formulation and implementation of this method can be found in the cited papers.
The aSDG method's unique approach to adaptive spacetime meshing is sufficiently powerful that it can track dynamic fracture networks with spacetime element boundaries without introducing mesh-dependent constraints or artifacts. Thus, we address the main drawback of explicit geometry models, albeit at the cost of some algorithmic complexity, while avoiding the drawbacks of implicit representations.
Other beneficial features of the aSDG solution scheme include We use the dynamic contact model introduced in [33] to model dynamic crack closure and sliding with friction, the main focus of this paper. As explained below, this method uses Riemann solutions to formulate dynamic aSDG jump conditions that preserve characteristic structure across fracture interfaces for open, stick, and slip crack conditions. Thus, in contrast to the essentially quasi-static contact models in most other fracture models, it enforces dynamic contact conditions for crack closure.
The organization of the rest of this paper is as follows. The next two sections review, respectively, an interfacial-damage fracture model with provisions for crack closure and frictional slip and models for stochastic nucleation, extension, and branching of fracture surfaces. We then present numerical examples that demonstrate the proposed model's capabilities in problems that involve crack closure and frictional sliding. We close with concluding comments that summarize our findings and indicate directions for continuing research.

Interfacial-damage fracture model with crack closure and frictional slip
This section reviews the formulation of the dynamic contact model, first introduced in [33], that we use to model crack closure in this work. Consistent with our sharp-interface aSDG fracture model, we use interfacial damage to model dynamic fracture processes and Riemann solutions to formulate suitable jump conditions across fracture interfaces. The Riemann solutions for crack closure include one for contact-stick and one for contact-slip governed by a Coulomb friction model.

Hierarchy of subscale interface states
We adopt a sharp-interface fracture model and assume that the macroscopic response of fracture interfaces can be approximated as a linear combination of the responses of a finite set of subscale interface states, each of which has a well-defined Riemann solution. This model reflects the fact that we generally observe a mixture of distinct interface conditions at microscopic length scales in the vicinity of any macroscopic location. Photomicrographs show evidence in dynamic fracture process zones of discrete debonding around nucleation sites as well as growth and coalescence into through cracks; cf. [34,35]. Numerical regularization is a separate, non-physical motivation for blending interface states. For example, separation-to-contact transitions in crack closure generally produce discontinuous response that can frustrate the convergence of nonlinear solvers [36][37][38]. Our goal in this work is to approximate the macroscopic interface response generated by mixed interface states without computing and homogenizing subscale solutions based on explicit representations of microscopic geometry. At the same time, we seek a model that preserves the characteristic structure of macroscopic waves impinging on fracture interfaces.
The binary tree in Fig. 1 depicts a hierarchy of subscale states at time, t, in the neighborhood of some macroscopic location, x, on a fracture or potential fracture interface. Three mix parameters describe the ratio of interface states at each level of the hierarchy, or in some cases, serve as regularization parameters for non-smooth transitions between interface states. An interfacial damage parameter, D(x, t) ∈ [0, 1], describes the debonding process at subscale level I, such that D = 0 indicates that the neighborhood is completely  The bonded state has no substates, but a contact parameter, η(x, t) ∈ [0, 1], describes at subscale level II the mix of separation and contact substates within the debonded condition such that η = 0 indicates complete separation and η = 1 indicates full contact in the neighborhood of x. Partition-of-unity interpolations, 1 − η and η, describe the mix of separation and contact for intermediate values of η. These interpolations may regularize discontinuous separation-to-contact transitions for interfaces with perfectly smooth faces or represent the physical effects of microscopic asperities in the surfaces of rough interfaces. Values of η are generally determined by the macroscopic opening across an interface, as described in [33].
The separation condition has no substates, while contact has two substates at subscale level III, slip and stick. A stick parameter, γ (x, t) ∈ [0, 1], describes the mix of slip and stick within the contact state. Here, γ = 0 indicates complete slip, γ = 1 indicates full stick, and interpolation functions, 1 − γ and γ , may be used to describe the mix of slip and stick states for intermediate values of γ , where γ is generally determined by an assumed friction model. In the special case of ideal homogeneous surfaces and isotropic friction response considered in this work, transitions between stick and slip response are continuous and smooth, so no regularization is required. In this case, it is sufficient to treat transitions as instantaneous with no regularization, i.e., γ ∈ {0, 1}. See [33] for more detail on the computation of γ .
The set of basic interface states, S := {B, SE, SL, ST}, contains the four leaf nodes highlighted in Fig. 1,

Riemann solutions for the basic interface states
This section summarizes Riemann solutions from [33] for each of the interface states in S in which Riemann values are decorated by a superposed˘. Let d ∈ {1, 2, 3} be the spatial dimension of the problem. Figure 2 depicts a spacetime fracture interface with a local coordinate frame at an arbitrary spacetime location P for d = 2. The local coordinates Fig. 2 Local coordinate frame at location P on a spacetime fracture interface for a problem in two spatial dimensions (d = 2). Traces of the stress and velocity fields from opposing sides, (S ± , v ± ), determine the Riemann values, (S,v ± ), on are (ξ 1 , ξ 2 , t) in which ξ i and t are, respectively, spatial and temporal coordinates. We use superscripts + and − to denote quantities associated with the regions on opposing sides of , as indicated in the figure. The frame is oriented such that the ξ 1 -direction aligns with e 1 , the unit outward spatial normal vector for the region on the − side of . Traces of the velocity and stress fields from the opposing sides, v ± and S ± , define the initial data for the Riemann problem.
We solve the Riemann problem by preserving characteristic values of the elastodynamic solution across while enforcing balance of linear momentum (BLM) and suitable kinematic constraints on the velocity. If we assume that there are no external forces acting on , then BLM requires continuity of the normal components of the Riemann stress tensor across . That is,S 1i+ =S 1i− for i ∈ {1, . . . , d}. On the other hand, the in-plane components of the Riemann stress tensor (e.g.,S 22± ,S 33± , andS 23± for d = 3) need not be continuous across the interface [6]. However, the in-plane components do not contribute to linear momentum or energy balance across , so we do not need to compute them. Accordingly, only the normal components of the Riemann stress solutions, defined byS 1i :=S 1i− =S 1i+ , are presented below. The kinematic conditions vary with the basic interface states as follows. Velocity is continuous across for the bonded and stick cases; i.e., v :=v + −v − = 0. For the slip case, the impenetrability condition requires continuity of the normal velocity components, v 1 :=v + 1 −v − 1 = 0, while possible slip admits discontinuous (unconstrained) tangential velocity components. In separation mode, all velocity components may be discontinuous, so there are no constraints betweenv + and v − .
For isotropic materials in linear elastodynamics, the characteristic trajectories in all directions on both sides of the interface are determined by the local dilatational and shear wave speeds, c ± d and c ± s , in which, where ρ is mass density and λ, μ are Lamé parameters. The wave speeds determine the impedance components, in which the index i corresponds to the spatial directions in the local frame for up to three spatial dimensions, cf. Fig. 2 for the case where d = 2. The Riemann solutions for the bonded and stick interface states are identical, and their velocity and normal stress components must both be continuous across . Thus, we writȇ in which no summation is implied for repeated indices i. In separation mode,v + andv − are independent, while the normal components of the Riemann stresses are determined byS SE (e 1 ) =s wheres :=s − = −s + in which s ± are tractions acting on the regions on the + and − sides of . For example, if the separated interface is unloaded, we haves = 0 as in [33]. Alternatively, in hydraulic fracturing applications,s is the traction induced by hydraulic or explosive pressure, as in the explosively loaded wellbore problem in the "Numerical examples" section. The separation Riemann values for anys arȇ A Mohr-Coulomb friction model governs transitions between stick and slip states in this work, although alternative friction models could be substituted. Thus, a debonded interface in contact enters the slip state when the normal components of the Riemann stress for assumed stick conditions satisfy the Coulomb condition, in whichS 11 ST e 1 andτ ST are the normal and tangential vector components of the traction vector induced by the Riemann stress tensor acting on the -side of for assumed stick conditions, k is the friction coefficient, . + is the positive Macaulay bracket, and When the slip state prevails, the magnitude ofτ SL is given by k −S 11 ST + , and its direction is determined by the interfacial slip velocity. However, the slip-velocity direction suffers a discontinuity when the slip velocity vanishes, and this is known to cause problems at stickslip transitions in numerical simulations [36][37][38][39]. We showed in [33] that, for isotropic friction relations, the interfacial slip velocity andτ ST agree to within a positive scaling, and this relation holds even when the slip velocity vanishes. Therefore, we replace the slip-velocity direction with e τ :=τ ST /|τ ST | to determine the direction ofτ SL .
Enforcing the stick Riemann solution for the normal part and the Mohr-Coulomb friction law for the tangential part, we obtain the slip-mode Riemann solutions, in whichv ST 1 is the normal velocity of the Riemann solution for assumed stick conditions and there is no sum on i.

Macroscopic response model
We approximate the macroscopic response as a weighted average of the basic interface state responses. To this end, we introduce a vector of weights, a := (a α (D, η, γ )) ; α ∈ S, based on partition-of-unity interpolations at each level of the state hierarchy. We combine the partition-of-unity interpolations across all three levels to obtain, This use of partition-of-unity interpolations ensures that (i) our macroscopic approximation exactly matches the subscale conditions when all three mix parameters assume either extreme value and (ii) the weights sum to unity for arbitrary combinations of the mix parameters: Balance of linear momentum, (9), and preservation of characteristic structure across interfaces leads to simple weighted average expressions for the macroscopic Riemann stresses, S, and velocities, v ± ; cf. [33].
where only the normal components of the stress tensors in equation (10a) need be considered. The Riemann solutions for the bonded/stick (B/ST), separation (SE), and slip (SL) subscale interface states appear, respectively, in (3), (4), and (7). The form of (10b) accommodates velocity (and displacement) jumps across fracture interfaces in separation and slip modes; cf. (4b) and (7b).

Evaluation of the mix parameters
In this subsection, we summarize methods for evaluating the mix parameters, (D, η, and γ ), that control the participation of the subscale interface states in determining the macroscopic response.

Interfacial damage parameter, D
We initialize D = 0 on new, undamaged interfaces and, following the model in [40], we adopt the damage evolution equation, in whichτ is a relaxation time, andD is a target damage value. In general, the function H has unit value at zero and decreases monotonically to 0 at infinity. As in [40], we use H(x) = exp(−ax), a form that enforces a maximum damage rate of 1/τ . We focus on mechanical damage processes and assume that Riemann stresses for the undamaged, bonded state drive the interfacial damage evolution. As in [3,19,41], we introduce a scalar effective traction,š, defined as a function of the normal and tangential vector components of the traction induced by the bonded Riemann stress (cf. (3a) and (6)),š := S 11 in which the Macaulay bracket ensures that only tensile normal traction components drive damage evolution, and the shear factor, β, adjusts the influence of the tangential component. The target damage value is then written as a function ofš, in which 0 < s <s, and s ands denote, respectively, thresholds for the onset of additional damage evolution and for attainment of the maximum damage rate, 1/τ . Numerous experimental studies demonstrate that fracture strength and fracture energy are rate-dependent properties [42][43][44]. Typically, increasing loading rates cause both properties to grow. Various cohesive fracture models attempt to capture these rate dependencies; cf. [45,46] and the references therein. The time-delay format of the damage evolution equation (11) introduces rate dependency for strength and energy, similarly to the model in [47]. Fracture strengths and energies in our model vary with position on fracture surfaces according to the loading history at each location and stochastic variations in the damage evolution parameters, as described in the following section. Nonetheless, the fracture energy for various loading rates and for given fracture parameters has been shown to be proportional to a fracture energy scale,G =τs 2 /Z d [48].

Contact parameter, η
The contact parameter, η, can serve either as a regularization of discontinuous response across separation-contact transitions in the case of ideally smooth interfaces or as a physically motivated macroscopic model for gradual transitions in the case of rough interfaces. In this work, we focus on the former case in which transitions between contact and separation occur instantaneously. A binary range, η ∈ {0, 1}, would suffice to describe the physical transition in this idealized model. However, the physical response predicted by the Riemann solutions for separation-to-contact (either slip or stick) transitions is discontinuous, and this can cause convergence problems in numerical simulations. On the other hand, contact-to-separation transitions do not suffer this problem. We therefore use η ∈ [0, 1] as a regularization parameter, but only in the case of separation-to-contact transitions. We use the regularization described in [33] in which η is computed as a function of the normal displacement jump across the interface and the normal component of the traction induced by the Riemann stress for assumed stick conditions; cf. (3a) for i = 1. Please refer to [33] for a more complete discussion of the contact-separation model and algorithm.

Stick parameter, γ
Our use of (7a) to determine the direction ofτ SL for isotropic friction models circumvents the discontinuity and ensuing numerical problems suffered by methods that use the interfacial slip velocity to determine the direction. In fact, the response predicted by (7) is continuous across stick-slip transitions in either direction, so there is no need to introduce a regularization. We therefore restrict the stick parameter to a binary range, γ ∈ {0, 1} and allow instantaneous changes governed by the Coulomb condition, (5).

Stochastic nucleation, crack extension, and crack branching
Cracks typically nucleate in macroscopically homogeneous material at microscopic flaws on interfaces, such as grain boundaries and material interfaces, or in the bulk at preexisting voids, inclusions, and microcracks. Failure to account for random distributions of the severity and orientations of these flaws generally leads to significant over-estimates of fracture resistance. This section summarizes a stochastic model for fracture surface nucleation and adaptive aSDG methods for modeling dynamic crack extension and coalescence that were first presented in [1]. The reader should consult that publication, in particular the implicit realization of microscopic flaws, for a more complete development.
The term fracture surface has a special meaning in this work that is consistent with our use of the interfacial damage parameter, D, to model debonding processes, as described in the preceding section. Fracture surfaces are initialized with D = 0 and represent material interfaces on which D might (or might not) evolve toward full debonding. Thus, they have a fully bonded condition at initialization and their insertion into bulk material does not immediately alter that material's response. A fracture surface only attains the behavior of a physical crack in regions where the evolution rule, (11), drives its damage to D = 1. Regions where the damage parameter takes intermediate values, 0 < D < 1, are partially debonded and correspond to active fracture process zones.

Nucleation of fracture surfaces
Fracture surface tips (FSTs) are vertices in the spacetime mesh from which new fracture surface segments may be extended. An FST is not equivalent to a crack tip because we set D = 0 when FSTs are created and fracture surfaces typically extend to a new FST before D reaches unity at an existing FST. Nucleation of a new fracture surface involves designation of a new isolated FST and subsequent addition of a fracture-surface segment that emanates from that point. Extension of existing fracture surfaces involves adding a new segment emanating from one of its FSTs, and replacing the old FST with a new one at the end of the new segment.
A directional effective stress,š(θ ) cf. [3], controls both nucleation and fracture surface extension. For two spatial dimensions, we havě s(θ) := s 1 (θ) 2 + + β 2 s 2 (θ ) 2 (14) in which θ is an angle that describes the continuously variable orientation of the interface associated with the traction vector, s(θ) direction. We use stress solutions on the interiors of elements impinging on the candidate FST vertex to evaluate s(θ ) which, in view of our DG model, is continuous almost everywhere for θ ∈ [0, 2π]. We use a probabilistic criterion to create FSTs that nucleate new fracture surfaces. We use a Weibull-type probability distribution function (PDF) to model microscale flaws with random strengths, and an inverse CDF method to sample the PDF to determine angleindependent flaw strengths,s, at each spacetime-mesh vertex, as explained in [1]. An isolated FST is created at any vertex where any of the sampled strengths satisfyš(θ ) >s for any angle θ ∈ [0, 2π]. Nucleation of a new fracture surface is completed by extension of a new fracture-surface segment emanating from the isolated FST.
The procedure for generating a new fracture-surface segment that emanates from an FST is essentially the same, whether the FST is isolated, as in nucleation, or the FST is an endpoint of an existing fracture surface. We generate new segments from any FST where max θ ∈[0,2π ]š (θ) ≥s, and the direction of extension is determined by argmax θš (θ ).

Tracking extensions of fracture surfaces
Our choice of a sharp-interface representation for dynamically propagating cracks presents several daunting challenges in numerical implementation. In general, we must align spacetime element boundaries with crack trajectories that nucleate and extend dynamically and whose paths are unknown prior to solution. In addition, numerical solutions for crack trajectories should depend only on the elastodynamic solution and on the given continuum model for crack extension. Thus, in the absence of stochastic nucleation effects, fracture surface trajectories should converge in the limit of mesh refinement (or vanishing error tolerances for adaptive meshing).
We focus on three requirements for meeting these objectives in the context of our proposed algorithm. First, the aSDG elastodynamic solutions must converge in the vicinities of FSTs and crack-tip process zones. We depend on the aSDG method's powerful hadaptive meshing capabilities to meet this requirement. Solution convergence is ensured and enhanced through the use of newest-vertex refinement, spacetime edge flips, and mesh-smoothing via tilted tent poles to continually preserve and improve mesh quality [32].
Second, the directions of incremental fracture-surface extensions must exactly match argmax θš (θ), irrespective of the current mesh layout. We accomplish this by invoking one of two remeshing options to provide an element edge in the specified direction. Edge division first creates a new FST vertex at the intersection of the extension direction and the opposite edge of a triangular facet impinging on the old FST. It then creates a new edge between the old and new FSTs to subdivide the original triangle into two. Alternatively, when the extension direction is close to an existing inter-element edge, we create the new FST by repositioning a nearest-neighbor of the old FST onto a line coincident with the extension direction. Details of these operations are presented in [1,32].
Finally, the lengths of fracture surface extensions must vanish in the limit of mesh refinement. This ensures that our incremental crack-path approximation converges to its corresponding continuum solution. Our algorithm satisfies this requirement automatically because the lengths of fracture-surface segments, generated either by edge division or repositioning, scale with element diameter. We demonstrated crack-path convergence with reducing adaptive tolerances in an example with stochastic nucleation disabled in [1]. Beyond its capabilities for accurately tracking fracture surface trajectories, our implementation properly handles fracture-surface extensions that intersect the domain boundary or another fracture surface; i.e., crack coalescence.

Modeling crack branching
As with many other aspects of fracture mechanics, dynamic crack branching is a multiscale phenomenon. The earliest experimental and theoretical studies of dynamic crack branching adopted a macroscopic perspective in which branching is viewed in the framework of linear elastic fracture mechanics (LEFM) as a form of instability in which a single crack suddenly bifurcates into two branches (or splits into multiple branches) [49][50][51][52][53]. In general, the goal of these studies was to identify critical values of the dynamic stress intensity factor, the crack-tip velocity, or other LEFM parameters that predict the onset of branching. Most numerical methods for modeling crack branching also rely on macroscopic representations. Most often crack branching arises naturally from the geometric flexibility of the representation without reference to critical LEFM parameters. Examples include methods based on intrinsic cohesive models [2], phase field models [22], peridynamics [54], and gradient damage models [24]. We also note molecular dynamics models [55] that predict qualitatively similar branching behavior.
Despite the emphasis on macroscopic resolutions, underlying micro-mechanical mechanisms were hypothesized and supported by experimental studies in some of the earliest research on dynamic crack branching [49,56,57]. The microscopic scenario for crack branching begins with nucleation of micro-cracks from voids, inclusions, or other flaws under the influence of the crack-tip field of a dynamically propagating crack. Micro-cracks that nucleate off the fracture plane propagate and coalesce with the main crack to form micro-branches. These may arrest without forming macroscopic crack branches. How-ever, at higher levels of the dynamic stress intensity factor, some of the micro-branches propagate away from the main crack to form new macroscopic branches. This micromechanical scenario, rather than macroscopic bifurcation, is the basis of our numerical implementation of crack branching.
Rather than model individual microscopic flaws explicitly, we use the probabilistic nucleation model to represent their collective influence on micro-crack nucleation. The method for extending fracture surfaces then handles micro-crack propagation and coalescence to form micro-branches and, in some circumstances, new macroscopic crack branches. Thus, our method's basic capabilities to model crack nucleation, extension, and coalescence suffice to model crack branching; cf. numerical examples in [1] and the "Numerical examples" section below. The success of this approach depends on numerical solutions that resolve the smaller length-scales associated with dynamic propagation of micro-cracks. The powerful adaptive capabilities of the aSDG method are essential to satisfying this requirement.

Numerical examples
This section presents numerical results that demonstrate the proposed method's capabilities. Three examples demonstrate the ability of the adaptive aSDG method, in combination with the interfacial-damage fracture model, to capture fine details of bulk elastodynamic response as well as crack initiation, arrest, opening, and closure for both stick and slip contact modes. The third example demonstrates stochastic nucleation and crack propagation along solution-dependent paths. All three examples employ linear plane-strain models of elastic response.
In all three examples, we use h-adaptive spacetime meshes consisting of tetrahedra, each with a complete cubic polynomial basis in 2d × time. We drive our adaptive spacetime meshing algorithm with three independent criteria to ensure the reliability of our solutions. Although the aSDG formulation conserves linear and angular momentum to within machine precision, it does not directly enforce energy balance. Therefore, similar to the methods in [29], we use distinct error indicators to limit numerical energy dissipation in the bulk and across fracture surfaces. We use a third residual-based error indicator to ensure accurate integration of the damage evolution equation (11) on fracture surfaces. We only accept patch solutions that satisfy all three error criteria.
The sequences of solution visualizations in Figs. 4, 6, and 8 were generated by the perpixel-accurate rendering procedure described in [58]. The log of the strain-energy density in these sequences maps to color, where blue indicates low energy density and violet indicates peak values. The height field depicts the modulus of the material velocity.

Planar fracture with crack closure under far-field cyclic loading
This example involves dynamic crack growth along a straight, predefined crack path driven by far-field cyclic loading. Its purpose is to test the proposed method's ability to model fracture initiation, propagation, and arrest as well as intermittent crack closure with separation-to-contact transitions. Figure 3 depicts the computational domain and boundary conditions for this two-dimensional problem. The rectangular domain has length L = 10 mm, width 2W = 6 mm, and a predefined fracture surface along the full length of its horizontal centerline. We assign an initial damage value, D 0 = 1, to represent the debonded condition of a pre-crack that extends a distance a 0 = 4.25 mm from the left edge of the fracture surface. We assign D 0 = 0 to the remainder of the fracture surface to indicate its initial fully-bonded condition.
We assign homogeneous displacement and traction boundary conditions along the bottom and left edges, respectively. We model transmitting boundary conditions along the right and top edges by prescribing piecewise uniform inflow characteristic values, ω(t), while leaving the outflow characteristic value unconstrained. This allows waves impinging on the top and right edges to exit the domain without reflections. We write ω(t) =ω(t) along the top edge to model cyclic far-field characteristic loading with a 2 μs period and maximum tensile and compressive characteristic values of 40 and 160 MPa, as shown in Fig. 3b.
The initial front mesh in this example consists of only eight triangular cells over the entire computational domain. Nonetheless, the aSDG adaptive meshing scheme automatically generates a strongly graded spacetime mesh as the solution proceeds to accurately resolve all moving wave fronts, crack-tip process zones, and contact-mode transitions. Figure 4 presents a sequence of solution visualizations. Figure 4a shows the solution at a time slightly later than when the first tensile wavefront reaches the fracture surface at t = 1.43 μs. A strong diagonal wake emanates from the left free boundary as the tensile wave passes. The tensile field around the pre-crack tip has grown sufficiently strong by the time of Fig. 4b to initiate crack propagation to the right of the pre-crack tip. We observe crack opening along the length of the pre-crack, except where it is intersected by the diagonal wake of the tensile wavefront which now approaches the bottom edge of the domain. Meanwhile, the first compressive wavefront appears after entering the domain through the top edge. Spikes in the height field around propagating crack tips in Fig. 4b, f, and especially in Fig. 4g indicate quasi-singular velocity response. This effect intensifies as fracture process zones shrink with increasing crack-tip velocity so that our nonlinear interfacial damage model better satisfies the small-scale yielding assumption of linear elastic fracture mechanics theory. See [29] for further discussion of quasi-singular response. Figure 4c shows a time just after the first compressive wave reaches the fracture surface, and Fig. 4d shows a slightly later time when the compressive wave has crossed the bonded fracture surface ahead of the crack and interferes additively with the inverted reflection of the initial tensile wave. In contrast, the compressive wave is reflected by the open crack on the left half of the fracture surface. The arrival of the compressive wave immediately arrests the propagating crack tip, and we observe an expanding circular quiescent zone that is concentric with the arrested crack tip, a signature feature of dynamic crack arrest. We also observe waves that emanate from expanding zones of crack closure on the fracture surface. These contact zones expand, to the left from the arrested crack tip and to the right from the diagonal wake, until they coalesce.
More of the fracture surface is in contact in Fig. 4e, and crack-closure waves have scattered further from their source contact zones. A subsequent tensile wave causes further crack propagation in Fig. 4f. The background wave pattern becomes increasingly complex as waves continue to enter through the top edge and reflect off the fixed boundary. These produce additional episodes of crack propagation, crack arrest, as well as crack closure and opening, as seen in Fig. 4g, h. Overall, this example demonstrates our model's ability to capture fine details of fracture initiation, propagation, and arrest as well as seamless transitions between separation and stick modes during crack closure.  We restrict brittle, fracture-like debonding to nucleate and propagate along the matrixinclusion interface. Initially, the interface is perfectly bonded (D 0 = 0) everywhere except for an interval from 17 • to 22 • , measured clockwise from the top of the inclusion, where we model an initial flaw with fully debonded conditions (D 0 = 1). Symmetry boundary conditions along the left and right edges of the domain simulate an infinite strip of cells with identical inclusions. In contrast to a similar example with non-cyclic loading in [29], we impose spatially uniform cyclic boundary conditions on the top and bottom edges where the prescribed normal velocity cycles between a maximum tensile value,v = 7 m/s, and a peak compressive value,v = −14 m/s, with a period of 4 μs, as shown in Fig. 5. Figure 6 shows a sequence of solution visualizations for this problem. Figure 6a depicts a time slightly after the first tensile wavefront arrives at the initial flaw where crack-tip fields begin to develop. The mismatch in elastic moduli between the matrix and the inclusion causes asymmetric response across the interface. Figure 6b shows a later stage where the entire interface is debonded and the arrival of a compressive wave begins to close the interface gap, first at the bottom of the inclusion and soon after at the top. These abrupt separation-to-contact transitions send higher-speed waves into the stiffer inclusion and slower waves into the matrix. The top and bottom contact regions expand along the interface, as seen in Fig. 6c, until the contact zone covers the entire interface in Fig. 6d. The entire interface remains in contact mode in Fig. 6e-i due to the long-duration far-field compressive loads. Meanwhile, a complex wavefront pattern develops due to reflections and partial transmission at the boundary of the inclusion and the top and bottom edges as well as waves entering the simulation domain from adjacent cells through the left and right edges.
The complete spacetime mesh up to the terminal time, t = 4.8 μs, contains 255.8 million tetrahedra arranged in 50 million patches. This h-adaptive discretization yields a total (spurious) numerical energy dissipation of 0.15% of the net energy inflow from the prescribed velocity loading. The total fracture energy error is 0.099% of the net energy inflow.

Explosively loaded wellbore in rock
This example demonstrates the full range of our model's capabilities, including stochastic nucleation, propagation, and coalescence of solution-dependent fracture surfaces as well as mixed-mode transitions between bonded, slip, and stick interface states. The adaptive spacetime meshing must now track dynamically evolving fracture surfaces in addition to capturing moving wavefronts and crack-tip fields-all while avoiding mesh-dependent influence on the solution. Figure 7 shows a wellbore in rock subjected to far-field confining stresses. The wellbore has a diameter of 30 cm and is subjected to hydrostatic in-situ confining stresses,   We select contact regularization parameters,δ = 0 and δ = 805 nm, so that δ = 0.1. An explosive compressive pulse acts on the wellbore walls, ramping from ambient pressure to 27.5 MPa in 750 ns and then back to ambient pressure. Figure 8 shows a sequence of solution visualizations. The in-situ confining pressure precludes mode-I failure in the rock, so we observe mostly mode-II fracture in which Fig. 9 Detail of adaptively refined front mesh (black) and damage values along fracture surfaces (color) near the wellbore at t ≈ 37.5 μs. Damage in the range, D ∈ [0, 1], maps to a blue-to-red color range with red indicating full damage tangential stress is the primary contributor to the effective stress; cf. (12). A circular pressure pulse expands from the wellbore at the dilatational wave speed in Fig. 8a, b while early stochastic nucleation of fracture surfaces breaks the solution symmetry. The radial expansion of the failure zone lags the circular pressure wave by a significant margin, as seen in Fig. 8c-f. This is expected because the Rayleigh wave speed limits the speed of crack propagation and because the propagation directions of the mode-II cracks are inclined relative to the radial direction. Crack branching and coalescence, inter-crack shielding, and the finite time required for interfacial damage evolution, cf. (11), further impede the radial expansion of the failure zone. As a consequence, many fracture surfaces that nucleate and begin to propagate as the pressure wave passes by never attain full damage. Figure 9 shows a detail of the asynchronous front mesh (in black) and the distribution of damage on fracture surfaces near the wellbore at t ≈ 37.5 μs. The circular pressure wave has moved beyond this figure's field of view at this stage, and the adaptive mesh refinement in this figure reflects the remaining network of fracture surfaces. Elsewhere, outside the fracture zone, the mesh has coarsened behind the pressure wave to improve solution efficiency. Fracture surfaces are displayed in colors that reflect the local damage level; blue indicates near-zero damage and red indicates full damage. We observe a complex network of curved fracture surfaces whose paths are free of mesh-dependent features. We observe higher damage levels close in to the wellbore, with lower levels further out. This confirms our expectation that damage evolution and fracture surface propagation slow to a halt as the circular pressure wave outruns the expansion of the failure zone.

Conclusions
We combined an adaptive aSDG solution method for elastodynamics and the interfacialdamage fracture model of [1] with the contact model of [33] to obtain an effective model for dynamic crack closure. We enforced specialized Riemann solutions to ensure preservation of hyperbolic characteristic structure across fracture surfaces for bonded, separation, slip and stick interface states. We used powerful spacetime adaptive meshing techniques to capture moving solution features at multiple length and time scales and track the free evolution of fracture-surface networks without imposing mesh-dependent constraints.
Numerical examples demonstrated the model's ability to reveal fine details of brittle fracture response in example problems of increasing complexity, ranging from dynamic crack initiation, growth, closure, and arrest along a pre-defined planar fracture surface to rock fragmentation by an explosively loaded wellbore including stochastic nucleation, free propagation, and coalescence of fracture surfaces.
All computations in this work were performed as serial calculations and did not exploit the intrinsic parallel structure of the aSDG solution scheme. An effort to implement fracture simulation within a parallel aSDG framework is in progress. Extension of our aSDG simulation method to 3d×time is essential to widen the range of scientific and engineering applications for this technology, and this is a current subject of active research and development.
The work reported here was intended to demonstrate numerical capabilities using simple constitutive relations. We made no attempt to calibrate our model to predict the measured response of real materials as in, for example, [14]. Research on improved material models and validation with experimental data are important goals for continuing work.