ANALYSIS OF A METHOD TO COMPUTE MIXED-MODE STRESS INTENSITY FACTORS FOR NON-PLANAR CRACKS IN THREE-DIMENSIONS

. In this work, we present and prove results underlying a method which uses functionals derived from the interaction integral to approximate the stress intensity factors along a three-dimensional crack front. We first prove that the functionals possess a pair of important properties. The functionals are well-defined and continuous for square-integrable tensor fields, such as the gradient of a finite element solution. Furthermore, the stress intensity factors are representatives of such functionals in a space of functions over the crack front. Our second result is an error estimate for the numerical stress intensity factors computed via our method. The latter property of the functionals provides a recipe for numerical stress intensity factors; we apply the functionals to the gradient of a finite element approximation for a specific set of crack front variations, and we calculate the stress intensity factors by inverting the mass matrix for those variations.


Introduction
The stress intensity factors, which characterize the stress singularities near the front of a three-dimensional crack, are important parameters for predicting the failure of engineering structures.In [17], a method to compute the stress intensity factors along the front of a three-dimensional crack was introduced; in this paper, we prove convergence and error estimates for such method.
In the literature, there are a variety of approaches to compute the stress intensity factors, which generally fall into one of two categories.The first category includes extrapolation-based methods.These methods sample the stress or displacement field around the crack and fit known asymptotic behavior (cf.[35]).The stress intensity factors are estimated from the fitting parameters [7,33].In the second category are methods where the stress intensity factors (or combinations thereof) are the outputs of certain functionals applied to the displacement field or its gradient.Among these are methods based on the -integral [8,28] such as the   -integrals [19,27], the Contour Integral and Cutoff Function Methods [4,31,32] and the Quasi-Dual Function Method [12,35], and the interaction integral [34].Most of these approaches originated in the study of two-dimensional crack problems, but have since been extended to three-dimensional cracks (e.g., [16] shows an approach based on the interaction integral for three-dimensional cracks).
The method in [17] is based on the solution of a variational problem involving a set of three functionals {  }  , which are derived from the interaction integral [34] and particularized for the elasticity problem of interest in the continuous setting.These functionals act on a virtual (normal) extension of the crack front  and a square-integrable tensor field.They have two important properties.First, when the tensor field is the gradient of the exact solution of the elasticity problem ∇, they output precisely integrals of the stress intensity factors {  }  along the crack front ℱ weighted by the virtual extension : [, ∇] = ∫︁ ℱ   d,  = I, II, III. (1.1) The functionals (and the interaction integral) provide an integral representation of the stress intensity factors in terms of the solution , which make them particularly suitable for numerical evaluation.Second, the functionals are continuous and affine with respect to their arguments.Two discretizations are needed for the method.First, we consider a fixed virtual extension  of the crack front, and introduce a numerical approximation of the exact displacement gradient with a discretization length scale ℎ  , termed ∇ ℎ  .In this paper and in the implementation [17], the approximate gradient is computed using the Finite Element Method (FEM, e.g., that of [18]) on meshes of the problem domain with mesh size ℎ  .Other methods, such as the Extended Finite Element Method (XFEM [26]) or Mapped Finite Element Method (MFEM [10]), may be used instead.For a fixed virtual extension of the crack front , convergence of the numerical gradient in the  2 -norm guarantees the convergence of the values of the functionals solely due to continuity.Second, we introduce a discretization of the virtual extensions of the crack front.In this case, we restrict ourselves to a finite dimensional subset with length scale ℎ  .For example, we may use piece-wise linear Lagrange finite elements over the crack front with mesh size ℎ  .Alternatively, we can build a spectral basis up to maximum order   (where   is treated like ℎ −1  ).By restricting the numerical stress intensity factors { ℎ  }  to the same space of the virtual displacements and endowing such space with the  2 -scalar product over the crack front, i.e., the right-hand-side in (1.1), each  ℎ  follows as the Riesz representative of the functional   [•, ∇ ℎ  ] in that space.This enables the computation of approximate stress intensity factors by solving a variational problem.Further, for each numerical stress intensity factor  ℎ  we prove an error estimate of the form for constants  1 and  2 independent of ℎ  and ℎ  .The first term emerges from an interpolation error, and the second one from a consistency error.Success of the method hinges on these two errors converging to zero as ℎ  and ℎ  do.What complicates this effort is that the consistency error grows with decreasing ℎ  .For low-order finite elements and scaling the discretizations like ℎ  ∼ ℎ  , this estimate does not guarantee that the method converges.This is the case, for example, when using the restrictions of the shape functions in the three-dimensional volumetric (or bulk) mesh to the crack front as the basis for the space of virtual extensions of the crack.In [17], such scaling between bulk and crack front meshes resulted in reduced rates of convergence of  ℎ  in the  2 -norm.Convergence is guaranteed by shrinking ℎ   more quickly than ℎ  , where  is the order of the approximation ∇ ℎ  (for FEM,  = 1/2).Because the meshes used in the bulk and along the front are different, we term our approach the Multiple Mesh Interaction Integral method (MMII).
In this work, we prove two properties of each functional   , namely continuity and (1.1), and we prove the error estimate (1.2).A critical ingredient for the proof of (1.1) is to show that any smooth part ( 1 ) of the displacement gradient belongs to the kernel of the interaction integral.To the authors' knowledge, such a result for the interaction integral in three dimensions is novel, though a sketch proof was provided for a similar result in two dimensions in [9].We then make use of the two properties of the functionals to derive (1.2), which stems from the variational formulation.This paper is organized as follows.In Section 2, we introduce some preliminary definitions of the geometry in the vicinity of the crack front.We present the linear elasticity problem of interest, and state a key assumption about the decomposition of the solution into a smooth part and a part containing the familiar 1/2 asymptotic behavior of Linear Elastic Fracture Mechanics (LEFM).In Section 3, we recapitulate the functionals and problem to define the approximate stress intensity factors, and we state as theorems the main results of the paper.Section 4 is devoted to proving the theorems, though ancillary and more technical results are relegated to the appendices.We finish with a numerical study of the convergence of the interaction integral in Section 5. Here, we address two key points.First, in computer implementations of the MMII, the functionals {  }  are approximated through numerical quadrature.The integrands of the functionals contain radial singularities like  −1/2 , and, depending on the choice of function space over the crack front, may also have discontinuities.These issues affect the convergence rates of standard quadrature rules under refinement of the bulk mesh, and we assess whether these errors interfere with the analytically-derived error estimate (1.2).Second, classical analysis of continuous linear functionals applied to finite element solutions (e.g., [3,4]) utilizes a duality argument to prove superconvergence (often double that of the finite element error).We recapitulate these arguments in detail, and we explain where analytical shortcomings may arise in their application to our present work.
Throughout this paper, we will make use of constants  (or  1 ,  2 , . . .if we wish to differentiate constants) whose value may change from line to line.We will also refer to the Sobolev space  , (Ω), the space of functions over the domain Ω with weak derivatives up to order  in   (Ω), as well as the Hilbert space   (Ω) :=  ,2 (Ω).We denote the norms for these spaces with ‖ • ‖ ,,Ω and ‖ • ‖ ,Ω := ‖ • ‖ ,2,Ω , respectively.For the space   (Ω), we will also make use of the inner product (•, •) ,Ω .When functions are vector-or tensorvalued, we will specify the range when writing the function space, e.g.,  , (Ω; R 3 ) or  , (Ω; R 3×3 ), though we omit the range when writing the norm or inner product.For a function  of a scalar variable , we denote the derivative by  ′ (),  , , or d /d.For a function  ( 1 ,  2 ), we write  ,1 or  / 1 to denote the partial derivative with respect to  1 .Lastly, we may express points  ∈ Ω ⊂ R 3 in terms of their Cartesian coordinates:

Preliminaries
In this section, we describe the geometry in the vicinity of a three-dimensional crack front, see Fig. 1.We then state the linear elasticity problem in a cracked domain.At this point, we make a key assumption about the elasticity problem -that the displacement field near the crack front may be expressed as a sum of a tip part containing the LEFM  1/2 asymptotes and a smooth part [11,35].We further assume that the tip part may be decomposed into the three stress intensity modes [22,37].

Near-front coordinates
A bounded domain Ω ⊂ R 3 with Lipschitz boundary contains a sharp crack  ⊂ Ω.We assume the crack is an orientable smooth manifold with boundary 1 .We select an orientation with unit normal  , and let  ± denote the crack faces where  points toward  + .We denote the crack front by ℱ = , which we assume is a closed, simple, regular, smooth curve.We assume dist(ℱ, Ω) > 0 to avoid cases where the crack front intersects the surface, such as a through crack.
Let  = len(ℱ) be the length of the crack front, and let   : [0, ] → R 3 be the arc length parameterization of ℱ from an arbitrary starting point.If  : [0, ] → R 3 is the unit tangent vector to ℱ and  :  → R 3 is the unit normal field over , then for any  ∈ [0, ], we introduce the orthonormal basis Figure 1.Left: three-dimensional domain Ω containing an internal crack  with crack front ℱ. Right: local description of the crack surface near the crack front.Points near the crack front may be described using coordinates (, , ) along with the crack front basis {  }.The crack surface is parameterized via the function   , from which we define the crack surface basis {  }.
The inclination angle  is used to specify limits of .This figure has been reproduced from [17].

Elasticity problem
Assuming linear elasticity theory [5], we seek the displacement field  throughout the cracked body Ω  = Ω∖ which satisfies equilibrium.In particular, the body is subjected to body force  in the interior of the domain, tractions  on the Neumann boundary   Ω  , and prescribed displacements  on the Dirichlet boundary   Ω  .We assume   Ω  ∪   Ω  = Ω  ,   Ω  ∩   Ω  = ∅, and  ± ⊆   Ω  (i.e., we do not prescribe displacements on the crack faces).We relate the Cauchy stress tensor  to the displacement gradient  = ∇ via the isotropic constitutive relation: where  and  are the Lamé parameters (though we will also refer to Young's modulus  and Poisson's ratio ).Here and throughout this manuscript, we make use of Einstein summation convention for repeated Roman indices.We now state the primal elasticity problem. (2.10) We make the following assumption regarding the solution to Problem 2.1.
Assumption 2.2 (Near-front decomposition of the elasticity solution).We assume that the solution to Problem 2.1 may be decomposed as for any  ∈    with tubular coordinates (, , ).The function   ∈  2 (Ω  ; R 3 ), while {  }  ⊂  2 (ℱ) are termed the stress intensity factors of  for modes  = I, II, III.The functions {   } , ⊂  ∞ (R) (provided in Appendix B of [17]) are of the form  1 cos , where the constants depend only on the elastic moduli.
In general, it is not known if the decomposition (2.11) always holds under the regularity assumptions of Problem 2.1.One result comes from Costabel et al. [11], in which the authors prove that, for an infinite domain Ω = R 3 containing a smooth crack with closed, smooth front, and with body force  ∈  ∞ (R 3 ; R 3 ), as  → 0 the solution has the decomposition for any integer  ≥ 0. The functions    () belong to  ∞ (ℱ), while the vector-valued functions    (, ) depend only on the linear operator and the type of boundary conditions prescribed on the crack faces.The smooth part  , belongs to  +1 (Ω; R 3 ).
We note the -dependency in the vector functions    (, ) in (2.12), which we have not assumed in (2.11).Instead, we have followed the common assumption in the literature for asymptotic expansions around a crack front that    (, ) =    ()  ().For example, Leblond and Torlai [22] assumed a similar decomposition for the stress field when analyzing the leading-order stresses near an arbitrary three-dimensional crack.Meanwhile, Yosibash et al. [36] computed an asymptotic expansion for the displacement field near cracks with a circular crack front of radius .Yosibash et al. assumed an asymptotic expansion of the form (cf. [36], Eq. (64)) For a penny-shaped crack with arbitrary loading ( [36], Eq. ( 94)), the lowest-order displacement terms in the asymptotic expansion coincide precisely with (2.11).Lastly, in [12], Costabel et al. derive the asymptotic structure for an infinite straight edge.Here, the eigenfunctions for the leading-order terms do not possess -dependency.Without a result that explicitly states the dependence of {   }  on , the definition of the stress intensity factors needs to be revised.Assumption 2.2 allows us to proceed.Of course, our results apply for displacements  of the form (2.11).
The form of  in (2.11) and (2.12) corresponds to the eigenfunction expansion in the vicinity of a threedimensional edge, which we briefly summarize here.A more complete treatment may be found in Yosibash [35].The functions {   } , are determined from the following problem around the crack front.
For such geometry, the tubular coordinates are given by the usual cylindrical coordinates.For each mode  = I, II, III, () =   We say that the functions {   } , are the angular variation of the asymptotic solution.These functions belong to  ∞ (R), and they are 2-antiperiodic (i.e.,    ( + 2) = −   ()).Starting in the sequel, we will refer to the related functions {Ψ   } ,, which define the angular variation of the asymptotic displacement gradient: More explicitly, by the definition of the gradient operator in cylindrical coordinates, (2.16) Hence, the functions {Ψ   } ,, inherit the regularity and periodicity of {  } , .As presented in [17], these functions take the form  1 cos , with the constants depending on the elastic moduli.

Main results
In this section, we recapitulate the problem-specific interaction integral functionals presented in [17] and the problem which defines the approximate stress intensity factors.We then state the main results of the manuscript.We conclude with a discussion of the significance of the theorems.
Before defining the problem-specific interaction integrals, we introduce some notation.For a function  ∈  1 (ℱ), we may define its extension from the crack into the -neighborhood via the closest point projection:  ∘ .We abuse notation by writing both with the symbol .For the extension we have where  ′ denotes differentiation of  with respect to .
Next, given the basis we define two auxiliary fields needed for the interaction integral -the auxiliary gradient for mode where {Ψ   } ,, are in (2.16) and the field where  : (0, ∞) → R is any continuously differentiable cutoff function such that () = 1 for  ≤  0 <  and () = 0 for  ≥ .The combination  is referred to as the material variation, and defines how the material domain changes when the crack front is perturbed by the extension .
Definition 3.1.The problem-specific interaction integral s are the functionals where the five terms are For any two tensors   and   , we define where, 1 is the second-order identity tensor, and if We present the main properties of the problem-specific interaction integrals in the following theorem.
Theorem 3.2 (Properties of the problem-specific interaction integral).The following hold: (1) For any  ∈  1 (ℱ) and  ∈  2 (   ; R 3×3 ), and (2) Let  be the exact solution of Problem 2.1.If Assumption 2.2 holds, then for any  ∈  1 (ℱ), where the constants   are given in terms of the elastic moduli: An immediate consequence of the term-wise bounds in the previous theorem is the following. (3.17) We now define the approximate stress intensity factors { ℎ  }  , which belong to a finite-dimensional subspace K ℎ  ⊂  1 (ℱ).The parameter ℎ  denotes the discretization level of K ℎ  ; e.g., the number of basis functions scales like 1/ℎ  .We do not specify K ℎ  , though we request the following.There exists an integer  ≥ 2 such that, if   ∈   (ℱ), then for  independent of   and ℎ  inf which defines the order of approximation in K ℎ  .Meanwhile, for any  ∈ K ℎ  , we have the inverse inequality 4 In indicial notation, these are where, in a Cartesian basis, for some  independent of  and ℎ  , which is a consequence of the equivalence of norms in finite-dimensional spaces alongside a scaling argument.The approximate stress intensity factors are computed for a vector field  ℎ  ∈  1 (Ω  ; R 3 ) which approximates ; namely, we assume that where  may depend on .
With K ℎ  and  ℎ  , we are ready to define the approximate stress intensity factors.
Definition 3.4.The approximate stress intensity factor  ℎ  ∈ K ℎ  for mode  = I, II, III is the unique solution of the variational problem Finally, we have the following convergence result.
Theorem 3.5.Let  ℎ  ∈ K ℎ  solve (3.21) for mode  = I, II, III and any  ∈ K ℎ  .Then, where the constants  1 and  2 are independent of ℎ  and ℎ  but may depend on .
We conclude this section with the following remark, which highlights the key challenge overcome by the method proposed in [17].
Remark 3.6.From the term-wise continuity bounds in Theorem 3.2(1), we observe that ℐ  [, ] is linear and continuous with respect to  ∈  1 ( ) for any arbitrary  ∈  2 (   ; R 3×3 ).Meanwhile, when  = ∇, identity (3.15) implies continuity of ℐ  [, ∇] with respect to  ∈  2 ( ) only.As we will discuss in the proof of Theorem 3.2(2), cancellations occur within the interaction integral when  = ∇, notably the domain integrands form an exact divergence.However, when we seek approximate stress intensity factors, ∇ is unknown a priori, and we must use ∇ ℎ  instead of ∇.Hence, we lose the cancellations that enable continuity in  2 (ℱ), and we instead settle for continuity in  1 (ℱ).If instead the functional ℐ  [︀ , ∇ ℎ  ]︀ were continuous with respect to  ∈  2 (ℱ), we would no longer require the method proposed in [17].
Remark 3.7 (Periodic cracks).While the analysis in this paper is particularized to cracks with closed front ℱ, we may also consider configurations where the geometry and loading are both -periodic in , which allows us to treat the problem in a single period.Examples include the semi-infinite flat crack which is growing around a periodic array of obstacles (cf.[15], Fig. 2), or a semi-infinite crack with a helical perturbation to the crack front (cf.[23], Fig. 2).For these crack configurations, the functionals defined in Definition 3.1 and the approximate SIFs in Definition 3.4 are unchanged, but the domains of integration (ℱ,    , and  ±  ) are restricted to a single period in .
For the subsequent analysis to hold in the periodic case, in particular Theorem 3.2(2), there are additional periodicity conditions that would need to be imposed on ∇, the virtual extensions  and the function space K ℎ  , and the crack geometry.These conditions are discussed later in Remark 4.8.

Proof of the main results
We prove Theorems 3.2 and 3.5.To proceed in certain locations, we introduce additional results; proof of these may be found in the appendix.

Properties of the problem-specific interaction integrals
We begin by defining a tensor space which is important for the definition of the interaction integral.Close to the crack front, tensors in this space behave like the asymptotic displacement gradient (2.15), rotated into the local basis { 1 ,  2 ,  3 }.Definition 4.1 (Space of asymptotic displacement gradients).Let where {Ψ   } ,, were introduced in (2.16).For any tensor  ∈ B  , {  }  are the stress intensity factors of .
We remark that B  ⊂  2 (   ; R 3×3 ).We next present a regularity result for the auxiliary gradient fields { aux, }  , namely that they are the direct sum of a tensor in B  and an  1 tensor field.The proof can be found in Appendix B. Proposition 4.2 (Regularity of  aux, ).The tensor field  aux, ∈ B  ⊕  1 (   ; R 3×3 ).Moreover, We are now ready to prove Theorem 3.2(1).
Second, equipped with the previous lemma, we show that where  is defined in (3.4).We now turn to the proof of Lemma 4. A similar identity is given by Gosz and Moran [16], and may be verified through direct calculation using the exact expressions for {Ψ   } ,, .Equation (4.4) amounts to an orthogonality result for the three modes of {Ψ   } ,, .The goal of our proof of Lemma 4.3 is to transform the left-hand-side of (4.4) to the right-hand-side of (4.2).We break this down into a number of steps.First, we use the integral over   () to obtain a test for virtual crack extensions by integrating over   .Proposition 4.5.Under the assumptions of Lemma 4.3, where  is extended to   via the closest-point projection, cf.(3.1).
Proof.Let ℎ be the stretch factor defined in (2.9).Because  ,  ∈ B  (and hence Here, the  −1 of Σ(   ,    ) has been exactly canceled by the length element d =  d.It thus follows that we may modify (4.4) to get that lim Integrating both sides over the crack front, we have Because of (4.6), we may apply the Lebesgue Dominated Convergence Theorem and Fubini's Theorem to the left-hand-side of (4.7), pulling the limit outside of the integral over ℱ and combining the double integrals, respectively, to yield the conclusion.
We next enlarge the set of tensor fields to which an expression like (4.5) is applicable from ).We require the following result for the trace of an  1 function over the boundary of a shrinking neighborhood; its proof may be found in Appendix C. While the result is stated for scalar-valued functions, it trivially holds for vector-and tensor-valued functions such as   ∈  1 (   ; R 3×3 ).Let us consider the term with IJ = TS.Fix 0 <  < .Then For the first inequality, we used the fact that |  | ∈  ∞ (  ).For the second inequality, as in the proof of Theorem 3.2(1), there exists a constant  such that |Σ(   ,    )| ≤ |   ||   |, which we lumped with the term outside of the integral.Using the exact form of    , we know that ‖   ‖ 0,  ≤ , where the constant depends only on {   }  and is independent of .Taking the limit as  → 0, we get The TS case of (4.10) follows from Lemma 4.6.Analysis of the ST and SS terms is handled similarly.
We are now ready to complete the proof of Lemma 4.3 Proof of Lemma 4.3.For notational convenience in this proof, we write Fix 0 <  < .We introduce the following domain:   ∖ (  ∪ ).This is a cut, hollow neighborhood of ℱ with four boundary surfaces: the outer wall   , the inner wall   , and the positive and negative sides of the crack ( ∩ (  ∖   )) ± .In this domain, we may apply the divergence theorem and we recover the integrand of term (3.6).Next, let us expand the divergence in the volumetric integral of (4.2): div We immediately recognize the first and third terms as the integrands of (3.8) and (3.10), respectively.Lastly, for the second term, we compute the divergence of Σ(∇,  aux, ), apply the compatibility of ∇, and equate −div((∇)) with  to yield These are the integrands of (3.9) and (3.7), respectively, and we reach the conclusion.
When we apply the divergence theorem in (4.11), we must consider integration over the end caps of    at  = 0 and  = , which we term    (0) and    (), respectively.By the periodicity assumption on , we have that the limits of integration of (, ) on    (0) and    () are identical.Meanwhile, the assumed periodicity of ,   ,   , and   cause  (0, , ) to hold for any (, ).The outward normals to    (0) and    () are − 3 (0) and  3 (), respectively, and hence the two surface integrals precisely cancel, which leaves the conclusion of Lemma 4.3 unchanged.
To apply this result to Theorem 3.2(2), we need the above assumptions to hold on ∇,  aux, , and  from (3.4).By the definitions of  aux, , and  in (3.3) and (3.4), respectively, the first and second assumptions are satisfied if holds for any  and for any ,  = 1, 2, 3.
As an example, these conditions allow us to consider cracks on a torus with  -fold azimuthal symmetry by integrating over only 1/ of the tubular neighborhood.Remark 4.9.In Assumption 2.2, we assumed that the stress intensity factors of  belonged to  2 (ℱ).Consequently, it was possible to show that ∇ ∈ B  ⊕  1 (   ; R 3×3 ).However, the space B  only requires stress intensity factors belonging to  1 (ℱ).It is possible to relax the restrictions on the input tensors, i.e., enlarge the space B  ⊕  1 (   ; R 3×3 ).For clarity, we opted not to do so in this paper.Nonetheless, a possible set of sufficient conditions for the enlarged space are as follows.First, to ensure div (thereby making valid (4.2)), we require tensors   and   such that ∇ , − (∇ , )  and div(( , )) are both square integrable on    .Second, so that they do not contribute to the final value of the general interaction integral, the parts of   and   not belonging to B  must satisfy (4.8).If the stress intensity factors of  belong only to  1 (ℱ), the first condition is still satisfied by ∇ through symmetry of second distributional derivatives and because  solves Problem 2.1.Via direct calculation, the second condition is also satisfied by the part of ∇ not belonging to B  ⊕  1 (   ; R 3×3 ) (which possesses  1/2 -dependency around ℱ).

Convergence of the approximate stress intensity factors
We conclude this section with a proof of Theorem 3.5.
Proof of Theorem 3.5.
For the first term, by the interpolation estimate (3.18) For the second term, we have where we have used the definition of the  2 -projection, Theorem 3.2(2), and the definition of  ℎ  (3.21).Via Corollary 3.3 Application of the inverse inequality (3.19) gives Dividing through by ‖ ℎ    −  ℎ  ‖ 0,ℱ , we get Finally, using the convergence estimate (3.20) for ‖ −  ℎ  ‖ 1,Ω  , we reach the conclusion.

Numerical study of error estimates
Here, we assess the error estimate in Theorem 3.5 through a numerical example.As shown in the proof of that theorem (see (4.12), combined with (4.13) and (4.14)), the error in the stress intensity factors is bounded as where  ℎ  :=  ℎ    −  ℎ  .As stated earlier, the two terms represent interpolation errors in the crack front functions space K ℎ  , and consistency error in our definition of  ℎ  .The functionals studied in this manuscript have no effect on the interpolation error; hence, in this section we will focus only on the consistency error, which for ease of notation we write as Err ()   .In [17], we report the full error (3.22) for several numerical examples.On a computer, ℐ  [, ∇ ℎ  ] is not computed exactly; rather, we use quadrature, which results in an operator ℐ  [, ∇ ℎ  ].In [17], we proposed applying the same quadrature rules on the finite element mesh that were used for computing  ℎ  .If we take into account the quadrature-evaluated interaction integral when defining  ℎ  , then the second term in the above error estimate becomes Since this term incorporates both inconsistency and quadrature errors, we use the shorthand Err (+)

𝛼
. By subtracting and adding ℐ  [ ℎ  , ∇ ℎ  ] in the numerator, we may partition the prior term in two, separating consistency and quadrature errors.
In the proof of Theorem 3.5, we bounded the consistency error using Corollary 3.3, the finite element error (3.20), and an inverse inequality (3.19) where we have used the fact that for standard FEM (e.g., [18]), the finite element error converged with order ℎ 1/2 .The suboptimal convergence rate results because the radial singularity in the exact elastic solution (2.11) is not well-approximated by piecewise polynomial basis functions [29].Equivalently, the radial singularity reduces the regularity of , causing the function to belong to a space such as  3/2 (Ω  ), which possesses suboptimal convergence of interpolation errors [14].
Estimating the quadrature error is more challenging.Here, standard quadrature estimates like Theorem 8.5 of [13] no longer apply, as the integrands of (3.6)-(3.10)contain radial singularities from the auxiliary field  aux, .Furthermore, if the function space K ℎ  is one with low continuity, for example the 1-D  1 Lagrange finite element space constructed over a mesh of the crack front, then the extension of test functions  into    also reduces the regularity of the integrand, thereby affecting quadrature convergence rates.These effects may be taken into account when constructing quadrature error estimates (e.g., for a function with a radial point singularity, see [25]).Rather than perform such (long) computations, we instead assess the quadrature error alongside the consistency error through a numerical example.

Example
In this section, we explore numerically the consistency and quadrature errors using an example adapted from [17].We consider the semi-infinite crack with straight crack front, subjected to 1-periodic loading along the  3 -direction.The geometry and loading for this problem are periodic in  =  3 , and we may verify that the conditions outlined in Remark 3.7 hold, so that our analysis is still valid for this problem.The crack front is straight and the crack surface is flat, hence for this geometry we have   (, ) =   () =   .Examples of non-planar cracks and cracks with curved fronts may be found in [17].
We construct an analytical solution with non-uniform mode II stress intensity factor by imposing appropriate tractions and body forces.The analytical displacement field around the crack front is The stress intensity factor  II is an even, 1-periodic function in  3 , taking value for  3 ∈ [0, 1].For this displacement field, the crack faces are traction-free and the required body force is In this way,  is the solution of Problem 2.1 with body force  and crack-face tractions  ≡ 0 on   Ω  =  ± .

Computation
We restricted our attention to the finite domain Ω = (−0.3,0.3) × (−0.3, 0.3) × (0, 1).We prescribed periodic boundary conditions on the faces with  3 ∈ {0, 1}, and Dirichlet boundary conditions on the faces with | 1 | = 0.3 and | 2 | = 0.3.As previously stated, the crack faces were traction free.We discretized the domain with a family of unstructured, tetrahedal meshes found through successive subdivision of each tetrahedron into eight smaller ones [24].The problem geometry and the coarsest mesh are shown in Figure 2.
For this example, the stress intensity factors were approximated using trigonometric polynomials with maximum order   ∈ {5, 10, 20, 40, 80}.For further discussion of these spaces, see [17].Unlike functions in the  1 Lagrange finite element space, trigonometric polynomials are smooth, thereby eliminating a potential source of quadrature error.Top row: variation in the errors under mesh refinement, fixing the maximum order of the spectral basis.Bottom row: variation in the errors with respect to the order of the spectral basis, fixing the bulk mesh refinement.Dashed lines indicate the computed errors using a four point quadrature rule in each tetrahedron, while solid lines show the errors when using 2048point quadrature (see text).Compared with the estimate in (5.3), the errors are superconvergent with respect to bulk mesh size and more slowly growing with spectral basis order.
Numerical integration was performed using a standard, second-order quadrature rule with four points in each tetrahedron (e.g., [30]).To isolate the effect of consistency error, we also used a 2048-point rule, which was found by subdividing each tetrahedron into eight smaller ones [24] three times, and applying the basic four point rule in each subdivision.With such a large number of quadrature points, we expected Err (+)  ≈ Err ()  .

Results
In Figure 3, we plot the Err (+)

𝛼
for the each stress intensity mode.For fixed   , the errors in the three terms converged with approximate order ℎ  , rather than the expected ℎ 1/2  from (5.3).The behavior of the error with respect to   was more complex.At fine values of ℎ  , the error grew more slowly than  1  .Rapid increase in the error was observed whenever   (ℎ  /ℎ 0 ) > 5.As seen in Figure 2, the coarsest mesh had ℎ 0 ≈ |ℱ|/10.Hence, the condition   ℎ  > |ℱ|/2 corresponded to cases where the highest wavenumber basis functions were poorly sampled on the given mesh.In practice, one would avoid such behavior, by ensuring that the highest frequencies in the crack front basis were adequately resolved on the bulk mesh.
We remark on quadrature.For each mode, quadrature errors became dominant for coarse bulk mesh size ℎ  and for large spectral basis order   .As expected in these cases, standard quadrature rules were insufficient to resolve both the radial singularities in the auxiliary fields and the rapid variation of the high wavenumber basis functions along ℱ.However, compared with consistency errors, quadrature errors converged more quickly with respect to ℎ  .

Superconvergence
We lastly comment on the appearance of superconvergence in the consistency error.As shown in Theorem 3.2(1), the problem-specific interaction integral is a continuous affine functional of the displacement field.It is well-known that continuous linear functionals applied to finite element solutions converge faster than expected from continuity (e.g., [3,4]).
The key step in classical estimates is to treat the functional (i.e., ℐ  [, ∇•] :  1 (Ω  ; R 3 ) → R) as the data for an adjoint problem (although, in our present case, the linear elasticity operator is self-adjoint).Then the functional error may be written in terms of the dual problem's solution   [] and its best approximation in the finite element function space V ℎ  : Note that the previous expression assumes that  −  ℎ  coincide on the Dirichlet boundary.If we assume that   [] possesses the same regularity as , then we get the so-called "rate doubling" behavior common in the finite element literature.The errors observed in Figure 3 may be indicative of such behavior.
It is trivial to show that   [] ∈  1 (Ω  ; R 3 ); hence the best approximation error will converge [13], and the functional error will be superconvergent.However, it is unknown a priori whether   [] possesses further regularity.The answer to this open question will rely on the properties of the functional ℐ  [, ∇•] and geometric features of the problem domain Ω  .Finally, in order to predict how changing the crack front basis K ℎ  affects errors, results for the regularity of   [] and its best approximation must also quantify the dependence on , a non-trivial task that requires further research.

Conclusion
In this work, we presented analysis of a method to approximate the stress intensity factors along the front of a three-dimensional crack.In particular, we proved that the functionals used in the method have two important properties, namely (a) that when applied to the exact displacement gradient, we recover a weighted integral of the stress intensity factors over the crack front, and (b) that the functionals are continuous.The latter property is essential for proving convergence of the method: for fixed , the functional error is guaranteed to converge if the gradient of the finite element solution converges.
We then presented the error analysis for the numerical stress intensity factors.We showed that this error was bounded by two terms.The first term corresponded to an interpolation estimate of the stress intensity factors in the finite-dimensional function space we construct over the crack front.The second was a consistency error in the interaction integral, which we estimated via the continuity results of Theorem 3.2(1) and Corollary 3.3.
A numerical example provided insights beyond the error estimate of Theorem 3.5.First, quadrature errors, a practical consideration in the implementation of the method, were faster converging than the consistency error in the case where  ℎ  was computed with standard finite elements (e.g., [18]).We caution the applicability of this observation to higher-order bulk numerical schemes (e.g., XFEM with singular tip enrichment or the Mapped FEM).Second, the consistency errors demonstrated superconvergent behavior, with a possible explanation being that the dual problem solution   [] possessed similar regularity to .Further analysis is needed to make more concrete statements on either point; however, we believe these two issues are exciting and challenging prospects for future work, as they may be crucial ingredients for ensuring rapid convergence of the numerical stress intensity factors in higher-order schemes.
divergence operator, which are used in Appendix B to prove Lemma 4.2.Additionally, for Appendix B, we prove a regularity result for the inclination angle  defined in (2.5).

A.2. Integration in tubular coordinates
In the proof of Theorem 3.2 and Appendix C, we perform integration over    ,   , and   in the tubular coordinate system.Here, we derive the expressions for the relevant Jacobians.

A.3. Differentiation in tubular coordinates
Because {  ,   ,   } is an orthonormal basis, for any differentiable function  , we wish to write where    is the directional derivative of  in the direction of .Computing the directional derivatives, we can show Note that this relationship holds for vector or tensor functions as well: for example, if  is a vector (or tensor) field, then

A.4. Curvilinear coordinate system
In certain situations in Appendix B, it is convenient to introduce the coordinate system { 1 ,  2 ,  3 } for the neighborhood   as in [16].Here, we set and define the map This is a diffeomorphism between Ξ and   .The curvilinear coordinates are related to the tubular coordinates: A.4.1.Differentiation in the curvilinear coordinate system The gradient operator in the { 1 ,  2 ,  3 } coordinate system, expressed in the { 1 ,  2 ,  3 } basis, is Later in Appendix B, when applying the above formula, we will rewrite the  3 -component as which makes use of the chain rule: A.4.2.Divergence in the curvilinear coordinate system In Appendix B, we also require an expression for the divergence of a tensor field, to which we build up starting from the divergence of a vector field.Let  be a vector field Then div() = tr(∇).Using the previous expression for the gradient, we may derive show Next, let  be a tensor field and  is any constant vector: Then, div() ]︃ . (A.4) The coefficients of  •   are the   -components of div().Again, when applying this expression, we will use the chain rule to rewrite

Appendix B. Properties of the auxiliary gradient fields
In this section, we prove Lemma 4.2, using the following results.
Proposition B.1.For each  = I, II, III, let Proof of Propositions B.1 and B.2 are performed via direct calculation.
For the first term, we expand inside the brackets: Thus, there exists  such that  Appendix C. Trace of an  1 function on a shrinking neighborhood This section is devoted to the proof of Lemma 4.6.For now we consider the case where  is scalar-valued, though extension to vector-or tensor-valued functions is trivial.Our proof relies on two ingredients: the construction of a diffeomorphism between    and a reference domain which is independent of , and a suitable trace inequality posed in the reference neighborhood.For a fixed , the combination of mapping to the reference domain, applying the trace inequality, and mapping back to the physical domain will allow us to bound the above surface integral by a positive power of , and the result follows by letting  → 0.
The regularity of , X, and   guarantees that   is a diffeomorphism between    and N   .By design of   , integrals over   map to integrals over   .Meanwhile, the map   stretches within orthogonal sections of   ; we next derive a trace inequality which depends only on the radial derivatives of the function f =  ∘   .This result follows from a similar argument to Brenner and Scott ([6], Sect.1.6).Proof.For ease of writing this proof, we will drop the "hat" from all symbols (e.g., we will write  instead of f ).We begin by splitting    in two.If Θ + = Θ ∩ {(, , ) :  > 0} (with similar definition for Θ − ), we let  ±  = (Θ ± ).In other words, we extend the crack across the ligament in the  1 -direction.
We now proceed for the "positive" half neighborhood.Proof for the "negative" half is nearly identical.Let  ∈  1 ( +  ).Then, By density, we extend the above result to  + ∈  1 ( +  ).Summing the results for  + and  − yields the conclusion.
Equipped with the mapping   and the above trace inequality, we are now ready to prove Lemma 4.6.
Proof of Lemma 4.6.Fix .We map to the reference neighborhood, 5 Under   , area elements of   change according to ℎ ĥ   . 6The Jacobian of  −1  is given by ĥ ℎ  2 2 , while one may show | f,r For the second term on the right-hand-side of (C.Letting  → 0 yields the desired conclusion.

Figure 2 .
Figure 2. Top: geometry for a straight, flat crack, crack geometry.Bottom: perspective view showing the coarsest mesh in the family of unstructured tetrahedral meshes used for the convergence test.

Figure 3 .
Figure 3. Consistency and quadrature errors (5.2) for the problem-specific interaction integral.Top row: variation in the errors under mesh refinement, fixing the maximum order of the spectral basis.Bottom row: variation in the errors with respect to the order of the spectral basis, fixing the bulk mesh refinement.Dashed lines indicate the computed errors using a four point quadrature rule in each tetrahedron, while solid lines show the errors when using 2048point quadrature (see text).Compared with the estimate in (5.3), the errors are superconvergent with respect to bulk mesh size and more slowly growing with spectral basis order.