An enriched phase-field method for the efficient simulation of fracture processes

The efficient simulation of complex fracture processes is still a challenging task. In this contribution, an enriched phase-field method for the simulation of 2D fracture processes is presented. It has the potential to drastically reduce computational cost compared to the classical phase-field method (PFM). The method is based on the combination of a phase-field approach with an ansatz transformation for the simulation of fracture processes and an enrichment technique for the displacement field as it is used in the extended finite element method (XFEM) or generalised finite element method (GFEM). This combination allows for the application of significantly coarser meshes than it is possible in PFM while still obtaining accurate solutions. In contrast to classical XFEM / GFEM, the presented method does not require level set techniques or explicit representations of crack geometries, considerably simplifying the simulation of crack initiation, propagation, and coalescence. The efficiency and accuracy of this new method is shown in 2D simulations.


Introduction
During the last decade, the phase-field method (PFM) has become one of the most popular techniques for the simulation of fracture processes. Starting from the variational formulation of Griffith's theory of fracture [25] and its regularised approach by Bourdin et al. [13], there now exists a large body of research, especially regarding phase-field formulations of brittle fracture. Comprehensive overviews of the method can be found in [1,18,64], additionally we refer to [5,14,37,46,47,61]. Proceeding from brittle fracture, many contributions branched out into the investigation of phasefield approaches to ductile fracture, e.g. [2,3,12,22,39,48], dynamics, e.g. [10,15,34,36,55,59] and more recently fatigue processes, e.g. [4,16,56,57]. The big advantage of the PFM compared to other popular techniques like the extended finite element method (XFEM) or the generalised finite element method (GFEM) is the ability to simulate crack nucleation, propagation, branching, and coalescence processes in 2D and B Stefan Loehnert stefan.loehnert@tu-dresden.de 1 Technische Universität Dresden, Institute of Mechanics and Shell Structures, August-Bebel-Str. 30, 01219 Dresden, Germany 3D without the necessity for crack tracking algorithms. However, one of the biggest disadvantages of the PFM is that a very fine discretisation of the domain is required and, already for linear elastic fracture problems, a highly non-linear equation system has to be solved, resulting in a high computational effort even for 2D simulations. Therefore, several techniques have been developed to improve the solution process and speed up simulations [1,29,32,35,46].
Unlike the PFM, the XFEM / GFEM [9,20,27,45,50] allows for the accurate simulation of fracture processes with rather coarse meshes. For linear elastic fracture mechanics, the resulting equation system remains linear, leading to a much lower computational effort compared to PFM simulations. However, within the XFEM / GFEM crack tracking techniques like the level set method [53,60], the fast marching method [17,58] or explicit crack advancement techniques [21] are required to reproduce the change of the crack geometry in space during propagation. Hence, phenomena like crack coalescence and branching involve significantly more difficulties, especially in 3D [7]. In addition to the necessity for crack tracking techniques, criteria determining whether the crack propagates, in which direction and how far it propagates need to be evaluated.
During the last years, several techniques were developed to overcome these difficulties by combining the XFEM with the PFM or a damage model, respectively. In [28] both methods are coupled in a multiscale context. The PFM is applied in a fine scale simulation in the vicinity of a crack tip to determine the local crack advancement. In a concurrent coarse scale simulation, the XFEM is applied to account for the crack behaviour away from the propagating crack tip(s). The so-called Xfield-method [31] is based upon two different grids: the development of the crack at the crack tip is simulated by the PFM on a moving fine grid, behind the crack front the crack path is represented by level sets in the context of the XFEM on a coarse background mesh. An alternative approach is the thick level set method [49,62]. Here, a damage model is applied within a small domain defined by a single level set function. Enrichment functions as they are used within the XFEM improve the reproduction of displacement discontinuities in parts of the domain that are fully damaged.
To take advantage of the knowledge that the solution of the phase-field differential equation has exponential quality, in [38] exponential shape functions for the phase-field in two dimensions are proposed. The presented ansatz can reduce the number of elements perpendicular to the crack enormously, however, parallel to the crack a rather fine mesh is still necessary. Another drawback of the approach is, that the crack path must be known a priori because the crack can only develop along element edges. This foils one of the main advantages of the PFM: the mesh independent development of cracks.
For heterogeneities, a new phase-field approach has been introduced by Finel et al. [23]. The basic idea of this socalled sharp phase-field method (SPFM) is to transform the ansatz space in such a way, that the solution can be captured more accurately for a given discretisation. The transformed ansatz functions directly depend on standard ansatz functions interpolating between nodes. This method was first used in [24] and has been successfully applied to grain growth in the context of the finite difference method [19].
In the present contribution, the respective advantageous properties of the PFM and the XFEM are combined to the extended phase-field method (XPFM) designed to simulate fracture processes. A transformed ansatz is used to reproduce the exponential character of the phase-field function. Additional enrichment functions are incorporated into the ansatz space of the displacement field to capture the discontinuity (i.e. the jump across the crack surface) behind the crack front and the transition zone from fully degraded material to sound material ahead of the crack front. These enrichment functions are directly determined from the degrees of freedom of the phase-field. No additional level set functions or explicit representations of cracks are necessary, and geometrical crack tracking algorithms are not required.
The article is structured as follows. In Sect. 2, the standard PFM for brittle fracture is summarised and discussed. Particular attention is given to some discretisation aspects.
Next (Sect. 3), the XPFM for 2D is derived from some preliminary investigations in one dimension. Furthermore, stabilisation and blending techniques are discussed. The algorithmic framework is introduced in Sect. 4. Numerical and analytical studies regarding quadrature and the different model parameters are conducted in Sect. 5. In Sect. 6 numerical examples are shown and compared to results from literature to demonstrate the applicability and efficiency of the XPFM. A conclusion and an outlook to future work is given in Sect. 7.

Phase-field model of brittle fracture
In this section, the basic equations of the classical phase-field approach to fracture are summarised. Some one-dimensional considerations and a discussion of discretisation aspects motivate the XPFM presented in the subsequent sections.

Governing equations
Following Griffith's theory of brittle fracture [33], a predefined crack in an elastic solid body B ⊂ R 2 will only propagate if the energy release rate is at least equal to the surface energy dissipated by the formation of new crack surfaces. Francfort and Marigo [25] reformulate this as minimisation problem of the total energy E tot . It is defined as The elastic strain energy stored in the body, with the elastic strain energy density function ψ 0 , dependent on the displacements u, and the external work due to external body forces b and tractionst are considered. ∂B t ⊂ ∂B describes the Neumann boundary for imposed tractions on the body and ∂B u ⊂ ∂B is the Dirichlet boundary where displacement boundary conditionsū are applied. Assuming small deformation theory and isotropic linear elastic material behaviour, the strain tensor ε is defined by The strain energy density function of an undamaged body is given by with λ being the first Lamé parameter and μ being the shear modulus. With σ 0 := ∂ψ 0 /∂ε the stresses of the intact material result in Here, 1 denotes the second order identity tensor. The formation of new surfaces due to crack growth evokes a certain energy dissipation, which is called the fracture surface energy E s . It is proportional to the newly created surface area Γ through the material-depended critical fracture energy release rate G c Due to its dissipative nature E s is only allowed to grow whereas the elastic part E e of the total energy can be recovered.
To treat the energy minimisation approach numerically a damage type parameter φ, the so-called phase-field parameter, is introduced in [13]. It interpolates smoothly between the fully broken state (φ = 1) and the intact state (φ = 0) of the material. Depending on this phase-transition, the strain energy density and the stresses can be degraded by a degradation function g(φ). This function must be monotonically decreasing and has to fulfil the properties g(0) = 1, g(1) = 0, g (0) < 0 and g (1) = 0. In this contribution, for simplicity the degradation function most frequently applied in literature is employed: resulting in and Herein, k g is a small stabilisation parameter which ensures that even in the fully degraded state a small residual stiffness remains.The influence of the parameter k g on the solution and its importance to the choice of discretisation are shown in Sect. 2.2 and in Sect. 2.3, respectively. The crack surface Γ is approximated by Γ l with where γ l is the crack density function. One possible choice for γ l is γ l (φ, ∇φ) := 1 2l φ 2 + l 2 ∇φ 2 (12) as proposed in [47] and in a similar way much earlier in [13].
The crack density function depends on the phase-field variable φ and its gradient ∇φ. The internal length parameter l governs the width of the transition zone between broken and unbroken state. The combination of this definition with the variation of the total energy of the phase-field problem (cf. Eq. (13)) leads to an inhomogeneous Helmholtz-type equation reminiscent of gradient enhanced damage (GED) formulations [54]. However, both approaches (PFM and GED) differ regarding their motivation and regarding the respective crack-driving source term. The PFM is motivated from an energetic approach and, therefore, the strain energy density serves as crack driving parameter. A comparison between the two methods is presented in [1,65]. Inserting the individual parts into Eq. (1) the total energy E l of the regularised problem reads (13) This total energy E l is minimised with respect to the primary variables u and φ. The additional penalty potential Π irr with prevents healing of a fully developed crack.Γ are the surfaces in the domain, where φ ≈ 1.0 applies. Here, irr is a small penalty parameter [30]. To inhibit non-physical crack propagation under compression, only certain parts of the strain energy density function should be degraded. Different splits of the strain energy density function have been developed, e.g. the volumetric-deviatoric split [5] and the spectral split [46,47]. In this contribution, the spectral split is used in combination with the hybrid model by Ambati

One-dimensional considerations
Miehe et al. [47] examine an infinitely expanded bar Ω = A × L with L = [−∞, ∞] where A represents the crosssection, L the length of the bar and x indicates the current position along the axis of the bar. The regularised crack sur-face Γ l of this one-dimensional bar reads with dv = A dx. For the fully broken state with a crack located at x = 0 (i.e. φ(0) = 1 and φ(±∞) = 0) Γ l = A must hold. For this case, the solution can be identified. This was derived in a similar way in [10] while taking into account the mechanical part of the problem and in [37] for a slightly reformulated definition of the crack surface density function. To the best of our knowledge, for all cases with only a partially developed phase-field with 0 < φ(x = 0) < 1 and φ(x = ±∞) = 0, no direct analytical solution for φ(x) could be found yet. Borden et al. [10] deduced an analytical solution, but had to evaluate the arising integrals numerically. They state, that only for the fully broken case (φ(x = 0) = 1) a kink occurs in φ(x = 0), for all other cases (i.e. the crack has not developed completely) the function φ(x) remains smooth. The degraded strain energy of a linear elastic bar reads The variation of this equation with respect to the displacements u results in Assuming φ(x) = exp (−|x|/l), k g = 0 and the boundary conditions u(−∞) = 0 and u(∞) =ū, the solution for the displacements yields where H (x) represents the Heaviside step function. If φ(x = 0) < 1, the material is not degraded completely at x = 0. In order to still fulfil Eq. (18), high strains u ,x correspond to the maximum values of φ(x). As a consequence, the shape of the displacement function u(x) is a smooth sigmoid function which converges to the Heaviside function in the fully broken state. To show this effect, different predefined phase-field functions representing different states in the process of a developing crack were applied to a one-dimensional bar with L = [−1, 1]. The regularising internal length was chosen to be l = 2. The mechanical part of the problem (Eq. (18)) was solved via FE-simulation for the displacement field u(x) across the crack, with the applied boundary If k g is chosen to be k g = 0, even for a fully developed phase-field, a small residual stiffness remains. Therefore, the displacement field u(x) will obtain the shape of a sigmoid function, the slope of which depends on the chosen parameter k g . This is shown in Fig. 2 for the same FE-simulation.

Discretisation aspects
Equation (13) is solved within the framework of the finite element method (FEM). For classical finite elements, the solution for the primary fields φ and u as well as their gradients can only be reproduced adequately using standard polynomial ansatz functions if the mesh size is sufficiently small or the ansatz order appropriately high. As shown in Sect. 2.2, high gradients along with non-polynomial and nonsmooth solutions will occur for u and φ. As an upper bound for the element size h, Miehe et al. [47] suggests l/h > 2 for first order quadrilateral elements in order to minimise the approximation error of the phase-field. Linse et al. [41]  note, that Miehe et al. only took the phase-field approximation and its gradient into account and neglected the coupled mechanical problem in their considerations. Therefore, they suggested an even higher threshold for adequate l/h-ratios than given in [47]. By means of a one-dimensional bar, it is shown in [44] that there exists a discrepancy between the theoretically proved Γ -convergence and numerical results if l → 0. They state that the ratio l/h should lie within a specific range in order to minimise the error. This range depends on the specific example. In [41] even ratios l/h ≥ 32 are applied to a onedimensional problem, where quadratic ansatz functions are employed. From a practical point of view, Wu [65] suggest a ratio from l/h ≥ 5 to l/h ≥ 10 for two-dimensional problems, which is also adopted in [43].
However, in studies addressing the l/h-ratio, the dependency on the chosen k g is not yet investigated. In a onedimensional setting, the parameter k g influences the slope of the sigmoid function the displacement field is tending to for a growing phase-field. Therefore, k g exhibits regularising properties. Choosing an l/h-ratio which is too large for a given k g leads to an overly stiff behaviour. Hence, for the standard phase-field approach, the mesh size h needs to be chosen not only dependent on the length scale l but also in relation to the artificial stiffness k g . This becomes evident when solving the mechanical part of the problem (Eq. (18)) for a fully developed, predefined phase-field, such as described in Eq. (16), with l = 2 and a chosen k g = 10 −5 again over a one-dimensional bar (cf. Sect. 2.2) for different meshes (10, 30 and 50 elements). Again, linear shape functions were used. The respective resulting displacement functions are depicted in Fig. 3. This overly stiff behaviour leads to a spreading of the phase-field in the coupled phase-field problem until at least one element is fully degraded.
Furthermore, independent of l and for k g = 0, the element size must fulfil h → 0 to be able to represent the displacement jump across the crack exactly.
Bourdin et al. [14] suggest a correction of the critical fracture energy release rate G c depending on the l/h-ratio readingḠ c = G c (1 + h/(4l)) to handle the localisation error due to the finite element dimensions. A detailed study on this effective critical fracture energy release rateḠ c can be found in [11]. Among others, this corrected energy release rate has been investigated in [44,51,61]. In this contribution such a correction is not necessary, since the presented numerical framework can handle localisation independently of the mesh.
In [43] the influence of a uniform mesh and the length scale l on the mesh bias of a developing crack is investigated. It is observed, that for undersized l/h-ratios, the phase-field profile follows the element edges. Adjusting the ratio to l/h = 15 the sensitivity towards the mesh bias vanishes. We believe, the authors neglect one important reason which is in accordance with their observations: in order to represent the discontinuity of the displacement field which appears if no stiffness remains along the crack (k g = 0) for low order elements, at least one fully degraded element is necessary. Otherwise, the displacement jump across the crack can not be approximated adequately. To illustrate this, a tension-test of a notched plate in a two-dimensional setup is examined (boundary conditions in Fig. 31a). The fully developed phase-field and the corresponding displacement field in y-direction is shown in Fig. 4 on the left. For the example, an internal length scale parameter of l = 0.008 mm has been applied. A detailed view of the phase-field and the displacement field for three different element sizes (h/l = 1, l/h = 2 and l/h = 5) is given in Fig. 4 on the right. The mesh, consisting of second order 6-node triangular elements, is regular but not parallel to the expected crack path (bias 7 • ). It can be seen that for an increasing l/h-ratio, the sensitivity towards the mesh size vanishes. The finer the mesh, the less a zigzag-pattern of the phase-field due to the mesh is observable on a global scale. Nevertheless, regarding the displacement field even for a bigger l/h-ratio, the crack still follows the mesh geometry. Also, the mesh orientation affects the crack position. Only with increasing l/h-ratio, the crack can be reproduced at the correct position. Without special treatment, it is not possible to capture the approximated jump in the displacement field at the correct position independently of the mesh.
In the context of the finite difference method, Finel et al. [23] consider a rather similar problem for a phase-field capturing the transition zone between two different material phases. There, the shape of the smeared transition zone conforms to a sigmoid-type hyperbolic tangent function (tanh) which is approximated by linear functions between the grid points. It can be observed that the approximated function can not develop independently of the position of the grid points of the underlying grid. In order to overcome this so-called grid pinning, the authors propose a transformed ansatz, which can reproduce a hyperbolic tangent function at any position of the grid independently of the grid point positions. Instead of a minimum of eight grid points necessary for the standard formulation, only one grid point is enough to capture the tanh-profile adequately. With this sharp phase-field method (SPFM) the number of grid points necessary can be reduced without any significant loss of accuracy.

Extended phase-field method (XPFM)
Motivated by the reduction of the computational cost of phase-field simulations of brittle fracture by reducing the l/h-ratio without loss of accuracy and with insensitivity to mesh bias, the new concept of the XPFM is presented in the following. It adapts concepts from the well known XFEM and the newly developed SPFM.

Derivation of the XPFM in 1D
For the derivation of the XPFM, it is reasonable to focus on the phase-field approximation ansatz first and to address the displacement field approximation afterwards. Particularly, the phase-field ansatz is based on the known analytical solution of the phase-field problem in 1D. Then again, the displacement field ansatz is based on the phase-field approximation.

Approximation of the phase-field
The FEM is able to reproduce every function which is included in the ansatz space exactly. Usually, linear or quadratic polynomial shape functions are applied. These are solely able to approximate the exponential shape of the phase-field if h → 0. Within one second order element, an arbitrary quadratic function f can be defined as where a, b and c, {a, b, c} ∈ R, are parameters defining the shape of f . This function f can be utilised to define a transformation of the phase-field ansatz: A possible choice for ϑ, following Eq. (16) is For a fully developed crack (c → 0), Eq. (16) can be reproduced exactly: Furthermore, the limiting cases φ and its basis function f both have one extremum at the position where the crack is located. The maximum in φ always corresponds to the minimum in f . Only for the fully broken state the transformation ensures a kink in ϑ( f ) at the position of the crack, for all other states ϑ( f ) and therefore φ remains smooth.
The transformation of the phase-field ansatz (Eq. (22)) presents two challenges: Firstly, the representation of φ(x) = 0 is not unique since φ(x) → 0 if f (x) → ∞ and, secondly, if φ(x) → 1 the derivative of the root function in ϑ( f ) tends to infinity for f (x) → 0. The first problem is resolved by applying the transformed ansatz only in regions where it is required (cf. Sect. 3.3). Secondly, the root function close to zero is approximated by a polynomial function with finite derivative. To ensure C 2 -continuity between the square root function and its approximation near zero, a third order polynomial function ς is defined as where k φ is a threshold value. Employing this regularisation function for f < k φ , the transformation function becomes The stabilised root function is depicted in Fig. 6. The solid graph corresponds to the final composed, stabilised root func-

Approximation of the displacement field
Besides being able to reflect the localisation behaviour for a developing crack, the displacement field must be able to reflect rigid body modes and all other deformation modes that can be captured by a classical at least first order finite element ansatz. As a consequence, we employ an XFEM ansatz introducing an enrichment function, which is able to capture the characteristics of the displacement-field as described in Sect. 2.2. Such an enrichment function can directly be constructed utilising the basis function f of the transformed phase-field ansatz. One possible enrichment function u φ reads This function converges to the sign function for the fully broken state (c = 0): For all other cases (c > 0), u φ is a smooth sigmoid function. In Fig. 7 the enrichment function (Eq. (28)) is presented for different maximal values of the phase-field.

Categorisation of element types in 2D
When expanding the presented one-dimensional method to 2D, different element types have to be considered, since not for all elements a transformed phase-field ansatz or an enriched displacement ansatz is required. Triangular elements with six nodes each are used. They are able to reproduce arbitrary quadratic functions, which are necessary for the transformed phase-field ansatz. The following four element types are introduced: The entity of all nodes is declared I . It must hold that as well as An exemplary element type distribution, the mapping of the nodes to their corresponding nodal subset and the ramp func-

Approximation of the phase-field
In standard elements (SE) the phase-field is approximated by the quadratic interpolation where 2 N i are the six Lagrangian shape functions of a second order 6-node triangular element and φ i are the phase-field values at the nodes. In all other elements (except the blending type TSb) a transformed ansatz referred to Eq. (27) is introduced: Note that the nodal values f i do not represent the nodal phase-field values directly any more. This transformed ansatz (Eq. (33)) is only introduced in elements where φ h exceeds a certain threshold value φ t . This guarantees, that in elements with phase-field values that are (nearly) zero, f h does not tend to infinity. Between elements with the transformed ansatz and elements with the standard ansatz, blending elements with the ansatz function need to be defined. This blending technique only guarantees continuity at the nodes, but not across the edges of an element. It becomes apparent that this limitation can be neglected if the threshold value φ t is small enough (cf. Sect. 4.2). It is emphasised that additional standard FEM approximation terms are neglected in the formulated transformed phase-field ansatz, since no constant phase-field modes (comparable to rigid body modes of a displacement field approximation) have to be considered due to the nature of the underlying phase-field differential equation. Besides, additional standard FEM ansatz functions would lead to linear dependencies between the standard degrees of freedom and the degrees of freedom of the transformed ansatz.

Approximation of the displacement field
The displacement field approximation is based on the XFEM in terms of enriching the ansatz space with special enrichment functions. Two enrichment functions are introduced to capture the displacement jump behind the crack and the crack opening at the crack front adequately. The approximation u h of the displacement field reads where u i are the standard displacement degrees of freedom and {a i , b i } the additional enriched degrees of freedom. The enrichment function χ J represents the crack jump behind the crack front. The front enrichment function χ F is similar in nature to the regularised jump enrichment function χ J . They both enable the ansatz function to reproduce a sigmoid function, however with different slopes. The two-dimensional extension of Eq. (28) for the enrichment functions is not straightforward, since additional information about the crack geometry is required. The crack position and geometry in the XFEM / GFEM are typically described by level set functions [60] or explicit crack representations [21]. Within the XPFM, however, the crack position and geometry is extracted from the phase-field and is not described directly. As a result, there is no need to introduce level set functions or explicit representations of the geometry of the discontinuity.

Crack front enrichment
At the crack front, the transition between the crack jump and the intact material must be represented by the enrichment function χ F . Therefore, Eq. (28) is rewritten, resulting in Here D r f h is the directional derivative of f h in the direction of r, where r is the normalised direction perpendicular to the crack surface (cf. Sect. 3.4.4). The additional parameter k reg corresponds to the regularisation of the enrichment function behind the crack front (cf. Sect. 3.4.2). This parameter k reg can be correlated to k g , since, as described in Sect. 2.3, it leads to a similar regularisation of the displacement function. In the following, however, the residual stiffness is chosen to be k g = 0 and k reg is therefore opted to be a small, but non-zero quantity to ensure numerical stability.
To avoid (almost) linear dependencies between the standard and the enriched term in Eq. (35) a stabilisation of the enrichment function similar to the SGFEM approach [8] is introduced. Furthermore, a ramp function ρ F in the FTbelements is used to ensure the fulfilment of the partition of unity in the blending elements [26,40,42]. The nodal values ρ F i of the ramp function take the value ρ i = 1 at all nodes belonging to FE-and FJb-elements and ρ i = 0 at the corner nodes of FTb-elements which are not also nodes of an FE-or FJb-element, respectively (cf. Fig. 8). In between, the values of the ramp function are interpolated linearly. Thus, the final front enrichment function reads

Regularised jump enrichment
Theoretically, the same enrichment function χ F could also be used to represent the regularised displacement jump behind the crack front. However, the values for f h as well as for D r f h tend to zero for the fully broken state. It turns out, that is a good choice for a robust enrichment function. The small parameter k reg controls the width of the regularised zone. For k reg = 0 the modified Heaviside step function commonly used in classical XFEM [50] can be reproduced. Applying the same stabilisation technique as before, the enrichment function for the jump enriched elements (JE) reads Note, that the directional derivative of the basis function perpendicular to the crack D r f h takes on level set like properties. Furthermore, along the edges of the jump enriched elements where the crack does not intersect, the jump enrichment is almost exactly zero, due to its stabilisation. Hence, no blending is necessary in the neighbouring elements behind the crack front (element type TE). Near the crack front, however, in element type FJb the ramp function ρ J is used.

Computation of the directional derivative D r f h
The directional derivative D r f h is defined as Since the transformed phase-field ansatz retains the C 0continuity of the standard ansatz, the gradient ∇ f h is discontinuous across the element edges and thus cannot directly be used as part of the enrichment function for the displacement field. Therefore, D r f h is approximated by a function ϕ h , which must fulfil the following properties: ϕ h has to be a C 0 -continuous, e.g. piece-wise linear function, and its iso-zero line must coincide exactly with the minimum of f h (and thus with the maximum of φ h ) for the fully broken state. The first two properties can be achieved by fitting ϕ h to D r f h using linear shape functions utilising a leastsquare-fit. The third property is guaranteed by introducing a Lagrange multiplier term. The resulting problem has to be solved on the limited subdomain B E around the crack, containing all elements of types JE, FE, FTb and FJb. The unknowns ϕ h and Λ h are approximated using linear shape functions 1 N i with i = 1 . . . 3: The Lagrange multiplier term is evaluated along the crack. Within one element, the transformed phase-field ansatz is merely able to represent a straight crack. Hence, the position of the crack inside an element can be identified by just two points. Advantageously, these two points per cracked element are localised at the element edges. Therefore, the directional derivatives of f h along the three element edges are calculated. If one of these directional derivatives takes the value zero within the bounds of the investigated element and the phase-field value φ is virtually 1 at this position, then this point is defined to be a crack point x c (cf. Fig. 8).

Computation of the direction r perpendicular to the crack
In Eq. (40) the direction r perpendicular to the crack has to be determined. A purposeful approach is to use the eigenvector corresponding to the greatest positive eigenvalue of the strain tensor ε. This principle strain direction is (nearly) perpendicular to the developing crack surface at the crack front and also in front of the crack, cf. Fig. 9b. Behind the crack front, the direction r is fixed in order to maintain the perpendicularity to the crack even in case of a mode II or a mixed mode crack. It is remarked, that this choice is physically meaningful, since the greatest positive principal strain at the crack front is primarily responsible for crack opening phenomena [52]. The strain field ε is discontinuous across the element edges when applying C 0 -continuous displacement-field ansatz functions. A least-square-fit of the form in conjunction with a lumped projection strategy and a piecewise linear approximation of the projected smooth strain field ε is used to obtain a C 0 -continuous strain field. Then, at each node i the eigenvalue problem (a) Gradient ∇f h in phase-field φ.
(b) Direction r in displacement field u y .

Algorithmic aspects
Solving Eq. (13) with respect to the unknowns {u, φ}, the convergence to the correct solution within a standard Newton-Raphson-scheme can not be guaranteed. Several approaches have been developed to overcome this problem, e.g. [29,32,35,36,47]. Most commonly used is the staggered solution process, proposed in [13,46] for phase-field calculations in the context of brittle fracture. Due to its simplicity and effectiveness, this approach is utilised in this contribution.

Staggered solution process
The idea of a staggered solution process is to solve the governing equations iteratively, alternating between the equations to solve as well as the unknowns to solve for. The variation of Eq. (13) is divided into and Note that in Eq. (48) ψ + drives the phase-field instead of ψ 0 to avoid crack propagation under compression [46]. Here, the strains ε are decomposed into a tensile part ε + and a compressive part ε − . Only ε + is supposed to contribute to driving the crack growth. This split of the strains reads with n being the spatial dimension, ε i the i-th principal strain, n i the corresponding principal strain direction and the definition • ± := (• ± |•|)/2. The parts of the strain energy density which are associated to tension (+) and compression (−) are defined as To circumvent the evaluation of the derivatives of ε ± with respect to ε which would ensue from introducing the split into Eq. (46), Ambati et al. [1] propose to degrade the whole strain energy density function if ψ + ≥ ψ − , otherwise, no degradation is considered, such that In an alternating manner, are solved for u and φ, respectively. Due to the transformed phase-field ansatz, Eq. (48) is non-linear with respect to the unknown degrees of freedom f i (but remains linear with respect to φ i ) whereas Eq. (46) remains linear in the unknowns u i , a i and b i even for the enriched state. Thus, Eq. (48) is solved by means of a Newton-Raphson procedure. Typically between 2 and 4 iterations are necessary to reach convergence without crack propagation and between 6 and 12 iterations if the crack propagates. The entire solution process is depicted in Fig. 10. The number of staggered iterations n stag as well as the load step size Δt are fixed and defined explicitly. Within every staggered iteration, first Eq. (52a) is solved for the displacement field u h and then Eq. (52b) is solved for the phase-field φ h . Subsequently, the enrichment functions u φ are recalculated from the phase-field values and the eigenvectors of the strains. If an adjustment of the enrichment scheme is detected (i.e. an adjustment of the element types), the current staggered iteration is repeated with the new enrichment scheme until no adjustment of the enrichment scheme is necessary anymore.

Update of the enrichment scheme
The enrichment scheme has to be updated if a crack propagates. If the phase-field φ h exceeds a certain threshold value φ t at an integration point, the element type is changed from standard (SE) to transformed (TE) type. Thereby, the interpretation of the phase-field degrees of freedom change. Hence, the degrees of freedom have to be shifted by if a node i ∈ I S is transformed to a node i ∈ I T . A reasonable threshold value φ t is in the range of φ t = 0.1 . . . 0.2. This choice prevents algorithmic difficulties within the blending elements of type TSb. Elements with a phase-field value φ h > φ e (φ e > φ t , e.g. φ e = 0.4 . . . 0.5) at an integration point are considered as candidates for displacement enriched elements (type JE or FE). In these elements, Eq. (41) is solved to compute the fitted directional derivative ϕ h for the enrichment functions u F φ and u J φ . Thereby, the positions of the crack points are determined as described in Sect. 3.4.3. In elements with exactly two crack points, the regularised jump enrichment function u J φ is used (JE). Elements with one or less crack points and non-constant enrichment function u F φ are regarded as front enriched elements (FE). Following this, the blending scheme and the ramp functions ρ F and ρ J are defined as depicted in Fig. 8.

Damped Newton-Raphson-procedure
For the calculation of φ n+1 i+1 a damped Newton-Raphson method including a line search algorithm is used to avoid nodal phase-field values φ i > 1 during the iteration process. Therefore, the updated solution vector φ n+1 i+1,k+1 = φ n+1 i+1,k + Δφ is checked for negative values f i (which would lead to nodal phase-field values φ i > 1). For those entries, a factor Δφ i being the change of the degree of freedom of the considered node. The increment Δφ is reduced by the factor The parameter α min can be interpreted as a line-search parameter which controls the increment in such a way, that no negative values f i occur.

Quadrature aspects and parametric study
The accuracy and numerical efficiency of the presented method depend on the choice of the stabilisation parameters k φ and k reg , on the l/h-ratio and, furthermore, on the choice of integration scheme. It can be seen from Sect. 3.1.1, that for a parameter k φ → 0 the phase-field tends to the expected solution for fully developed cracks, but, as shown in the following section, a small k φ leads to an increased numerical effort due to the need for a more accurate quadrature rule. Furthermore, as noted in Sect. 2.3, for standard PFMsimulations the l/h-ratio must be sufficiently large to enable an accurate approximation of the phase-field, the displacement field and their respective gradients. For the XPFM, the l/h-ratio can be chosen much smaller. However, the reproducible width of the phase-field and thus the applicability of the method for small l/h-ratios certainly depends on the numerical quadrature scheme, which can lead to a quite high computational effort on the element level in a few elements and in extreme cases. Naturally, the aim is to keep the computational effort low while maintaining accurate solutions. In the following, several studies for the right choice of the different parameters are conducted and presented. Furthermore, the choice of the numerical quadrature scheme is discussed and assessed.

Preliminary Investigation in 1D
To get a good understanding of the influence of the parameters k φ , l and k reg on efficiency and accuracy, two different aspects are preliminarily analysed in one dimension. Firstly, the error of the crack surface integral computed with the stabilised phase-field solution compared to the analytical solution is studied. Secondly, a usable quadrature scheme is presented and investigated for the crack density function γ ς l depending on the phase-field approximation and for the strain energy density function ψ 0 depending on the displacement field.

Error of the stabilised crack surface
The analytical solution φ (Eq. (16)) of the phase-field problem (Eq. (18)) and its stabilised counterpart φ ς (Eq. (27)) as depicted in Fig. 11 are employed to calculate the regularised crack surface. Note, that both functions are identical apart from the stabilised area around the phase-field peak (magnified) whose size depends on the parameter k φ (cf. also Fig. 6).  Fig. 12. The integral of γ l is the regularised crack surface Γ l and the integral of γ ς l , Γ ς l , should be a close approximation depending on k φ . The smoothing of the phase-field tip is necessary to make the transformed phase-field ansatz differentiable at the crack. However, this leads to significant deviations of the approximated crack density function γ ς l from γ l . These deviations originate in the gradient of the stabilised phase-field function, which is incorporated in the crack density function. With decreasing k φ , the deviations increase and tighten, which corresponds to a smaller error of the crack surface compared to the analytical solution of the phase-field but also to increasing difficulties concerning accurate numerical quadrature.
Taking into account that Γ l = A, the error of the stabilised crack surface in comparison to the analytical solution for different stabilisation parameter k φ (blue) and length parameter l (orange) is depicted in Fig. 13. It becomes obvious that e Γ l is solely dependent on k φ and independent of l. Furthermore, k φ ≤ 0.006 seems to be a natural choice to ensure an error of less than 1% for 1D problems.

Integration strategy
The integration within a one-dimensional reference element with the reference coordinate −1 mm ≤ ξ ≤ 1 mm is investigated. The phase-field and therefore the crack density function is prescribed in the whole element with the maximum of φ max at the location ξ = ξ max with a chosen basis function f (cf. Eq. (20)) with In Fig. 14 the stabilised functions φ ς (blue) and γ ς l (orange) are depicted firstly for a crack in the middle of the element (ξ max = 0 mm, solid lines) and secondly (dashed lines) with φ max = 1.0 (if not denoted explicitly different) in a pseudorandom location ξ max ≈ 0.3351 mm where it is ensured, that no integration point position coincides exactly with ξ max .
To reproduce the crack surface Γ ς l satisfyingly, a suitable quadrature scheme is necessary. The standard Gauss-Legendre-quadrature for polynomials is not very accurate, since depending on the position of the crack the underlying function is piece-wise defined and, in general, the location of the crack is not known a priori. Therefore, a subdivision of elements is used to be able to capture the function more accurately. The integration is then carried out over each subelement. A regular subdivision in evenly-sized elements seems to be the best choice to account for the unknown crack position. Aside from that, it should be noted, that this kind of subdivision does not imply an increase in degrees of  Fig. 15 Quadrature schemes in one dimension freedom, since the subelements are solely used for integration purposes. Other choices of integration are an interesting subject of further research due to the great potential of reducing numerical effort and increasing accuracy of the XPFM even more. In Fig. 15 the five tested integration schemes can be seen. Each subelement contains one, two, three, four, or six Gauss-integration points which allow for a different order polynomial Gauss-Legendre-integration within that subelement. The accuracy of numerical quadrature is evaluated by the means of the error e int  Fig. 16 (top) the influence of the position of integration points becomes obvious due to the widely different results for even numbers of integration points compared to odd numbers, which only congregate with increasing number. This behaviour is less obvious in Fig. 16 (bottom), but still existent.
To be able to make reliable assertions without knowing the position of the crack a priori, the envelope curves depicted in the diagrams serve as a threshold for the highest quadrature error possible, depending on the number of integration points. The fact that, as shown in Fig. 17, the envelope curves match for the different crack positions with increasing n i p , indicates that these curves serve as a reliable threshold independent of the crack position. For clarity and comprehensibility, in the following graphs solely the envelope curves are depicted for comparison.
Furthermore, it can be seen that the necessary number of subelements and therefore, the necessary number of quadrature points, is strongly dependent on the choice of k φ . By choosing a smaller k φ , the gain of accuracy is negligible (less than 0.5%, cf. Fig. 13) compared to the increase of computational effort due to the higher number of integration points (approximately factor two, cf. Fig. 17).
In Fig. 18 a similar comparison of envelope curves for k φ = 0.002 and different l/h-ratios shows a strong dependence of the quadrature error on the l/h-ratio as well, and it demonstrates the requirement to carefully choose the In Fig. 19 the integration errors for not yet fully developed cracks are compared. The l/h-ratio is chosen to be l/h = 1 and the stabilisation parameter k φ = 0.002. It can be seen that the necessary number of quadrature points and thus the computational effort reduces drastically for not yet fully developed cracks. Therefore, in the FE-simulation, it might be advantageous to gradually introduce integration points dependent on the phase-field growth. The comparison of the different integration strategies of the subelements (cf. Fig. 15) is shown in Fig. 20. For the centred crack as well as for the shifted crack, the quadrature scheme with one Gauss-integration point per subelement seems to be the most effective one, since compared to the alternative strategies the reliable error for this scheme decreases most rapidly with regard to the total number of quadrature points. This behaviour has been observed for different parameter k φ and l/h-ratios as well. Therefore, for the integration of the crack density function with the stabilised exponential ansatz (cf. Eq. (27)) equidistant and equally weighted integration points seem to be the best choice as they lead to the most accurate results with the lowest computational effort. Beyond that, this choice of quadrature scheme is advantageous with respect to simplicity and independence of the position of the crack. It also shows advantages for the identification of enriched elements by evaluating the phasefield value at each integration point (cf. Sect. 4.2).
In addition to investigating the integration of the crack density function, the integration of the driving force, the strain energy density function (cf. Eq. (17)), is analysed. As discussed in Sect. 3.4 for a fully developed crack, a Heaviside-function is introduced for the displacement field. For this case, an integration as presented in [63] would be suitable. For developing cracks a possible displacement function u(ξ ) and the corresponding strain energy density function ψ 0 incorporating the enrichment function derived from the phase-field are given by with f again being the basis function. In Fig. 21 both functions are depicted for the centered (solid lines) and the shifted (dashed lines) crack with φ max = 0.9.
with ξ i being the position of each integration point and w i being its respective weight. The l/h-ratio plays an important role for the slope of the displacement curve and thus for the width of ψ 0 which is crucial for the accuracy of the quadrature. As shown in Fig. 22 the computational effort for the numerical integration reduces drastically with increasing l/h-ratio.
Furthermore, it can be seen in Fig. 23 that, similar to the quadrature of the crack density function, the computational effort for the quadrature of the strain energy density function reduces significantly for lower φ max .
The comparison of the different quadrature schemes shown in Fig. 24 indicates that also for the mechanical part In general, the computational effort for the accurate integration of the strain energy density function is significantly lower compared to the effort required for the integration of the crack density function. Since crack propagation is driven by the strain energy density, it is especially important to place a sufficient number of integration points on the developing crack. This, however, is always ensured by the quadrature scheme chosen for the accurate integration of the crack density function.

Integration Scheme in 2D
Extending the integration strategy discussed in Sect. 5.1.2 to two dimensions is straightforward. The triangular reference element is subdivided again into congruent subtriangles which are provided with one integration point each. A sketch can be found in Fig. 25. This integration scheme, however, is only applied in the elements with the transformed phase-field ansatz, since in SE and TSb only a quadratic ansatz function needs to be integrated for which a three-point Gauss-integration is sufficient. To evaluate this integration scheme for different crack positions within one finite element, the plate depicted in Fig. 26a is tested for different subtriangulation orders. The crack is enforced in different positions via a penalty formulation which ensures φ(y = y c ) = 1.0. In Fig. 26b the respective crack position is shown in detail. The characteristic length is chosen to be l = 0.01 mm, which results in the ratio l/h = 1. Furthermore, k φ = 0.01 is chosen. The theoretical crack surface should match the width of the plate and should therefore be Γ = 0.1 mm.
The numerically attained crack surface Γ h l is once again evaluated by the error defined as In Fig. 27 e num Γ is plotted over the number of subtriangles within one finite element for the different crack positions. Most noticeable are the results for the crack positioned right on the element edge. The ansatz is able to reproduce the kink in the fully developed phase-field and can therefore reproduce the correct phase-field function better, than the stabilised ansatz. For the crack positions shifted more to the middle of the element and due to the restrictions of accuracy of the stabilised phase-field ansatz dependent on k φ , the error is not able to reduce beyond a certain threshold value, even with extremely accurate integration. This can be seen by the convergence of all the envelope curves to a constant value (except for y c = 0). Nevertheless, it can be seen that the error can be reduced to less than 1% with a sufficient number of integration points.

Investigation of stabilisation parameter k
The simulation as set up in Fig. 26 is tested again for different k φ and a crack position with y c = 0.004 mm. To ensure adequate integration accuracy in each case, 80 subtriangles in each direction are used. In Fig. 28 the error e num Γ (blue, cf. Eq. (61)) is plotted over the different k φ . It can be seen that the error is below 1% for k φ = 0.01 in this case. Apart from this, the curve is comparable to the related 1D investigation depicted in Fig. 13. Additionally, the error of the numerical solution compared to the stabilised analytical solution increases. This is due to the fact, that the stabilisation area has become too large and φ ς deviates too much from φ to be the best approximation the given ansatz is able to reproduce.

Convergence to the correct crack length and L 2 -error-convergence
To assess the properties of the XPFM in comparison with the classical PFM, the convergence to the correct crack length and the quality of the displacement-field approximation is examined for the plate under tension depicted in Fig. 26a. A predefined crack (by enforcing φ h = 1.0 at the crack) cuts the plate transversely (inclination of 11,3 • to the horizontal). The analytical solution of the crack surface Γ and the displacement field u is known for this case. For a constant internal length parameter l = 0.01 mm and a prescribed displacementū y = 5 · 10 −5 mm at the top of the plate the crack surface Γ h l and the L 2 -error of the displacements u h for different l/h-ratios (variable element size h of the regular mesh), defined as are evaluated numerically for both, the standard phase-field approach and the presented extended approach. The results are displayed in Fig. 29. In case of the XPFM simulation, a maximum number of 80 integration points per direction have been used in the required elements. Moreover, the parameters k φ = 5 · 10 −3 and k reg = 4 · 10 −7 were applied. Figure 29 reveals that the standard phase-field method is highly sensitive to the l/h-ratio. It can be seen that the crack  Fig. 29 surface Γ can not be reproduced adequately for l/h < 5. However, the XPFM does not show such a dependency on the element size for an arbitrary positioned crack. Furthermore, the L 2 -error in the displacements is constant (and thereby independent) of the element size for the XPFM as well, whereas an appropriate approximation of the displacements in the standard phase-field approach requires a large l/h-ratio. The quality of the approximation of the displacements within the XPFM only depends on the regularisation parameter k reg . For a varying k reg , the L 2 -error of the displacements is displayed in Fig. 30.ū

Crack propagation simulations
One of the advantages of the PFM for fracture problems is that crack propagation is calculated automatically, i.e. unlike for the XFEM / GFEM it is not necessary to evaluate a crack propagation criterion, a crack propagation direction and a crack increment length. In the following, the newly developed XPFM is compared to the standard PFM for the crack propagation examples of a standard tension test (mode I) and a standard shear test (mode II) of a square plate to assess whether the XPFM is able to capture the expected forcedisplacement curves as well as the expected crack paths. For both tests the stabilisation and regularisation parameters are set to k φ = 5 · 10 −3 and k reg = 2 · 10 −4 · l 2 .

Mode I: plate under tension
The plate under tension depicted in Fig. 31a is tested first. A linear elastic material for small strains with Lamé-constants λ = 121.15 kN/mm 2 and μ = 80.77 kN/mm 2 , a critical energy release rate of G c = 2.7·10 −3 kN/mm and an internal length of l = 0.008 mm is chosen.

Behaviour of the standard PFM approach
In order to contextualise the results of the developed method, results of the standard phase-field approach are presented first. A staggered solution process with a fixed number of staggered iterations n stag per load step is applied. The load stepping scheme is taken from [46]. For the first 500 load steps, an increment of Δu = 10 −5 mm of the prescribed displacement is applied. After that, the increment is reduced to Δu = 10 −6 mm until failure. The initial crack is modelled explicitly. Within the staggered solution process, no global convergence criteria is evaluated. This leads to slightly different results for different number of staggered iterations n stag and different load steps sizes [1]. For n stag = 8 the mode I It can be seen that with increasing ratio the curves converge to the same solution. The choice of a ratio of l/h ≈ 5 . . . 10 seems to be sufficient to obtain accurate results. For a small number of staggered iterations n stag convergence is not achieved. This can be seen in the forcedisplacement curves depicted in Fig. 33 for a ratio of l/h = 10. It becomes obvious, that the number of staggered iterations influences visibly how steep the curve declines after exceeding the peak load.
In the following, the results of the XPFM simulations are compared to a standard PFM simulation with l/h = 10 in the critical region and n stag = 8 to ensure a comparison to a (nearly) converged numerical solution.

Comparison of results obtained with the XPFM and the PFM approach
The mode I tension simulation using the XPFM is carried out on two meshes, obtaining a ratio of l/h = {1, 2}. Due to the brittle material behaviour and the abrupt load drop after exceeding the peak load, a Γ l -controlled adaptive load stepping is applied. A minimal possible displacement increment of Δu min = 10 −12 mm is allowed. Therefore, and by reason of the additional enrichment loop, only one staggered iteration n stag per load step is sufficient. Figure 34 shows the load displacement curves compared to the reference solution computed using the standard PFM with l/h = 10. The curves converge rapidly, i.e. the XPFM leads to accurate results even for quite coarse meshes. The resulting peak loads are summarised and compared to standard PFM simulations in Table 1. It can be seen that the error decreases to the reference solution for the PFM as well as for the XPFM with decreasing element size. However, even for the coarsest chosen mesh with a ratio of l/h = 1 the error of the XPFM solution is already less than 1%. This means that using the XPFM, sufficiently accurate solutions can be obtained at a significantly lower computational cost compared to the standard PFM.
The insensitivity against the mesh size and mesh orientation of the presented approach is demonstrated in Fig. 35. For the crack state directly before final rupture, three crack details are shown: the starting point of the crack, the middle of the crack and the crack tip. Independent of the position of the nodes and the orientation of elements, the phase-field can develop in the whole domain. The (slightly regularised) jump in the displacement field always corresponds to the phase-field and is also independent of the mesh. Due to the stabilisation (cf. Eqs. (39), (37)) the crack can also cut elements close or even through nodes. At the crack tip, the sigmoid-like enrichment function interpolates between the fully broken and the intact material.

Crack initiation in comparison to the analytical solution
For an ideal traction-controlled mode I tension test, the analytical solution of the critical external load σ c leading to unstable crack propagation is known. For plane strain conditions, the relation between the critical energy release rate G c and σ c is with a and Y being parameters depending on the geometry. For the geometry depicted in Fig. 36a the crack length a is a = 0.1 mm and the factor Y is Y = 1.1957 [6]. Figure 36b shows the analytical solution for the critical stress σ c for different critical energy release rates G c as well as the XPFM solution for the stress at which crack propagation occurs. It can be seen that the numerical results obtained with XPFM are able to reproduce the analytical results quite accurately.

Mode II: plate under shear
To demonstrate that the XPFM is able to predict the expected crack path, a mode II example is simulated. Here, the same material parameters as for the mode I tension test are used. The boundary conditions are depicted in Fig. 31b. A quasistatic displacement increment is applied with 140 steps and Δu = 5 · 10 −5 mm for each load step. After that, the displacement increment is set to Δu = 1·10 −6 mm until failure. To suitably compare the PFM approach and the XPFM approach, one staggered iteration per load step (n stag = 1) is applied in both cases. The standard PFM simulation is carried out on a mesh with l/h = 5 in the critical region, compared to l/h = {1, 2} for the XPFM simulation. As can be seen from the load-deflection curves in Fig. 37 the XPFM can reproduce the expected curve, even for a mesh with l/h = 1.
Besides the final phase-field φ and the corresponding displacements u x in x-direction of the XPFM simulation with l/h = 1, Fig. 38 shows detailed views of the crack for the XPFM and the PFM simulation. It can be seen that the resolution of the mesh for the PFM simulation is still not sufficiently fine to allow for a smooth crack geometry approximation. The crack follows the mesh in a zig-zag-pattern. In contrast, the XPFM is able to capture the smooth crack path independently of the mesh orientation and node positions.  Both examples demonstrate the great advantage of the developed XPFM compared to the standard PFM. The crack can be represented entirely independent of the mesh orientation. Much coarser meshes can be used without any need of an a priori or adaptive mesh refinement in regions where the crack is expected to grow. In two dimensions an XPFM simulation using a mesh with l/h = 1 compared to a comparably accurate PFM simulation using a mesh with l/h = 5 reduces the number of degrees of freedom approximately by the factor of 25. For the extension to three dimensions, this factor is expected to be approximately 125. Despite the higher number of iterations due to the nonlinear transformed phase-field approach, a huge potential in reducing the computational effort, especially in three dimensions, can be seen.

Conclusions and outlook
In this contribution, an extended approach to the phase-field method for brittle fracture was proposed. It is based on a transformation of the phase-field ansatz and an enrichment strategy of the displacement field ansatz. The method offers the possibility to resolve the phase-field profile and the corresponding jump in the displacement-field on much coarser meshes compared to standard phase-field approaches. Good results can already be obtained for a ratio of l/h = 1. A possible robust staggered solution algorithm with enrichment scheme update and an effective integration scheme were introduced and tested. Furthermore, numerous tests were performed to give an insight into the different aspects of the method with regard to parameter setting and to enable the comparison to approaches and results described in literature. The carried-out simulations exhibit the great potential of this method due to the tremendously reduced computational effort compared to classical phase-field approaches. Furthermore, the extension to three dimensions is expected to be straightforward.
Nevertheless, since this contribution serves solely as an introduction into the XPFM, possible algorithmic improvements and extensions to the method are numerous such as the application to ductile fracture, dynamic fracture and fatigue processes. With this method, the global number of degrees of freedom and thus the main computational effort can be reduced enormously, which is especially advantageous looking ahead at 3D simulations. However, the computational effort on the element level increases compared to the standard phase-field method, since the currently used integration scheme is much more cost intensive. This still offers room for improvement, which can increase the efficiency of the approach even more. Furthermore, the research about convergence criteria and auto-time stepping of phase-field problems in general and of the extended approach in particular poses an interesting starting point for more in-depth investigations.