A quasi‐static discontinuous Galerkin configurational force crack propagation method for brittle materials

This paper presents a framework for r‐adaptive quasi‐static configurational force (CF) brittle crack propagation, cast within a discontinuous Galerkin (DG) symmetric interior penalty (SIPG) finite element scheme. Cracks are propagated in discrete steps, with a staggered algorithm, along element interfaces, which align themselves with the predicted crack propagation direction. The key novelty of the work is the exploitation of the DG face stiffness terms existing along element interfaces to propagate a crack in a mesh‐independent r‐adaptive quasi‐static fashion, driven by the CF at the crack tip. This adds no new degrees of freedom to the data structure. Additionally, as DG methods have element‐specific degrees of freedom, a geometry‐driven p‐adaptive algorithm is also easily included allowing for more accurate solutions of the CF on a moving crack front. Further, for nondeterminant systems, we introduce an average boundary condition that restrains rigid body motion leading to a determinant system. To the authors' knowledge, this is the first time that such a boundary condition has been described. The proposed formulation is validated against single and multiple crack problems with single‐ and mixed‐mode cracks, demonstrating the predictive capabilities of the method.


INTRODUCTION
Recently, there has been significant interest in the numerical prediction of crack propagation. However, despite numerous frameworks being proposed, accurate and efficient simulation of crack probation is still one of the most challenging problems in solid mechanics. The fracture mechanics community requires algorithms that can predict the evolution of cracks from initiation through to large-scale propagation. In this paper, we present an algorithm based on the concept of configurational forces (CFs) combined with a discontinuous Galerkin (DG) numerical framework that allows for efficient brittle crack propagation in 2 dimensions.
The work of Eshelby, 1,2 Rice, 3 and Irwin 4,5 are fundamental to all work in this field. The local variational formulations in previous studies [6][7][8][9][10][11][12] use a CF acting at a crack tip to describe the propagation of a crack. The CF values have been determined numerically at static fracture fronts by other studies. [13][14][15][16] Using the CF to describe a moving fracture front was initially attempted by Mueller and Maugin 17 within the conventional finite-element context and Larsson and

Continuous formulation in space and time
Consider the homogeneous evolving material body in Figure 1, ℬ⊂ ℛ 3 , which has an exterior boundary ℬ, crack lips Γ, and crack tip Γ. The subset ℬ Γ = ℬ∖Γ ∪ Γ is defined with the crack lips Γ + and Γ − and the surface C encircling the crack front. In the limit Γ + → Γ, Γ − → Γ, and C → Γ. Finally ℬ Γ = ℬ ∪ Γ − ∪ Γ + ∪ C is a set containing all the boundaries of the problem. Here, ℬ= ℬ D ∪ℬ N ∪ℬ T where ℬ D , ℬ N , and ℬ T are the Dirichlet, Neumann, and roller boundary conditions, respectively.
The material points are defined x∈ ℬ Γ . They evolve from a set of reference coordinates ∈ Ω, where Ω⊂ ℛ 3 is the reference domain. This is achieved by using the time-dependent mapping t ( ) = x, which allows for a change in material structure that reflects the propagation of a crack. The small strain displacement at time t∈ ℛ + is defined as Next, we define the global power postulate of the mechanical form of the second law of thermodynamics for power dissipation, , is the power applied to the boundary ℬ by a traction t,̂is the free energy function, and is the infinitesimal strain Following the work by Miehe et al, 21 (2 1 ) becomes where =̂( ) is Cauchy stress, = ( ) −( u t ∕ x) T is the Eshelby stress, and is an identity matrix. The boundary conditions for the spacial velocity field v are which has a prescribed valuev on a Dirichlet boundary ℬ D . The boundary conditions for the material velocity field, V, are with a material velocityȯ at the crack tip. Given that v and V are arbitrary in ℬ Γ and have boundary conditions (5) and (6), respectively, the following statements of equilibrium can be defined · = 0 in ℬ Γ , · n = t on ℬ, · n = 0 on Γ, and · = 0 in ℬ Γ .

Discrete formulation in space
A discrete formulation of the power dissipation also exists. This is taken directly from (4) using isoparametric shape functions N I , which act on a node I of element K. K ∈  , where  is a subdivision of the polygonal domain ℬ Γ ⊂ R 2 into disjoint triangular elements, with its coordinates mapped into the reference domain using −1 t (x). The nodal material and spacial velocities existing on element K are respectivelyḊ I andḋ I . The derivative of the shape functions in the material domain B I (x) = ∇ x N I also exists. Equation 4 can be discretised into the form where n I is the set of all element nodes in the mesh. The conventional force components in (8) are defined as and the CF is A is the usual finite element summation operator. Additionally, the material and spacial velocities in element K take the discrete form, denoted by the superscript h, with all interpolation occurring in the material domain, wherė is a prescribed displacement on the boundary and, whereȯ I is the crack tip material velocity. A consequence of (12) is that the following is true This means that the reduced discretised global power dissipation (8) at the crack tip is g I is the CF ; it determines crack growth and direction and is calculated in post processing procedure, once the linear elastic system has been solved.

Discrete formulation in time
The final step is determining how the crack will propagate. Here, we will use a quasi-static crack propagation framework as presented in Miehe et al, 21 Miehe and Gürses 22 to perform a quasi static analysis. First, we integrate the discrete dissipation power at the crack, (15), over the time period [t n , t n+1 ] This gives an incremental constant increase in the crack surface length, Δo I , over the time period [t n , t n+1 ]. It has the form where g c is a Griffith material failure criteria. h o is the increase in crack length defined as h o = g I ∕|g I | · L c , L c is the original length of the most aligned element face with g I ∕|g I |. A crack will propagate the entire reorientated element face associated with L c . Δ I is subject to the Karush-Kuhn-Tucker conditions Motion of nodes can be permitted in the material configuration, except motion that would change the shape of the boundary. We recognise that is possible to dissipate power by moving nodes in the material configuration 60,61 and thus achieve a minimal energy solution to the problem by increasing the total power dissipated by the term ∑ n I I=1 g I ·Ḋ I in (8). However this is a highly non-linear problem and therefore computationally expensive. We therefore do not solve for it here, consistent with the works of previous studies 16,[21][22][23]42,43,45 and many others but instead recognise that it could potentially improve our solutions. Here, we only consider power dissipation in the form of surface generation, or crack propagation, when the Griffith failure criterion |g I | > g c is satisfied at a crack tip.
The key equations for modelling brittle fracture propagation based on CF have now been outlined. Their values are calculated in a post-processing procedure once the linear elastic system for small strain problems has been solved. The linear elastic system is cast within a DG framework presented in the next section. It should also be stated that it is possible to simultaneously solve for the CF and material velocity as in Kaczmarczyk et al, 24 based on the works of Kuhl et al. 62 However, there is debate in the literature on the validity of linearising the CF in these approaches. 63

DISCONTINUOUS GALERKIN FINITE ELEMENT APPROXIMATION
We consider the following model problem on a bounded Lipschitz polygonal material domain ℬ Γ in R 2 , with ℬ N ∪ ℬ D ∪ ℬ T = ℬ Γ . In this section, three different boundary conditions are described within the DG formulation. On ℬ D , Dirichlet boundary conditions are applied; the prescribed displacement on this boundary has a value g D . On ℬ N , Neumann boundary conditions are applied; here, the traction has a value g N . Last, a roller boundary condition is applied on ℬ T ; here, the tangential component to the boundary surface has a value of zero, and the displacement normal to the surface has the value g T · n, where n is the unit vector normal to the boundary. The strong form of the problem with the boundary conditions is defined as where n || is the tangential unit vector to the boundary. The prescribed values on the respective boundaries, g D , g T , and g N , are the prescribed displacements (fully prescribed g D and roller g T ) and the traction vector g N . Each element K of the mesh  , where  is in general irregular, is the image of the reference triangle under an affine elemental mapping F K ∶K → K. We denote by  (K) the set of the 3 elemental faces of an element K. If the intersection F = K ∩ K ′ of 2 elements K, K ′ ∈  is a segment, we call F an interior face of  . The set of all interior faces is denoted by  I ( ). Analogously, if the intersection F = K ∩ ℬ Γ of an element K ∈  and ℬ Γ is a segment, we call F a boundary face of  . The set of all boundary faces of  is denoted by  B ( ), and it is the union of the 3 sets  N ( ),  D ( ), and  T ( ) of faces on the 3 boundaries ℬ N , ℬ D , and ℬ T , respectively. Additionally the internal crack is represented by tractions on K ∩ ℬ N ∩ (Γ + ∪ Γ − ) as zero and the crack tip being represented by an element node existing at Γ. Moreover, we set  ( ) =  I ( ) ∪  B ( ). For each element K ∈  , we define p K to be the order of the element. We also define the vector For any mesh  of ℬ Γ with the degree vector p, we then define the hp-version DG finite element space by represents the functional space of the 2 components of the function w and  p K (K) denotes the set of all polynomials on the triangle K of degree no more than p K . The basis functions chosen are hierarchical, 64 with the test function and displacement described respectively as w = N s w s and u h = Nu s . w s , and u s are the hierarchical basis function coefficients, and N s are the hierarchical functions in the reference elementK.
We now introduce the SIPG method in the bilinear form, for the approximation of the model problem (19): for all w ∈ W p ( ), where where is a penalty term for linear elastic SIPG In this paper, = 10, its range is defined by Hansbo and Larson, 59 E Y the material's Young's modulus, and as Poisson's ratio. Further, h F is the element face length and where {·}, [[·]], (·, ·) and ⟨·, ·⟩ are defined in Arnold et al. 65 The first term on the right hand side of (21) describes the virtual energy in the material bulk. The second is a face stiffness term, which averages the jumps in tractions existing between elements. This is followed by its transpose to make the global stiffness matrix symmetric. The fourth term is required to stabilise the method. The last 3 terms impose the conditions on the boundary ℬ T weakly, the first of the 3 terms is the weak implementation of the boundary, and the next is the transpose followed again by a term to stabilise the method.
The right hand side of (20) has the form The first term is the weak implementation of the Dirichlet boundary condition, followed by its stabilising term. The third term is the implementation of the Neumann boundary condition. The fourth is the boundary condition applied on ℬ T followed by its stabilising term. For more information on the implementation of problem (20), we invite the reader to refer to Hansbo and Larson. 59

NUMERICAL CALCULATION OF CONFIGURATIONAL FORCE
Here, we consider the tip 21,22 and domain methods 14 for calculating the CF at the crack tip within the SIPG framework. The tip and domain methods are shown in Figure 2A and 2B, respectively. The tip method considers only the CF value at the crack tip node, marked as white in Figure 2A. The choice of this node is a result of the power dissipated from crack growth being only associated with CF at this node, as discussed in Section 2.3.
The domain method, shown in Figure 2B, considers the CF on the interior set of nodes, marked as white, of elements within the boundary r d . The motivation for this choice is that the CF values should be zero on the interior set of nodes, except the crack tip node at Γ. Denzer et al 14 concluded that spurious CF values exist at interior nodes around the crack tip, which should exist at the crack tip. It was found in summing these to the value at the crack tip gave a more accurate solution for the CF. As shown by Irwin, 5 the stresses close to the crack tip are a function of r −1/2 , where r is a distance from the crack tip. The finite elements within the discretised space struggle to capture the stress singularity, and so the CF, a function of stress, is poorly represented. When using elements, which can capture the stress singularities at the tip, 66   The equation for calculating the CF for the tip and domain method is For the domain, method A is the set of elements K, which have all their nodes within r d . For the tip method, the set A is the set of all elements K, which have a node on Γ. Finally, n b is a list of all nodes in A, which do not lie on the exterior. The white shaded nodes in Figure 2A and 2B are the interior nodes and crack tip node in the set n b .
Evaluation of (25) is performed in a post-processing procedure after (21) has been solved.

RP-ADAPTIVITY
In this section, we describe a method for propagating a crack in quasi-static rp-adaptive procedure using CF fracture mechanics cast within the presented SIPG finite element framework using hierarchical shape functions. The benefit to using SIPG is the flexibility available to switch off face interactions between elements by removing the DG face stiffness terms from the global stiffness matrix. This creates new surfaces and is used to propagate a crack. No dof are added to the data structure to propagate a crack whilst only minimal manipulation is required to enable a p-adaptive scheme. The data structure is arranged such that all the dof corresponding to first-order components of all elements are numbered first. The labelling of all these dof is unchanged throughout a simulation, as this is the minimum requirement for a finite element discretisation to exist. All subsequent higher order dof are numbered greater than their first-order counterparts. An example of the data structure is shown in Figure 3. The numbers on the rows and columns of the matrices correspond to element numbers in the mesh

Hierarchical shape functions
In this paper, we use shape functions defined for a referenceK triangular element in Solin et al 64 with improved interior bubble functions. An advantage of the adopted shape functions is that values for an arbitrary high-order shape function can be created in a hierarchical algorithm. This means that it is not necessary to hard code equations for high-order shape functions and an element of an arbitrary high order can be generated.

rp-Adaptivity algorithm
In a crack propagation scheme, the CF, g I , is evaluated at each crack tip using either the tip or domain method. If |g I | ≥ g c , then the crack will propagate in the direction g I and the rp-adaptivity method will be applied as given in Algorithm 1.
An example of a crack propagating through a mesh, using Algorithm 1, with its corresponding changing global stiffness matrix is shown in Figure 3. The mesh is constructed from 6 elements. For the simplicity of this example, only elements sharing a node at the crack tip having a polynomial order, p K , greater than 1. It is possible to have a group of elements with p K > 1 about the crack tip; these elements reside within the radius r p . An element is considered inside r p if at least one of its nodes are inside r p . To propagate a crack, first, the linear elastic system is solved producing a stress field. g I is then calculated from the stress field in the material domain, then following This is equivalent to applying homogeneous Neumann boundary conditions on the new crack surfaces. Furthermore, the DG face stiffness matrix calculations for this face are also removed from any further calculations to prevent any face interaction reappearing. This removes any direct interaction between elements along the face creating a new surface, which extends the boundary of the domain, and propagating the crack.
• Step 3, as only elements on the crack tip have a polynomial order greater than 1, and the crack has moved, all rows and columns of the global stiffness matrix associated with the higher order dof of elements no longer at the crack tip are removed.
This is highlighted by the black lines through the final rows and columns of the second stiffness matrix. Additionally, as the geometry of elements, which share a node with the new crack tip, have also changed, all values associated with these elements' local stiffness are removed. This is represented by the solid black blocks in the second matrix.
• Step 4, the updated local stiffness matrix components of elements are added back into the matrix. This corresponds to element with a changed geometry or increase in polynomial order. All new values are highlighted with black boxes in the third matrix in Figure 3.
The specific detail of the rp-adaptive method, which Figure 3 follows, is provided in Algorithm 1. The last stage of Algorithm 1 is recalculating the DG area and surface local stiffness matrices for E r and E p and adding them back into the global stiffness matrix. When adding new higher order dof, the only constraint is that their numbering is unique to an element, as with all dof in DG.

Average boundary conditions
This type of boundary condition is used for nondeterminant 2D systems. This is the case for linear elastic problems where the displacement is not constrained in both axes. To the author's knowledge that this is the first of its kind implemented for a linear elastic problem. These boundary conditions average the displacements in x and y (u and v), and the rotation, to be 0 and Here, u K and v K are dof corresponding to u and v within element .
Equations 26 and 27 are incorporated into the global stiffness matrix and thus form part of the solution.
where ne is the number of elements, f n is l(w) evaluated at each node, U contains all displacement dof and BC is a set of arbitrary unknowns that form part of the solution vector.

NUMERICAL EXAMPLES
In this section, several crack propagation problems are investigated to demonstrate the accuracy of the SIPG for predicting crack paths against analytical numerical benchmark tests and experimental results. All problems considered in this section assume linear elastic material behaviour. Only small strains in plane strain are present with the free energy where is strain and C is the fourth-order stiffness matrix for linear isotropic elasticity. The accuracy of the numerical scheme for predicting a crack propagation path can be evaluated by considering the CF value and direction for mode 1, mode 2, and mixed mode problems.

Single edge notched static tensile test
This static single edge notched (SEN) test is used to show that SIPG method produces CF values within the range of accuracy obtained in literature. 14,21 The value of the CF is dependent on • the characteristic mesh size at the crack tip h cF ; • the domain around the crack tip where elements have p K > 1 defined by the radius r p ; and • the radius containing all elements used in the domain evaluation of the CF r d , (25).
Additionally, r h defines the region around the crack tip where elements are of a different length scale to the rest of the mesh. These variables are also shown graphically in Figure 4A.
The geometry of the test is taken from the benchmark provided by Fagerström and Larsson 19 and is shown in Figure 4A. Here, the crack length a = 0.1 m, the width of the plate b = 0.5 m, the half height of the plate H = 1.0 m, and the uniaxial tensile stress applied is = 10 MPa. The plate has a Poisson's ratio of = 0.3 and a shear modulus = 80 GPa. Zero average displacement and rotations boundary conditions were applied using (26) and (27). Last, the mesh was generated using the unstructured mesh triangle generator, Triangle, 67 and is shown in Figure 4B. All h-refinement occurred in a homogeneous manner uniformly across the entire mesh, as in Figure 5.
An analytical solution of the energy release rate for the plane stress problem, in Figure 4A, for the mode 1 crack path is derived from the J-integral, where the empirically corrected SIF, K I , can be determined from multiple empirical equations provided by Tada et al. 68 So not to be biased, we choose to define a range in the SIF and therefore also the expected energy release rate. The range is provided by the empirical equations, which give the smallest and largest K I values in Tada et al, 68 respectively:  and Using the SIF, K o = √ a for an infinite plate with a crack length of 2a. 5 The energy release rate for the mode 1 SEN problem is in the range G ∈ [256.7, 261.8] N/m. As the finite element solutions converge from below for h and p refinements, we will consider G = 261.8 N/m as our reference solution and comment all percentage differences from this value, unless stated otherwise. Considering virtual work in the direction of mode crack propagation, the reference mode 1 CF value is g = 261.8 N.

SEN static tensile test results
The first investigation was conducted to verify that using SIPG, it was possible to obtain accuracies within the range obtained in literature for CF using CG methods, see Table 1. Given that the most accurate results obtained in both Denzer et al and Miehe et al 14,21 were obtained using the domain method, this method was used to compute the CF for the problem described in Figure 4A, Section 6.1.
First, we investigated how varying r d around the crack tip affected the accuracy of the CF. For this test, r h = 0.1 m and p K = 1 for all elements. The initial element length within r h was set to 0.04 m and graded to 0.34 m outside r h . The mesh refinement was uniform, as in Figure 5.
The changes in CF magnitude when considering element rings to compute the domain evaluation of the CF (25), Figure 6 is consistent with Denzer et al and other works. 14,[69][70][71][72] We demonstrate that the CF domain integral 14 is poorly represented by all refinement levels when only considering elements at the tip (ring 0). This is followed by a large increase in accuracy, ≈ 20% for all refinements. After this initial jump, an overall average ≈ 1% increase in accuracy between 1 and 4 rings of elements occurs; however, it is noted that increasing the ring size in this region does not guarantee an improvement, and indeed, even with refinement, an improvement for an element ring size is not guaranteed. For instance, Thoutireddy and Ortiz 60 exploited the fact that the CF is dependent on the mesh for a mesh improvement strategy. If an increase in the element rings includes a section of poor mesh, the domain evaluated CF could be perturbed leading to a worse result.
The lack of CF angle convergence, both when increasing element ring size and refining the mesh for a ring size, is demonstrated in Figure 7A. However, when considering a fixed area, which is dictated by the element ring size on the most coarse mesh, Figure 7B, monotonic convergence is achieved with uniform homogeneous mesh refinement. The domain integral in its continuous form 14 is evaluated by considering an area around the crack tip. In the finite element  formulation if the area for computing the domain evaluated CF is fixed, and uniform homogeneous refinement occurs, the stress solution will improve in this area. Therefore, the domain integral, represented here as the summation of CF nodal values, will improve monotonically.
The next investigation demonstrates how varying element polynomial order around the crack tip affected the accuracy of the CF magnitude with homogeneous mesh refinement. The initial mesh is displayed in Figure 4B.
In Figure 8A, r d = 0.08 m and r p = 0.05 m, both are kept constant. The polynomial order p K of elements within r p at the crack tip was varied.
To demonstrate convergence rates for this problem, a value g = 256.9 N for the CF was obtained using a structured mesh with greater than 10 6 dof with h cF = 9.8 × 10 −5 m. This value is within the range presented in Section 6.1 by the empirical Equations 32 and 33. Figure 8A demonstrates that refinement in either h or p converges to a CF value of g = 256.9 N. For all polynomial orders, monotonic convergence was achieved with uniform h-refinement ( Figure 5). h-Refinement converges more efficiently than p-refinement as the stress is singular at the crack tip. 5 This agrees with the analytical convergence studies obtained in Pin and Pian 73 and Schwab 74 where it was shown for problems with a singularity h-refinement is more efficient than p-refinement. Overall a minimal error of 1.1% for p k = 7 against the upper bound of g.
Last, in Figure 8B, the rings of elements around the crack tip where p K = 5 is varied whilst r d remains constant at 0.08 m. It demonstrates that the accuracy of the CF solution is more dependent on the rings of elements around the crack, which have p K = 5 rather than r p . By having 2 element rings of high-order elements around the tip. The CF error reduces by 7.32% whilst the radius for considering 2 element rings decreases by 8.75 times.
To be consistent, the error values in Table 1 are found using the same empirical solution for this problem found in Miehe et al and Miehe and Gürses 21,22 ; here, g = 259.1 N. The SIPG obtains results in the range found in literature for CF values.

Single edge notched quasi-static crack propagation test
The SEN quasi-static crack propagation test, presented here, has the same geometry as the previous section. Here, weakly imposed heterogeneous Dirichlet boundary conditions are applied for a displacement of 0.01 m at the top and bottom of the specimen. These act in place of the uniaxial tensile stress shown in Figure 4B. The specimen had a Griffith failure criterion g c = 1000 N/m. The mesh is shown in Figure 9C, with an element length graded from 0.04 m around the expected crack path to 0.35 m. For this test, r p = 0.08 m and r d = 0.1 m, the domain of elements with a higher polynomial order moved with the crack tip.
The results for the instantaneous CF deviation from the expected crack direction, of 0 • , are shown in Figure 9A and 9B. Each figure shows a total of 8 element face separations. This equates to a total crack path length across the plate of 0.31 and 0.30 m, for Figure 9A and 9B, respectively. The figures demonstrate the improvements gained by using the domain approach. The domain method obtained a maximum difference, from 0 • , of 0.09 • compared to the tip evaluation, which achieved 0.55 • . The figures show how the path direction is governed by integration scheme more than the polynomial order of elements around the tip. The average difference between p K = 1 and p K = 7 for the domain approach was 0.029 • , and for the tip approach was 1.86 • .

Double notched 2-holed quasi-static crack propagation test
This benchmark, taken from Bouchard et al, 34 is a tension test of a double edge notched specimen with 2 holes for mixed mode problems; the geometry and loading conditions are shown in Figure 10A.
This test is necessary to show that SIPG can produce accurate results, compared to CG, for the mixed mode crack propagation problems. The plate acts in plane strain and has a shear modulus = 8 GPa, Poisson's ratio = 0.3, and a Griffith failure criterion of g c = 1000 N/m. In Figure 10B, the mesh is refined around the crack tips as in Miehe and Gürses 22 ; this is to ensure a more valid comparison.  Figure 11B.
A slight variant of the experiment, shown in Figure 11A, was performed with the same material properties, loading conditions, radii r d and r p , and polynomial order distribution. The experiment analysed the precision of the crack propagation. As the problem is symmetric, the 2 cracks paths in the problem should have the same normalised position relative to their starting point. The difference in the crack paths is a measure of the precision. The results are displayed in Figure 12A; however, it is important to note that the difference in the crack paths is empathise by the x and y axes being different scales.
The precision was measured as the maximum percentage difference from the mean of the 2 crack paths. The coarser mesh, element length of 0.25 m, obtains a precision of 20%. The refined mesh, element length of 0.123 m, achieved a precision of 2%. The lack of precision is caused by 2 features. First, for a coarse unstructured mesh, the stress field is poorly represented. This means on the first load step the configuration force is unlikely to be the same at both crack tips, and so the 2 cracks will propagate in slightly different directions. Secondly, as the increase in crack length Δo is larger for the coarse mesh, the error in crack path is magnified. This results in a diverging crack path, Figure 12A, and different stress fields at the tips, Figure 13B. The locations of the new crack tips and the stress fields, Figure 13a, now contrast more than if a finer mesh was used, Figure 12A. Ultimately, the difference between the stress fields and the error in crack path compounds the inaccuracy as the crack propagates through the specimen.

Mixed mode single crack propagation in a beam
In this section, a mixed mode single crack propagation problem in a double cantilever beam is analysed; the problem was first described by Belytschko and Black. 75 The problem is described in Figure 14A; it considers a fixed beam with a centre The crack step is dictated by the side length of elements in the mesh; Figure 14B elements are graded from 1∕400 m at the crack tip to 1∕200 m at the far field. The resultant crack path using (17) with the SIPG face splitting method is compared against the XFEM methods using the MCSC, 75,76 and against the cracking particle method, which also uses MCSC. 77 The resultant crack paths are displayed in Figure 15. It is important to note that the axes are different scales and the axes are also normalised with respect to the maximum length of the plate in the axes' orientation. Very similar results are obtained between CF SIPG, XFEM MCSC, 76 and the cracking particle method MCSC. 77 The only outlier is Belytschko and Black 75 ; however, the general shape of their crack path is similar. In this paper, we choose the CF method, 21 as this method maximises the power dissipation through crack growth, derived from the second law of thermodynamics.

FIGURE 15
Comparison of crack paths, for the crack in a beam problem 14A, using the SIPG CF crack path criterion (17), XFEM MCSC, 75,76 and the MCSC cracking particle method. 77 The axes x and y correspond to the same orientation as in Figure 14A (A) (B)

FIGURE 16
Cross crack: A, problem geometry where 2w = 1 m, the cracks at the centre have a length a = 0.1 m with the cracks labelled to 1 to 4. B, the deformed structured symmetric mesh with element face length h cF = 1∕30 m before propagation, with the displacement scaled by a factor 10 3

Quasi-static cross crack propagation in a finite plate
The final numerical example is taken from Sukumar and Belytschko, 78 shown in Figure 16A. The example is used here to investigate the consistency in CF values within a numerical experiment for 2 different boundary conditions, the average boundary conditions presented in Section 5.3 and a fixed bottom left corner node. The plate has a shear modulus = 8 GPa, Poisson's ratio = 0.3, and a Griffith failure criterion of g c = 1000 N/m. The traction acting perpendicular to the plate surface has a value of 10 MPa. The problem is a symmetric plate biaxial loaded, acting in plane strain, with 4 crack tips orientated in a cross. The cracks are propagated using the r-adaptive method combined with the face splitting algorithm presented in this paper. The increase in crack length for each crack propagation step is the same as the characteristic element size of the mesh, 1∕30 or 1∕50 m and corresponds to splitting of one element face. p K = 1 everywhere, and the tip method was used to investigate the effect the different boundary conditions had on the CF value at the crack tip node. This ensures that only the boundary conditions can cause any disparity in CF values at the crack tips. It is possible for the domain method to capture a different number of elements around the crack tip depending on the element distribution. Each crack was propagated at the same rate. The mesh generated was structured and symmetric; the mesh before crack propagation is shown in Figure 16B.
Each crack should undergo mode 1 propagation; the CF vector should therefore only have one non-zero component, which, in this example, is parallel with the crack lips. The first investigation is designed to isolate the effect the different boundary conditions have on the CF value at each crack tip. Therefore, a symmetric structured mesh is used as in Figure 16B. The magnitude of the deviation of the CF away from its expected direction for each crack tip is represented the by angle between its vector components as shown in Figure 17. Figure 17 shows that refining from h cF = 1∕30 m to h cF = 1∕50 m shows no improvement or even a worse representation of the crack path. The fixed corner makes the problem nonsymmetric by introducing a reaction force. Therefore, the boundary conditions are not sufficient for pure mode 1 crack propagation; hence, it is unachievable, even with refinement. However, the boundary condition does conserve symmetry in one respect as the deviation angle is the same for all cracks, within computational accuracy.
The range of the deviation angle for the 4 cracks using the average boundary conditions is shown Figure 17 by the greyed regions in the graph. The crack path deviation is within the order of computational error for both h cF = 1∕30 m and h cF = 1∕50 m. The example demonstrates that the average boundary conditions are suitable to obtain pure mode 1 fracture. The crack path deviation is at least an order 8 times more accurate than that obtained with the fixed corner boundary condition.
Finally, the norm of the CF values at each of the 4 cracks, see Figure 18B, is compared using an unstructured mesh with average boundary conditions and an average h F = 1∕30 m, Figure 18A. The cracks undergo 8 element face splits corresponding to a total crack length increase of ≈ 0.27 m. The percentage difference of the CF, from the average value, for each crack ranges from ±0.39% to ±1.24%. For the symmetric mesh in Figure 16B, the CF values were found to be the same for each crack. The range here is therefore generated by the nonsymmetry of the unstructured mesh. At a = 0.1 m, the mean CF value is 105.6 N; this is comparable to the XFEM solution calculated in Sukumar and Belytschko 78 where a value of 102.9 N was obtained, via a J-integral domain evaluation.

CONCLUSION
In this paper, we have used the thermodynamically consistent framework presented by Miehe and coworkers [21][22][23] to model brittle fracture for small strain problems in a SIPG finite element method. The proposed method exploits the element specific dof and the weak interaction between elements and existing as stiffness terms in the global stiffness matrix to propagate a crack in a fashion that is independent of the original element interface orientation. These benefits have not been taken advantage of in the past. In the proposed numerical method, a crack is propagated through (1) moving an element face in line with the CF, (2) removing the DG face stiffness values, in the global stiffness matrix corresponding to the reorientated face, (3) recalculating the local stiffness matrices of elements with changed geometry or polynomial order, and (4) updating the values global stiffness matrix. This r-adaptive procedure in DG adds no dof to the data structure, only the change of local element stiffness values in the global stiffness matrix is required. Additionally we have presented an algorithm to implement arbitrary high-order elements around a moving crack tip and produced an rp-adaptive scheme for SIPG methods. Good agreement for CF values and crack paths was obtained for static and quasi-static, single, and mixed mode problems using the proposed formulation. Mixed mode crack paths using the CF method were also compared against the MCSC method within the numerical framework of XFEM and the cracking particle method. The comparisons were made from results found in literature using existing finite element methods for the same problems. Finally, for nondeterminant systems, we introduce an average boundary condition that restrains rigid body motion leading to a solvable system. The proposed average boundary condition method has been shown to achieve machine accuracy for mode 1 crack propagation irrespective of element size for a cross crack problem.