Energy consistent framework for continuously evolving 3D crack propagation

This paper presents a formulation for brittle fracture in 3D elastic solids within the context of configurational mechanics. The local form of the first law of thermodynamics provides a condition for equilibrium of the crack front. The direction of the crack propagation is shown to be given by the direction of the configurational forces on the crack front that maximise the local dissipation. The evolving crack front is continuously resolved by the finite element mesh, without the need for face splitting or the use of enrichment techniques. A monolithic solution strategy is adopted, solving simultaneously for both the material displacements (i.e. crack extension) and the spatial displacements, is adopted. In order to trace the dissipative loading path, an arc-length procedure is developed that controls the incremental crack area growth. In order to maintain mesh quality, smoothing of the mesh is undertaken as a continuous process, together with face flipping, node merging and edge splitting where necessary. Hierarchical basis functions of arbitrary polynomial order are adopted to increase the order of approximation without the need to change the finite element mesh. Performance of the formulation is demonstrated by means of three representative numerical simulations, demonstrating both accuracy and robustness.


Introduction
The serious consequences of cracks in materials and structures mean that the computational modelling of crack propagation continues to be a critical area of research, and a major challenge. This paper presents a finite element based computational framework for modelling brittle crack propagation in elastic solids, based on the concept of configurational mechanics. The focus of this paper is on both the mathematical formulation and the computational framework to model and continuously resolve propagating cracks in a robust and computationally tractable manner. This paper is concerned with quasi-static problems where the influence of inertia is ignored. A sequel to this paper will demonstrate the extension of this work to dynamic fracture.
Within the context of the finite element method, the modelling of crack propagation was traditionally undertaken in a smeared sense, whereby the problem is solved within a continuum setting, without the need for approximation of discontinuities or changing mesh connectivity. However, as strains localise, numerical difficulties arise, and regularisation is required. In contrast, discrete approaches are able to directly approximate the macroscopic crack geometry and describe cracks in a more natural and straightforward manner in terms of displacement jumps and tractions. One of the first methods to simulate discrete crack propagation, without predefining the crack surfaces, was the embedded discontinuity method [1,2]. Since the enhancement in this method was local, it was computationally efficient, but the piecewise nature of the resulting crack surface led to numerical difficulties. Partition of Unity methods, such as the extended finite element method, represented a significant advancement to overcome this, see [3,4].
Tracking the crack path, irrespective of the methodology for resolving the resulting crack surface, is also a significant challenge. Jäger et al. [5] reviewed different discrete crack tracking algorithms for quasi-brittle fracture and undertook a systematic comparison, motivated by the challenge of undertaking crack propagation in threedimensions. They demonstrated the difficulty in predicting and representing crack surfaces of arbitrary shape and which are continuous across element boundaries, and identified a global tracking algorithm by Oliver [6] as the most algorithmically robust, yet simple and easily implemented into existing finite element software. The criterion for crack propagation was the classical principal stress-based Rankine criterion. However, in this paper, the need for a crack tracking methodology is obsolete, since the crack front is simply advanced by an amount and in a direction based only on physical laws.
The variational formulation for fracture by Francfort and Marigo [7] represented a watershed in the formulation and numerical analysis of crack evolution, recasting and generalising the engineering approach to fracture proposed by Griffith into a continuum mechanics framework as an energy minimisation problem. This removed some of the constraints of Griffith's original approach, for example any assumption as to the actual geometry of the crack path or its restriction to stable crack propagation. Numerous developments have emerged from this seminal work, many focusing on the regularisation of the complex energy functional to make it more suitable for a numerical modelling. Most notably, the popular phase-field methods for fracture [8] have emerged where cracks are represented by a smooth damage variable leading to a phase-field approximation of the variational formulation for brittle fracture [9,10]. Although the crack surface is not represented by a discrete discontinuity, phase-field methods can be implemented in a straightforward manner using coupled multi-field finite element solvers.
This paper is focused on the concept of configurational mechanics that dates back to the original work of Eshelby and his study of forces acting on material defects [11,12]. The concept of configurational (or material) forces is now a well established method to evaluate defects in a material providing a unified framework for the analysis of material imperfections and has been adopted by, amongst others, Maugin [13,14]. Steinmann [15,16] developed a computational strategy for the assessment of fractured bodies. Gurtin and Podio-Guidugli [17] derived the equations that govern the motion of the crack front, applied to dynamic problems in two dimensions, establishing a link between Griffith's fracture criterion and configurational forces in two dimensions. Miehe et al. [18,19] and Kaczmarczyk et al. [20] build on this work to establish a finite element methodology for crack propagation.
To formulate the crack propagation problem within the framework of configurational mechanics, two related kinematic descriptions are defined in the spatial and material settings. In the former, the classical conservation law of linear momentum balance is described, where Newtonian forces are work conjugate to changes in the spatial position, at fixed material position (i.e. no crack propagation). In the material setting, which represents a dual to the spatial setting, an equivalent conservation law is described, where configurational forces are conjugate to changes in material position but with no spatial motion. This decomposition of the behaviour is proven to be a simple but powerful methodology for describing crack propagation. The authors' previous paper [20] described the mathematical formulation for crack propagation and a methodology for resolving the evolving crack path within the context of the finite element method, and represented an advancement of the work of Miehe et al. [18,19]. The current paper utilises the same underpinning mathematical formulation but introducing new computational strategies to continuously resolve the evolving crack path.
This paper brings together several key developments, establishing a new methodology for simulating 3D brittle crack propagation, representing a significant improvement on the authors' previous work [20]. Griffith's fracture criterion is expressed correctly in terms of configurational forces. An expression for equilibrium of the crack front is established, balancing the configurational forces on the crack front with the resistance of the material. This is exploited so that the crack front can advance continuously, moving nodes on the crack front to the physically correct position without recourse to any kind of crack tracking algorithm. There is no need for mesh splitting, or changes to mesh connectivity in order to resolve the crack front advancement, as with other approaches, e.g. [18][19][20]. To maintain mesh quality, a mesh smoothing strategy, with surface constraints, is presented as a continuous process as part of a problem-tailored Arbitrary Lagrangian Eulerian formulation. This is supplemented by face flipping, node merging and edge splitting where necessary. Post-processing of the stresses is not required to determine if the crack should propagate and the crack front shape is calculated based purely on the physical equations. The spatial and material displacement fields are both discretised using the same finite element mesh, although we adopt different levels of approximation for the two fields. The resulting discretised weak form of the two conservation equations represents a set of coupled, nonlinear, algebraic equations that is solved in a monolithic manner using a Newton-Raphson scheme. In addition, an arc-length method is adopted to trace the dissipative load path for brittle fracture propagation, using crack area rather than displacements as a control. Three numerical examples are presented that demonstrate the ability of the formulation to accurately and robustly predict crack paths without bias from the original mesh. Fig. 1 shows an elastic body with an initial crack in the reference material domain. As a result of loading, the crack extends and the body deforms elastically. It is convenient to decompose this behaviour into a purely configurational change, i.e. crack extension, which is described by the mapping from the reference material domain to the current material domain (Ξ ), followed by elastic deformation only, described by the mapping from the current material to spatial domain (ϕ). We utilise these mappings to independently observe the evolution of the crack surface in the material domain B t and the elastic deformation of solid in the spatial domain Ω t .

Body and crack kinematics
The material coordinates X are mapped onto the spatial coordinates x via the familiar deformation map ϕ(X, t). The physical displacement is: The reference material domain describes the body before crack extension. Ξ (χ , t) maps the reference material coordinates χ on to the current material coordinates X, representing a configurational change, i.e. extension of the crack due to advancement of the crack front. Φ maps the reference material coordinates χ on to the spatial coordinates x. The current material and spatial displacement fields are given as H and h are the gradients of the material and spatial maps and F the deformation gradient [20], defined as:  Given that the physical material cannot penetrate itself or reverse the movement of material coordinates, we have: In addition, the time derivative of the physical displacement u and the deformation gradient F (material time derivative) are given as [20]: The crack surface, comprising two crack faces, is denoted as Γ and the crack front is denoted as ∂Γ , see Fig. 2. In [20], a kinematic relationship between the change in the crack surface areaȦ Γ and the crack front velocityẆ was derived that is given as: where A ∂Γ is a dimensionless kinematic state variable that defines the current orientation of the crack front. In deriving this expression, it was recognised that any change in the crack surface areaȦ Γ in the current material space can only occur due to motion of the crack front.

Dissipation of energy due to creation of new crack surfaces
The first law of thermodynamics can be expressed as: where the left hand side is the power of external work, the first term on the right hand side is the rate of the crack surface internal energy and the last term is the rate of internal body energy. γ is the surface energy [N m −1 ] and Ψ is the volume specific free energy. Making use of Eqs. (5) and (6), and given that dV = ∇ X ·ẆdV , this expression can be rewritten as ∫ where are the first Piola-Kirchhoff and Eshelby stress tensors, respectively. The former is the familiar driving force for elastic deformation in the spatial domain, whereas the Eshelby stress is its material counterpart and the driving force for local configurational changes.
In order to get a local form of the first law, the Gauss divergence theorem is applied to the last integral in Eq. (8) resulting in the following expression To simplify this equation, we recognise that, in the limit, the surface C ( Fig. 2) collapses to the crack front ∂Γ , and integrals over the crack front can be expressed as In addition, the first and second terms on the right hand side of Eq. (10) vanish since the spatial and material conservation laws of linear momentum balance are expressed as follows: Moreover, considering only admissible velocity fields and stress fields in equilibrium with external forces, Eq. (10) can be expressed as: The configurational force is defined as Thus, the local form of the first law (Eq. (10)) is expressed as: This equation represents the equilibrium condition for the crack front. The term γ A ∂Γ can be considered the material resistance. The exploitation of this expression, to describe a formulation and computational framework for crack propagation, is the essence of this paper's novelty.
In the case of an evolving crack front, we can deduce thatẆ ̸ = 0 and γ A ∂Γ = G. It should be noted that for a continuous elastic body comprising a homogeneous material, the configurational forces G within the volume, away from the crack front, should be zero.
Since all energy dissipation is restricted to creation of new crack surfaces, it follows that the local form of the second law is given as where D is the dissipation of energy per unit length of the crack front. This inequality restricts evolution of the crack to positive crack area growth at each point of the crack front. Although the first law defines if the crack front is in equilibrium and the second law places restrictions on the direction of crack evolution, it does not determine how A ∂Γ orẆ evolves. A crack growth criterion is presented in the next section.

Evolution of the crack front
A straightforward criterion for crack growth, equivalent to Griffith's criterion, is proposed: where g c = 2γ is a material parameter specifying the critical threshold of energy release per unit area of the crack surface Γ , also known as the Griffith energy. For a point on the crack front to satisfy the crack growth criterion, either To determine evolution of a point on the crack front, we adopt the principle of maximum dissipation, that states that, for all possible configurational forces G * satisfying the crack growth criterion φ(G * ) = 0, the dissipation D attains its maximum for the actual configurational force G. Therefore, we have This can be interpreted as an unconstrained minimisation problem, for which the Lagrangian function is: withκ the Lagrange multiplier. The Kuhn-Tucker conditions become Thus, for a point on the evolving crack front, the crack front orientation is colinear to the configurational force (Eq. (15)), i.e. γ A ∂Γ = G, and the crack extension is given asẆ =κA ∂Γ .κ has the dimension of length and represents the kinematic state variable for a point on the crack front, which can be identified as:

Spatial and material discretisation
Finite element approximation is applied for displacements in both the current material and physical space where superscript h indicates approximation and (·) indicates nodal values. Three-dimensional domains are discretised with tetrahedral finite elements. In the spatial domain, hierarchical basis functions of arbitrary polynomial order are applied, following the work of Ainsworth and Coyle [21]. This enables the use of elements with variable, nonuniform orders of approximation, with conformity enforced across element boundaries. In the material domain, linear approximation is adopted, as this is sufficient for describing the crack front. The discretised gradients of deformation are expressed as The residual force vector in the discretised spatial domain is expressed in the classical way as: where τ is the unknown scalar load factor and f h ext,s and f h int,s are the vectors of external and internal forces, respectively, and defined as follows where ATET and ATRI indicate the standard FE assembly for tetrahedral elements and triangular faces, respectively.
The discretised version of Eq. (15) establishes the material counterpart to Eq. (24), expressed as The nodal configurational forces,G h , are the driving force for crack front evolution and, for the purposes of establishing equilibrium of the crack front, are calculated only for nodes on the crack front as f h res is the vector of nodal material resistance forces, given as: where g c is a vector of size equal to the number of nodes in the mesh, with zero for all components except those associated with nodes on the crack front, where the value is g c . The matrixÃ h Γ , defined in [20], has dimensions of length and describes the current orientation of the crack surface. Each row of matrix f h res can be interpreted as containing the components of a node's material velocityẆ which lead to a change of crack area. Since the crack surface area can only change due to material displacements of the nodes of the crack front,Ã h Γ has relatively few non-zero rows, equal to the number of nodes on the crack front. The nodal forces of material resistance f h res have dimensions of force and are work conjugate to the material displacement W on the crack front.

Arc-length control
The global equilibrium solution for the spatial and material displacements is obtained as a fully coupled problem using the Newton-Raphson method, converging when the norms of the residuals are less than a given tolerance. To trace the nonlinear response resulting from the dissipative behaviour, an arc-length technique is adopted. The system of equations for conservation of the material and spatial momentum is supplemented by a load control equation that imposes a constraint that controls the crack area increment during each load step. This load control equation takes the form: n+1 ) I − ∆A Γ = 0, for I ∈ {I : N I is crack front node} (29) where ∆A Γ is the prescribed increment of the crack area, i + 1 is the current iteration of load step n + 1.

Resolution of the propagating crack and mesh quality control
In the previous work of the authors [20], the crack front equilibrium equation (15) was established but not fully exploited. This is addressed in this paper, establishing an implicit, continuous crack front evolution formulation. Previously, a discrete face-splitting methodology was adopted, whereby the mesh was changed to create the new crack surfaces and crack front. The new crack front was generated by identifying element faces ahead of the current crack front, aligning them to the direction of the configurational forces and finally splitting these faces if the crack criterion was violated, creating a displacement discontinuity. This process was continued until equilibrium was achieved. The approach in the current paper is quite different, whereby the crack front advances in an implicit, continuous manner, to establish equilibrium of the crack. Mesh nodes are subsequently moved (the mesh is not split or the connectivity changed) in order to resolve the new crack geometry.
However, in the process of moving the mesh nodes to resolve the advancing crack front, the mesh can become distorted, potentially creating poor quality elements, leading to numerical errors. To mitigate this effect we adopt several strategies: 1. Edge splitting is applied to elements behind the crack front that have become too elongated. This action will result in the creation of new elements. It is enforced if, for a given node, there exists an adjacent edge with a length greater than 1.5 times the average edge length of all adjacent edges to the node. 2. Node merging is applied to elements ahead of the crack front that have become too contracted. This action will result in the removal of elements. It is enforced if, for a given node, there exists an adjacent edge with a length less than 1/3 of the length of the longest edge adjacent to the node.
3. Face flipping is applied to elements in the vicinity of the crack front to ensure that a 3D Delaunay triangulation exists, with optimal internal angles. This is described in more detail in Section 6.1.
These procedures are utilised, if necessary, at the beginning of each load step, before the Newton-Raphson iterations begin, when the solution is already out of equilibrium. Furthermore, in the case when new nodes are added, variables are transferred to the new mesh based upon approximation of the variables using the old mesh. Fig. 3 demonstrates how the crack front evolves using the example of a three-point bending of a beam with an initial corner notch. Only the projection of the FE mesh onto the crack surface is shown. Crack surface A shows an equilibrium solution, with the configurational forces identified. Crack surface B shows a subsequent configuration, where the crack front has advanced to a new equilibrium position, the nodes have moved to accommodate the crack geometry but there has been no changes to the mesh connectivity. Crack surface C represents a further equilibrium configuration. The crack front has continued to advance but, to maintain mesh quality, the mesh behind the crack front has evolved, with new elements being created using the edge splitting procedure as required.

Face flipping
At the beginning of each load step, when the solution is out of equilibrium, a patch of elements around the crack front is checked to ensure that it represents a 3D Delaunay triangulation. Fig. 4 demonstrates the idea in 2D. Considering the two elements on the left, edge i − k is prohibited because it lies in the interior of the circle that intersects the nodes of element i − j − k (and element i − k − l). Flipping edge i − k will address this problem, redefining the two adjacent elements, without affecting the rest of the mesh. Thus, edge i − k is removed and replaced by edge j − l, and two new adjacent elements are formed that represent a Delaunay triangulation.

Mesh quality control
In addition to the mesh adjustments described above, a global mesh quality control procedure is also adopted in the spirit of an Arbitrary Lagrangian Eulerian (ALE) method. During the global Newton-Raphson procedure, constraints are imposed on the shape of elements to ensure good mesh quality, but without influencing the physical response. Consequently, mesh smoothing is a continuous process in the current material domain, i.e. the domain in which the topology evolution occurs.
The authors have proposed a volume-length measure of element quality [20,22]. This measure does not directly determine dihedral angles; however, it has been shown to be very effective at eliminating poor angles, thus improving stiffness matrix conditioning and reducing interpolation errors [23,24]. As the volume-length measure is a smooth function of node positions, and its gradient is straightforward and computationally inexpensive to calculate, it is ideal for the problem at hand. The volume-length quality measure is defined as where q 0 , V 0 and l rms,0 are the element quality, element volume and root mean square of the element's edge lengths respectively, in the reference configuration. H h is the material deformation gradient, b is a measure of element quality change, relative to the reference configuration, and dl rms = l rms /l rms,0 is the stretch of l rms,0 . Element quality q 0 is normalised so that an equilateral element has quality q 0 = 1 and a degenerate element (zero volume) has q 0 = 0. Furthermore, b = 1 corresponds to no change and b = 0 is a change leading to a degenerate element. An element edge length in the current material configuration is expressed as where ∆χ j is the distance vector of edge j in the reference configuration. Thus l rms is calculated as To control the quality of elements, an admissible deformation H h is enforced such that b(H h ) > γ for γ ∈ (0, 1).
In practice, 0.1 < γ < 0.5 gives good results. This constraint on b is enforced by applying a volume-length log-barrier function [25] defined, for the entire mesh, as where B is the barrier function for the change in element quality in the current material configuration and N el is the number of elements. It can be seen that the log-barrier function rapidly increases as the quality of an element reduces, and tends to infinity when the quality approaches the barrier γ , thus achieving the objective of penalising the worst quality elements. In order to build a solution scheme that incorporates a stabilising force that controls element quality, a pseudo 'stress' at the element level is defined as a counterpart to the first Piola-Kirchhoff stress as follows where matrixQ is defined as followŝ It is worth noting that Q should be a zero matrix for a purely volumetric change or rigid body movement of a tetrahedral element.
It is now possible to compute a vector of nodal pseudo 'forces' associated with Q as

Shape preserving constraints
The continuous adaption of the mesh must be constrained to preserve the surface of the domain being analysed, including the crack surface. Nodes on the surface can only be allowed to slide along the surface and not deviate from it. In this work we start from the observation that a body's shape can be globally characterised as a constant: where V is the volume of the body and A is its surface area. This can be expressed in integral form as ∫ where X s are the coordinates on the surface S in the current material domain. Applying Gauss theorem we obtain ∫ ∂B 0 where N is the outer normal to the surface. The local form of this equation is given as N where 3C is a constant given for the reference geometry. If the above constraint is satisfied, the body shape and volume are conserved locally. An equivalent form this constraint equation is N where χ s are the original positions of the surface.
To enforce the mesh quality control described in the preceding section, subject to the surface constraints described in this section, the method of Lagrange multipliers is used, with the following functional, where H h and Q are functions of material positions,X and are current mesh nodal positions defined in the previous sections. ξ and η are parametrisation of surface ∂B 0 . Φ λ and Φ X are the shape functions for Lagrange multipliers and material coordinates respectively. Φ λ are piecewise continuous functions with order of approximation equal to that of the material coordinates.
Calculating the stationary values of the Lagrangian results in the first ( L(X,λ) ,X = 0 ) and second ( L(X,λ) ,λ = 0 ) Euler equations. The first is given as Taking a truncated Taylor series expansion after the linear term of this nonlinear equation results in where the differential operators are defined as Linearising the second Euler equation leads to ∫ These surface constraint equations are applied for each surface patch independently, i.e. where two surfaces meet, two independent sets of equations with Lagrange multipliers are applied. Moreover, these geometry preserving constraints are not applied to the crack front, since configurational forces drive the material displacement of those nodes.

Linearised system of equations and implementation
A standard linearisation procedure is applied to the residuals r h m , r h s , f h q & r τ . For the moment, the surface constraints defined in Section 6.3 are excluded. The global equilibrium solution is obtained using the Newton-Raphson method, solving for the iterative changes in spatial displacements, current material displacements and the load factor as a fully coupled problem. Since the material residual is non-zero only for nodes on the crack front, it is convenient, for the purposes of presentation, to decompose the material nodal positions into those at the crack frontX f , those associated with surfacesX s (crack surface and body surfaces) and the rest of the meshX b , i.e.X =X f ∪X b ∪X s . The resulting linear system of equations for iteration i and load step n + 1 is given as: . The vectors f h q b and f h q s are the components of f h q associated withX b andX s respectively. This can be simplified for presentation purposes as: whereq a is the vector of all unknowns excluding the nodal coordinates of the surfaces,X s , f h a the corresponding terms on the right hand side of Eq. (49).
In order to preserve the surfaces of the body during the analysis, it is necessary to impose the constraints described in Section 6.3 on the coordinates of the surface nodes,X s . Thus, the above system of equations is augmented as follows: and When solving this nonlinear system of equations, convergence is quadratic and typically requires 3-4 iterations per load step to achieve convergence. We adopt a Total ALE approach; at the beginning of each new load step the current material mesh becomes the new reference mesh. Discretisation is undertaken using 3D tetrahedral elements with hierarchical basis functions of arbitrary polynomial order [21] in the spatial domain. A linear approximation is adopted in the material domain.
The solution strategy presented in this paper is implemented for parallel shared memory computers, utilising opensource libraries. MOAB, a mesh-oriented database [26], is used to store mesh data, including input and output operations and information about mesh topology. PETSc (Portable, Extensible Toolkit for Scientific Computation [27]) is used for parallel matrix and vector operations, the solution of linear system of equations and other algebraic operations. MOAB and PETSc are integrated in MoFEM [28], which is an open source code for finite element analysis, developed and maintained by our group.

Numerical examples
Three numerical examples of crack propagation in brittle materials are presented. Two examples consider the fracture of unirradiated nuclear graphite. The third example considers the fracture of PMMA.

Graphite cylinder slice test
This numerical example considers a slice of a graphite cylindrical brick, placed in a loading rig and loaded as shown in Fig. 5. The red box indicates the part modelled in the numerical analysis. The loose key adaptors on the left are fully fixed along their left hand side. The numerical load is applied to the mid-point of the crosshead beam. The brick slice is 25 mm thick, the Young's modulus is 10,900 MPa, Poisson's ratio is 0.2 and Griffith energy is 0.23 N/mm. The specimen is loaded via the key adaptors. Contact between the key adaptors and the graphite slice only occurs on one edge of the loose keyway and this is modelled by tied degrees of freedom on those edges.
The mesh for this example is shown in Fig. 6(a) and consists of 6033 tetrahedral elements. The numerical analyses were undertaken using one mesh but repeated for 1st, 2nd and 3rd-order approximations. The numerically predicted crack path is shown in Fig. 6(b). It can be seen that the crack propagates from the keyway corner with a curved trajectory to the free surface of the inner bore. These compare well with the experimentally determined crack paths shown in Fig. 7. The elastic energy versus the crack area and the force versus displacement plots are shown in Figs. 8(a) and (b) respectively. The displacement in Fig. 8(b) is known as the generalised displacement and does not represent a particular point on the structure, but its value is work conjugate to the applied forces and is calculated as u = 2Ψ /τ f , where f = 1N is the reference force, Ψ is the elastic energy, τ is the load factor and u is the generalised displacement. The arc-length control, described earlier, is used to trace the nonlinear response.  Snapshots of the brick slice at three points S1, S2 and S3, shown in Fig. 8 are shown in Fig. 9(a), (b) and (c) respectively. It should be noted that Fig. 8 demonstrates that, as compared to the fluctuating results for 1st-order analysis, the results for 2nd and 3rd-order analyses are very smooth. The fluctuation in the results for the 1st order analysis is due to shear locking.

Mixed mode I-III loading in three-point bending test
The second example assesses the ability of the proposed approach to simulate the mixed mode I-III experiments of [29] on PMMA pre-notched beams loaded in a three-point-bending configuration. This problem was also simulated by Pandolfi and Ortiz [30] (although they did not model the entire specimen, instead restricting calculations to the central section of the specimen, applying equivalent static boundary conditions representative of three-point bending). The specimens are 260 mm long, 60 mm deep, and 10 mm thick, pre-cracked at 45 • , 60 • , and 75 • to the front face,    9. Snapshots for points S1, S2 and S3 shown in Fig. 8 for graphite cylinder slice test. The experiments demonstrated that for homogeneous and linear-elastic isotropic materials the growing crack, loaded in mixed I-III modes, reorientates towards a pure mode I situation. The meshes are shown in Fig. 11 and consist of 8696, 7824 and 8705 tetrahedral elements for 45 • , 60 • , and 75 • initial notch angles respectively. All three cases are solved for 1st, 2nd and 3rd-order of approximation. Elastic energy versus crack area plots are shown in Fig. 12(a). Consistent with the previous example, the response for all three cases fluctuates for the 1st-order analysis but is smooth for the 2nd and 3rd-order analyses. Fig. 12(b) compares the load-displacement response for γ = 60 • , for different orders of approximation. Fig. 12(c) compares the load-displacement response for 2nd-order approximations, for all three initial notch angles. The ultimate load for γ = 45 • , γ = 60 • and γ = 75 • is 0.46 kN, Fig. 10. Mixed mode I-III crack growth in PMMA specimen loaded in a three-point bending configuration [29]. The specimen length is L = 260 mm, L e = 120 mm, depth w = 60 mm, thickness t = 10 mm, notch height a = 20 mm and notch inclination γ = 45 • , 60 • and 75 • .  0.4 kN and 0.38 kN respectively, decreasing with increasing initial notch angle. For the γ = 45 • case, snapshots of the evolving crack for the 2nd-order analyses associated with points S1, S2 and S3 (see Fig. 12(c)) are shown in Fig. 12(d). Furthermore, the final crack surfaces for γ = 75 • , γ = 60 • and γ = 45 • are shown in Fig. 13(a), (b) and (c) respectively, which clearly shows the re-orientation towards a mode I situation. Table 1 compares the experimentally measured average kink angles, obtained by averaging the kink angles on the front and back surfaces of the specimen, and the corresponding numerical values. The numerical results capture the general trend of increasing kink angle with increasing pre-crack inclination and mode III component. The difference between the numerical and experimental results is very small.

Torsion test
An experimental study by Brokenshire [31] of a torsion test of a plain concrete notched prismatic beam (400mm × 100mm × 100mm) has been repeated for nuclear graphite. The experimental procedure and full details of      the boundary conditions and dimensions for the original study are described in Jefferson et al. [32] and illustrated in Fig. 14(a). The notch is placed at an oblique angle across the beam and extends to half the depth. The beam is placed in a steel loading frame, supported at three corners and loaded at the fourth corner. The material properties used for this example are the same as those used in the graphite cylinder slice test example. The beam and steel frame are discretised using tetrahedral elements and are shown in Fig. 14(b). The mesh consists of 11890 elements and the problem is solved with 1st, 2nd and 3rd-order approximation. Fig. 15(a) and (b) shows the elastic energy-crack area and load-displacement response for the three numerical tests respectively. Good numerical convergence is observed with increasing order of approximation and the arc-length control method is able to track the dissipative loading path. The ultimate load for 1st, 2nd and 3rd-order approximation is 0.326 kN, 0.294 kN and 0.291 kN respectively. The snapshots of evolving crack at points S1, S2 and S3 (see Fig. 15(b)), are shown in Fig. 16(a), (b) and (c) respectively. The front and top views of the simulated crack surface are also shown in Fig. 17(a) and (b) respectively, which clearly shows the complicated crack path. It is worth noting that the same problem was discussed in our previous paper [20] but the material was concrete rather than graphite. In this previous case the crack (a) Snapshot at S1.
(b) Snapshot at S2. (c) Snapshot at S3.  path showed excellent agreement with experiments but over predicted the experimental ultimate load by approximately 2.5 times. This was a consequence of assuming linear elastic fracture mechanics for a problem where the size of the fracture process zone is significant compared to the size of the problem. This is not an issue in the current situation because graphite's microstructure is significantly smaller than for concrete and the assumption of linear elastic fracture mechanics is valid.

Conclusions
A novel formulation for brittle fracture in elastic solids within the context of configurational mechanics has been presented that incorporate several advancements on the authors' previous work. The local form of the first law of thermodynamics provides a condition for equilibrium of the crack front and the direction of crack propagation is given by the configurational forces on the crack front, maximising the local dissipation. These are exploited such that the advancing crack front is continuously resolved by moving the nodes of the finite element mesh, without the need for face splitting or the use of enrichment techniques.
A monolithic solution strategy has been described that simultaneously solves for both the material displacements (i.e. crack extension) and the spatial displacements. In order to trace the dissipative loading path, an arc-length procedure has been developed that controls the incremental crack area growth. In order to maintain mesh quality, smoothing of the mesh is undertaken as a continuous process, together with face flipping, node merging and edge splitting where necessary. Hierarchical basis functions of arbitrary polynomial order are adopted to increase the order of approximation without the need to change the finite element mesh.
Three numerical examples have been presented to demonstrate both the accuracy and robustness of the formulation. Convergence studies have been undertaken in all cases. All three problems demonstrate the ability to simulate experimental crack paths.