Energetically motivated crack orientation vector for phase-field fracture with a directional split

The realistic approximation of structural behavior in a post fracture state by the phase-field method requires information about the spatial orientation of the crack surface at the material point level. For the directional phase-field split, this orientation is specified by the crack orientation vector, that is defined perpendicular to the crack surface. An alternative approach to the determination of the orientation based on standard fracture mechanical arguments, i.e. in alignment with the direction of the largest principle tensile strain or stress, is investigated by considering the amount of dissipated strain energy density during crack evolution. In contrast to the application of gradient methods, the analytical approach enables the determination of all local maxima of strain energy density dissipation and, in consequence, the identification of the global maximum, that is assumed to govern the orientation of an evolving crack. Furthermore, the evaluation of the local maxima provides a novel aspect in the discussion of the phenomenon of crack branching. As the directional split differentiates into crack driving contributions of tension and shear stresses on the crack surface, a consistent relation to Mode I and Mode II fracture is available and a mode dependent fracture toughness can be considered. Consequently, the realistic simulation of rock-like fracture is demonstrated. In addition, a numerical investigation of Γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Gamma $$\end{document}-convergence for an AT-2 type crack surface density is presented in a two-dimensional setup. For the directional split, also the issues internal locking as well as lateral phase-field evolution are addressed.


Introduction
A basic principle of structural design is the integrity and stability of load bearing components. Here, an important aspect is the prevention of crack evolution. From a continuum mechanical viewpoint, cracks are a separation of previously sound bulk material. In consequence, the material's capacity to transfer tensile stresses is lost and the capability to transmit shear stresses is reduced to interlocking and frictional effects.
From the computational point of view, the continuous representation of a discrete crack's topology, see Fig. 1a, by the phase-field, see Fig. 1d, is a rather novel approach. Previous approaches often involve the explicit representation of the crack schematically visualized in Fig. 1b, i.e. by means of an adaptive remeshing strategy and the duplication of nodes, e.g. Ingraffea and Saouma (1985), Miehe and Gürses (2007) and Ortiz and Pandolfi (1999). A comprehensive overview on the theory, implementation and analysis of this method in combination with a material force approach to evaluate dynamic crack propagation criteria is given in Özenç (2016). Another approach is denoted as extended finite Here, an enrichment of the displacement field ansatz functions is employed to account implicitly for the jump of displacements across the crack face without an explicit representation of the crack surface, see e.g. Belytschko et al. (2001) and Song and Belytschko (2009). A variational description of brittle fracture is investigated in Francfort and Marigo (1998), considering the variation of the crack domain. Through the regularization of the fracture energy, inherent mesh dependency and challenging adaptive meshing strategies are overcome. An approach for the regularization is presented in Schmidt et al. (2009) as eigenfracture. In Pandolfi and Ortiz (2012), the method is implemented into an element deletion scheme for the finite element method (FEM) called eigenerosion, see Fig. 1c. The kinematics of a crack topology is obtained via the degradation of the stiffness of individual finite elements. Recently, an eigenfracture approach based on the framework of representative crack elements is shown in  and a field crack mechanics approach in Morin and Acharya (2021), which also make use of a crack orientation vector.
Similar to eigenerosion, the phase-field method for fracture introduces a regularization for the fracture energy, which is related to the surface area of the crack where G c denotes the fracture toughness. The numerical approximation of Eq.
(1), based on the ideas of Mumford and Shah on image segmentation in Mumford and Shah (1989), is discussed in Bourdin (2007) and Bourdin et al. (2000). In essence, the crack surface density γ l is introduced in order to provide an approximation of the crack surface by The crack surface density γ l ( p) is a function of an additionally introduced field variable-the phase-field p, which lends his name for the approach in general. In analogy to standard phase-transition problems, e.g. the evolution of austenite and martensite in metals or the transition of chemicals between the different physical states solid, fluid and gas, the crack is assumed to be a different "phase" of the material. More precisely, a sound solid is associated with p = 0 and the cracked state (or phase) is represented by p = 1, see Fig. 2. The numerical treatment of discontinuities is challenging and can lead to mesh dependency, localization and convergence problems, see e.g. Bažant and Pijaudier-Cabot (1988) and Pijaudier-Cabot and Bažant (1987) for a similar problem within the theory of damage and softening. A common approach to resolve the issue is gradient enhancement, see e.g. Zreid (2018). The question, whether the phase-field approach should be considered as a special type of a damage formulation, is discussed e.g. in de Borst and Verhoosel (2016) and Steinke et al. (2017). The paper at hand is based on the developments, findings and discussions published in Steinke (2021). It aims at the further development of the phase-field method for crack approximation towards a realistic and reliable approach, not only for the prediction of crack patterns due to excessive loading, but also for the accurate simulation of the post fracture structural response. In Sect. 2, the basic theory of a phase-field method for brittle fracture is outlined. Section 3 is focused on the question, how the spatial orientation of the crack surface can be defined for the directional split. Section 4 is dedicated to numerical examples and the final section closes the paper with conclusions.
2 The phase-field approach to brittle fracture in the framework of the finite element method

Regularized crack surface
The phase-field method provides a regularized approximation l for the discrete crack in terms of the volume integral given in Eq.
(2) with a crack surface density of AT-2 type (Ambrosio and Tortorelli 1990). An additional field variable p as well as the regularization length parameter l are considered. A fundamental feature of the crack set's approximation by this elliptical functional is the continuous or smooth transition between the sound and broken regions of the solid shown in Fig. 2, that is governed by the gradient term in Eq. (3). In analogy to the gradient enhancement of damage, the regularization length parameter l governs the width of the transition zone. The relation of the regularization length to the discretization is discussed for instance in Hofacker and Miehe (2011), Linse et al. (2017) and Pandolfi et al. (2021). AT-2 is a common choice in the computational mechanics community, see e.g. Borden et al. (2012), Hofacker et al. (2009) and Kuhn and Müller (2010), due to the straight forward implementation and onset of phase-field evolution at external loading. Nevertheless, the additional numerical effort to constrain the phase-field evolution in AT-1 is rewarded by an initial elastic phase before the phase-field evolution, which is a better approximation for brittle fracture than the softening behavior of AT-2 prior to the crack, see Alessi et al. (2018) and Tanné et al. (2018) for a discussion. Furthermore, higher order functionals can be used for the approximation of materials with anisotropic fracture toughness and an enhanced approximation of l at the expense of higher order gradients, see e.g. Borden et al. (2014).
By the continuous crack approximation, an explicit representation of the actual crack topology is lost. This results in additional effort necessary, if the position of the crack tip, the location of branching or the angle of kinking are required. A common strategy to visualize the crack surfaces is their representation by isosurfaces at a threshold value for the phase-field variable p c , e.g. p c = 0.95, combined with the blanking out of elements where p > p c . Numerical approaches to evaluate the crack tip propagation speed are discussed e.g. in Steinke et al. (2016).

Crack driving force
According to Griffith, the formation of crack surface requires energy, that is dissipated from the strain energy stored in the deformation of the structure. In the initial considerations, a tensile force on the atomic structure is assumed. As soon as the distance between the atoms exceeds a critical value, the atomic bond is released and a separation of the initially continuous body, i.e. a crack, is obtained. In the framework of continuum mechanics, the counterpart to force and separation in the sense of Griffith are stresses and strains, respectively. Furthermore, the energy associated with forces acting along the separation distance has its local continuum mechanical counterpart in the strain energy density at the material point. Hence, in a similar way as γ l is the representation of the regularized crack surface l at the material point level, the strain energy density ψ ε is the local basis to obtain the phase-field analogue of the structural level quantity energy release rate G. Assuming linear elasticity, the strain energy density is written in terms of the Lamé coefficients λ and μ, the second order identity tensor 1 and the small strain tensor ε. The additive decomposition of the strain energy density into crack driving and persistent components, ψ + and ψ − , respectively, is called the phase-field split. The phase-field split allows for a definition of the specific energetic components that drive the crack evolution. Furthermore, a thermodynamic consistent formulation of the phase-field results in a strong link between the split of the strain energy and the post fracture behavior of the phase-field crack. When proposing variational fracture for the phasefield method in Francfort and Marigo (1998), Francfort and Marigo considered the total strain energy density to be available for the formation of the crack surface. In a strict sense, the proof of -convergence for the phase-field functional, i.e. ≈ l for l → 0 is available for this split only. However, a -convergence proof for eigenfracture exists in Schmidt et al. (2009), which takes a decomposition of the strain energy into account. The decomposition of the strain energy has turned out to be an essential part of phase-field and eigenfracture models in order to avoid nonphysical predictions for deformation kinematics of cracks, e.g. crack face interpenetration and lateral deformations in the postfracture state. Thus, the split implicitly defines the traction transfer through a crack and the contact state of the crack surface. The energetic split further defines the contributions to the crack driving force. With such a split, nonphysical crack propagation at uniaxial or hydrostatic pressure is overcome, which is observed in the early phase-field formulations, see Bourdin et al. (2000).
The question on the correct split has appeared with its introduction to phase-field fracture. The most common volumetric deviatoric split and spectral decomposition are adopted and modified in many phase-field approaches. An overview on energy splits for phasefield can be found in Storm et al. (2020). The violation of the crack boundary conditions, reported for instance in May et al. (2015), Schlüter (2018) and others, inspired the stress based directional split published in Steinke and Kaliske (2018). The starting point is the set of basic crack kinematics. These are defined via The key feature of the stress based directional decomposition is the information about the orientation of the crack surface. This information is provided in a unique way by the crack orientation vector r perpendicular to the crack surface, see Fig. 4. In combination with vectors s and t, both perpendicular to r and to each other, a Cartesian crack coordinate system (CCS) is established. The decomposition of a linear elastic ground stress tensor based on Eq. (4) with respect to the CCS is visualized in Fig. 4. In order to approximate the crack kinematics visualized in Fig. 3, the crack driving stress tensor includes the shear components on the crack surface as well as the stress components perpendicular to the crack surface with tensile magnitude by • ± = (• ± | • |)/2. Furthermore, the continuum description of the crack kinematics originating from a discrete crack behavior requires an additional correction term to account for Poisson effects during crack opening, see Steinke and Kaliske (2018). It can be shown, that Eq. (6) is identical to the recent development of a representative crack element theory in case of isotropic linear elasticity at small strains, see Storm et al. (2020). The persistent stress tensor, being the counterpart to the crack driving stress tensor, reads containing the remaining components of the total stress tensor as well as the correction term with negative sign. This ensures, that linear elastic behavior is recovered in case of no phase-field crack is present. In reverse manner, the crack driving and persistent strain energy densities are obtained from the stress and strain tensors by This ensures, that the amount of dissipated energy is directly linked to the post fracture behavior of the phase-field crack, which results in a thermodynamic consistent formulation.
With the measures of the crack surface as well as the crack driving strain energy in terms of their densities at hand, the energetic link to describe the dissipation of the latter into the evolution of the former is provided by the total energy of the solid domain , where G c denotes the fracture toughness and the external forces are omitted for the sake of simplicity. This leads to the governing partial differential equations for the two-field problem, a standard quadratic degradation function is used in the following.

Crack irreversibility
From the macroscopic point of view, crack formation is an irreversible process and the irreversibility of fracture is a basic assumption in linear elastic fracture mechanics. A fundamental challenge of macroscopic irreversibility in the context of the phase-field method for fracture is the continuous approximation of the crack surface itself. While irreversibility of a discrete crack can be formulated aṡ in a straight forward manner, the question remains, whether this results in a similar constrainṫ on the phase-field degree of freedom as well as postulated in . Alternatively, a restriction on the global evolution of the regularized surfacė could be required. It should be noted, that fulfilling Eq. (14) implicitly results in Eq. (15), but not vice versa. Furthermore, both conditions are not a strict representation of the original constraint of Eq. (13). The implementation of the constraint Eq. (14) can be based on a history variable H, that is used in Eq. (11) instead of the crack driving strain energy density ψ + . While this approach is common in the theory of damage, its application in the context of the phase-field method for fracture is initially introduced by . The implementation of this approach is straight forward, as the history variable can be formulated in terms of the material point level value of the crack driving strain energy density ψ + . However, the analysis of Linse et al. (2017) indicates, that the amount of crack surface area is overestimated by this method.
An alternative approach is published by Bourdin et al. (2000). They propose an additional boundary condition on the phase-field degree of freedom, as soon as a critical value is exceeded. In consequence, the local "healing" of the phase-field value is allowed and the phase-field profile at the end of the evolution closely resembles the analytical solution, see Linse et al. (2017). Therefore, the final amount of crack surface obtained in their setup is very close to the expected value. However, this result can be obtained only by the violation of Eq. (14). In consequence, the energetic counterpart of the crack surface, i.e. the strain energy, that is dissipated locally for the higher phase-field value, can increase due to the increase of the degradation function for the decreasing phase-field value. From the energetic point of view, strain energy is created from crack surface locally. While the global energy balance is not violated in this case due to the equivalence of energetic contributions from strain and crack surface in Eq. (9), a different problem arises. The energy is a scalar value without direction. In contrast, the degradation function affects a second order tensor in the definition of the effective stress tensor in Eq. (10). As the crack evolution is reversed locally, the energy of the crack surface, that is without direction, is transformed into an increase in the effective stress tensor. Let us assume a volumetric deviatoric split. Furthermore, assume volumetric expansion triggered the phase-field evolution in the first place. Now, local healing permits the transformation of this initial "volumetric" strain energy into deviatoric components. For a directional split, the same applies to the interchange of tension and shear as well as a change of the crack orientation during evolution. However, a study on this issue is a non-trivial challenge, that may have not been tackled yet, as the phase-field community seems to be totally unaware of this hardly realistic process. Nevertheless, this issue is not investigated in this paper any further.
In this work, the combination of both approaches, that is already published in Steinke and Kaliske (2018), is applied in order to benefit from the elegant simplicity of a formulation with a damage-like history while obtaining a proper profile for fully evolved cracks. A critical value for the phase-field p c is defined. Based on this value, a history variable H(t) is introduced, that depends on the current solution time t n as well as on a previous step of the simulation at time t n−1 . The history variable reads and is used in the modified evolution equation for the phase-field

Finite element aspects
The solution of multifield finite element equations systems based on the governing equations Eqs. (10) and (17) is obtained within the parallel version of FEAP 8.5 (parFEAP) documented in Taylor (2017). ParFEAP's solution is based on a METIS Karypis and Kumar (1998) partitioning of the discretization and the application of parallel solvers provided by PETSC Balay et al. (2019). Within this paper, the generalized minimal residual algorithm published in Saad and Schultz (1986) with a Jacobi preconditioner and a modified conjugate gradient approach published in van der Vorst (1992) are applied. A multifield equation system can be solved in a monolithic or staggered scheme. In the monolithic scheme, the increments of all DOF are obtained simultaneously. Therefore, the stiffness coupling terms are of major importance. The staggered solution scheme is based on an alternated solution for specific DOF. In terms of the phase-field method, the increments of the displacements are obtained at frozen phase-field values and the increments of the phase-field are computed for a frozen state of deformation. The alternated solution is performed iteratively, until a convergence criterion is met. Convergence criteria are e.g. no further changes in the increments are obtained for all DOF or the norm obtained for a post-processed monolithic residuum is below the tolerance limit.
The application of the staggered solution scheme to phase-field simulations is required for static crack propagation, as the monolithic solution cannot obtain convergence in the case of brutal phase-field evolution, see . Instead, a staggered solution is proposed in . Originally, a single solution for each increment of the mechanical field and the phase-field is motivated by very small load steps during the crack evolution, which is replaced by the staggered iterations in later phase-field approaches. In transient simulations, the inertia, which is not affected by the phase-field, acts as a stabilization to the solution and the monolithic solution scheme is applicable, see e.g. Borden (2012) and Steinke et al. (2016).

Crack orientation vector
A major ingredient of the directional split is a realistic definition of the crack orientation vector r. The explicit definition of the crack orientation forms further the basis of the generalization of the directional split in the framework of representative crack elements, see e.g. Storm et al. 2020, Fig. 5 Representation of the crack orientation r at the crack tip by the gradient type definition r = ∇ p/|∇ p| 2021; Yin et al. 2021). The definition r = ∇ p/|∇ p| proposed in Strobl and Seelig (2015) seems an intuitive approach. It results in a realistic orientation in the special case of the through crack discussed by Strobl and Seelig. However, this definition fails both at the crack tip as well as within fully degraded elements. At the crack tip, the gradient of the phase-field does not represent the actual orientation of the crack, see Fig. 5. While this is the critical region, where crack propagation is expected, the phase-field evolution is affected by a wrong decomposition of the stresses. Furthermore, the stress free boundary of the crack surface requires an approximation by at least one row of connected fully degraded elements, see Steinke et al. (2016). However, in a fully degraded element, i.e. all nodes having p = 1, the gradient of the phase-field is not specified properly. While this region is critical for the proper approximation of the crack kinematics, again the gradient definition fails to deliver a realistic approximation of the crack.
In the following, two approaches to define the orientation of a phase-field crack are proposed. The focus during the development is on a realistic approximation of the actual crack orientation of the represented discrete crack, especially during crack propagation. At first, a definition based on an evaluation of the principal stresses is proposed, that resembles closely to the principal stress failure criterion of early fracture mechanics theory. The second approach is based on an analytical evaluation of the crack driving strain energy density at the material point level.

Principal direction
The evaluation of the principal direction of stresses or strains is well established in early approaches to model fracture. Furthermore, it sustains a significant influence on damage and failure simulation, see Gross and Seelig (2011). Indeed, the formation of a crack surface per-pendicular to the largest tensile strain or stress is an intuitively realistic approach based on everyday life's observation.
In an isotropic linear elastic material, the principal directions of strains and stresses are identical. A similar conclusion is obtained for the strains and the effective stress tensor for special cases, e.g. the spectral and volumetric deviatoric split, see Steinke (2021) for a detailed investigation of the spectral split on this issue. However, an analogous derivation is hardly possible for the directional split. The reason is the composition of the stress tensors in Eqs. (6) and (7), which depend on the crack orientation vector r, as well as the dependence on the specific value of the degradation function. In consequence, analytical conclusions on the principal directions of the effective stress tensor of the directional split, i.e. whether or not they are identical to the principal directions of the strain tensor, are not available. Rather, the principal directions of the strains are evaluated and employed.
A general definition of the crack orientation vector by a principal (strain) direction criterion can be provided by where n 1 is the direction of the largest principal strain, i.e.ε 1 >ε 2 holds. In combination with the irreversibility approach, the change of the crack orientation is formulated in an explicit scheme by This enables the change of the crack orientation during the evolution of the crack, while keeping the orientation of the evolved crack in unloading and reloading situations.

Energetic crack orientation in two dimensions
The fracture toughness G c is a measure for the material inherent energetic barrier against crack evolution.
In combination with the crack surface density γ l , the resistance against crack propagation is quantitatively specified at the material point level. In the phase-field method, the crack driving component of the local strain energy density is the energetic counterpart to that resistance. Crack evolution is a result on the global level, that is based on the local formulation and evaluation of driving force and resistance by Eq. (17). The specific amount of crack driving strain energy density depends on the state of strain. Furthermore, it strongly depends on the orientation of the crack surface. In the previous section, an intuitive alignment of the crack orientation vector, i.e. along the direction of the largest principal strain, is postulated. It can be shown, that in the majority of possible strain states, this approach yields the maximum driving force for isotropic elastic material due to maximization of the crack driving strain energy density with respect to the crack orientation vector. However, certain states of strain exhibit the maximum in the crack driving strain energy density for a crack orientation, that is at a specific angle to the principle direction. Also recent experimental observations in Rozen-Levy et al. (2020) support the hypothesis that the principal of maximum energy dissipation also applies to the propagation of cracks.
In this section, an analytical approach to obtain the maximum crack driving component of the strain energy density is proposed and evaluated for a linear elastic material in a 2D setup. At the basis of known eigenvectors n 1 and n 2 of the 2D strain tensor ε, a principal direction coordinate system (PCS) is established. The representation of both eigenvectors in the RCS and the PCS reads and where the RCS and the PCS are indicated by subscripts x and n, respectively. The representation of the base Fig. 6 Relation between the reference coordinate system (RCS), the principal direction coordinate system (PCS) and the crack coordinate system (CCS) vectors of the CCS in terms of the PCS reads and Here, an angle α between the principal direction n 1 and the crack orientation vector r is considered according to Fig. 6. In analogy to the base vectors, the strain tensor is rewritten as A relation between the principal strainsε 1 andε 2 is introduced based on the local strain magnitudeε bŷ By definition, holds, i.e.ε 1 is the major principal strain, n 1 is the major principal strain direction and the spatial orientation of r is defined in an unique manner, see Fig. 6. The linear elastic ground stress tensor σ 0 of a linear elastic material with strain energy density provided in Eq. (4) is given in PCS representation by The tensor bases specified in Eqs. (6) and (7) are given with respect to the representation of r and s in the PCS by and The stress tensor components, relevant for the crack driving strain energy density, read σ rr = σ 0 : (r ⊗ r) = λ ε 1 +ε 2 + 2μ ε 1 cos 2 (α) +ε 2 sin 2 (α) and The crack driving stress tensor σ + contains components related to tension on the crack surface and contributions related to shear on the crack surface in an additive decomposition, that reads where k = λ/(λ + 2μ). For brevity of notation, the trigonometric functions are abbreviated by s = sin(α) as well as c = cos(α).
Beside the comprehensive study of all possible states of strain, the concept of a modified F-criterion, see e.g. Shen and Stephansson (1994), is considered as well, i.e. a mode dependent fracture toughness is analyzed. Therefore, also the crack driving energetic counterpart has to be decomposed into components associated with Mode I and Mode II. An application in the framework of the phase-field method is published in Zhang et al. (2017), where a volumetric deviatoric split is employed. Thus, Mode I is associated with volumetric expansion and Mode II is related to a deviatoric deformation. While this results in a realistic approximation of wing cracks and secondary cracks in pre-notched sandstone specimens, the connection of the continuum measures of volumetric and deviatoric strains to the representation of Mode I and Mode II deformation is ambiguous. In contrast, an association of the components in a directional decomposition to the modes of crack opening is much more straight forward. Here, tensile components in r ⊗ r direction can be associated with Mode I and shear components in r ⊗ s and s ⊗ r direction are related to Mode II. Thus, ψ + is decomposed into the Mode I crack driving strain energy density and a Mode II crack driving strain energy density A reference fracture toughnessĜ c is introduced in order to specify the relations Furthermore, the range of relations n investigated is restricted to materials characterized by Then, in analogy to the modified F-criterion, mode dependent driving forces are introduced via and The total driving force for phase-field evolution reads For short notation, two additional material parameters are introduced as and which are bounded due to the range of Poisson's ratio 0 ≤ ν < 0.5 by 0 ≤ q and 0 ≤ k ≤ 1, respectively, see Fig. 7 for a visualization of the valid ranges for q and k in logarithmic scale. The individual crack driving strain energy density components are weighted by their corresponding fracture toughnesses. Therefore, the phase-field evolution equation Eq. (17) is rewritten as The mode dependent driving force F for the phasefield evolution equation according to Eq. (42) is a function of the angle α. The application of the modified F-criterion within the framework of the phase-field method is based on the postulate "phase-field evolution due to local maximum of driving force at the material point level governs the global energy dissipation realistically". Thus, an analytical discussion of the function is necessary to find its maximum value and to obtain the proper driving force for phase-field evolution at each material point. To this end, the condition for an extremum is analyzed in order to determine all possible α 0 . Then, the second derivative condition is used in order to evaluate, whether F(α 0 ) is a maximum. The discussion of equations, that contain trigonometric functions, results in an infinite number of possible solutions. However, in terms of the orientation of the crack orientation vector r, the range for α 0 can be restricted to −90 • ≤ α 0 ≤ 90 • . In other words, α 0 = 0 • represents the case, where the evolving crack is aligned to the major principal strain and α 0 = ±90 • represents the case, where the evolving crack is aligned to the minor principal strain. Detailed information on the derivations are provided in Steinke (2021). In the following, the results of the discussion are summarized for the main aspects.
The investigation of small strain linear elasticity with isotropic material behavior and fracture toughnesses results in the noteworthy aspect, that the energetic level of the driving force for phase-field evolution according to Eq. (42) at an angle α 0 is exactly the same as for its negative counterpart, i.e.
(49) In this case, an energetical branching point is obtained, i.e. the kinking angle of the crack is not well defined due to the continuum mechanical approximation of the discrete kinematics of the crack. At the level of the crack driving stress tensor, the sign of the shear components in its PCS representation is the opposite, i.e.
This is a direct consequence of a basic assumption of the directional split, i.e. the decomposition of the stress tensor with respect to the orientation of a single crack. In fact, having identified a point of branching energetically, also the kinematics of a branched configuration should be considered. First ideas to develop a suitable approach are based on the framework of the representative crack element introduced in Storm et al. (2020). However, a final solution is not available yet. Instead, the intrinsic capability of the phase-field method to obtain crack branching on a global level is exploited and a unique local orientation of the crack is chosen. In general, it should be noted, that the update of the crack orientation vector is based on an explicit scheme in this context. Furthermore, an interesting yet undiscussed aspect is the actual meaning of a change of crack orientation during its evolution. Due to the continuous representation of the phase-field, intermediate states between sound ( p = 0) and broken ( p = 1) are possible at the material point level. The local evolution of the phase-field is a contribution to the global crack approximation l . However, it is always linked to a local crack orientation vector r.
On the one hand, this results in the possibility, that the orientation of the crack, which is a uniquely defined property of every point on the surface of a discrete crack, can vary over the profile of its continuous representation. This could be interpreted as an actual violation of a proper approximation scheme. However, if two or three dimensions are considered, it is hardly possible to associate every point in the phase-field profile to the discrete point on the crack surface it represents. Furthermore, a major role for the realistic approximation of the crack kinematics is played by the fully degraded state at p = 1 and the influence of the transition zone is of only minor degree.
On the other hand, the local orientation may change during the local evolution of the phase-field. The fundamental question is, whether the energy dissipated for the local evolution in a specific direction can be transferred to the updated direction without modification. From the viewpoint of a discrete crack, the answer is definitely no, as a different crack orientation at the same location more likely represents another crack at another orientation, e.g. a point of crack branching. Therefore, in the case of a discrete crack, tearing apart the material in a different direction requires additional energy because a different atomic bond has to be broken. It could be concluded that the energy dissipated for the initial direction in the first place cannot contribute to the second crack with a different orientation.
Nevertheless, as already stated, kinematics for a branched crack at the material point is not available yet and only single cracks are considered. Furthermore, the number of possible crack orientations at a single material point is infinite and the computation and storage of individual driving forces for each orientation is not feasible.
Therefore, the level of energy dissipated in the initial direction is considered as the starting point for the evolution in the new direction. If an energetic branching point is obtained, e.g. in the sense of Eq. (49), then the new crack orientation is chosen such that the previous orientation has a minimum angle to the new crack orientation. The same holds, if multiple energetically equivalent orientations are obtained. Furthermore, the change of crack orientation is allowed for the transition zone only, i.e. as long as the phase-field is below a critical value p c . Together with the modified history variable approach proposed in Eq. (46), a realistic approximation of the crack kinematics is obtained at the global level with this approach.
Instead of a general analysis of Eq. (42) for an arbitrary state of strain, the strain state is categorized into characteristics based on the magnitudes of the principal values and the analysis is realized in terms of a case study. A positive local strain magnitudeε > 0 specifies a strain state, where the major principal strainε 1 is a tensile strain. Then, based on the value of m, the strain state can be verbally described as Principal strains for a local strain magnitude ofε = 0 cannot be determined and the relation according to Eq. (25) cannot be applied. Due toε =ε 1 ≥ε 2 , only uniaxial compression in n 2 direction is possible in this case.
A negative local strain magnitudeε < 0 specifies a strain state, where the major principal strainε 1 is a compressive strain. Then, based on the value of m, the following strain states are possible -ε 2 =ε 1 → m = 1 -volumetric compression and -ε 2 <ε 1 → m > 1 -biaxial compression.
Volumetric tension A state of volumetric tension yields a constant driving force of for phase-field crack evolution, which is independent on the angle α, see Fig. 8.
Biaxial tension A state of biaxial tension with dominant tension in n 1 direction yields a maximum driving force of for phase-field crack evolution at an angle of α = 0 • , see Fig. 9. Uniaxial tension A state of uniaxial tension in n 1 direction yields the maximum finite amount of driving force of for phase-field evolution at an angle of α = 0 • , see Fig. 10 with n = 1.

Biaxial tension compression with dominant tension
A state of biaxial tension compression with dominant tension in n 1 direction yields a maximum driving force for phase-field evolution of In addition, a maximum driving force for phase-field evolution of is obtained at an angle of in the case of m < m dt ∧ n < n dt . The visualization of the single and double maxima, that depend on the combination of m and n, are given in Figs. 11 and 12, respectively.
Pure shear A state of pure shear yields a maximum driving force for phase-field evolution of at α = 0 • ∧ n > n dt , see Fig. 13. For n = n dt , a constant driving force identical to Eq. (59) is obtained for the whole range of 0 • ≤ α ≤ 45 • , see Fig. 14. In the case of n < n dt , which is visualized in Fig. 15, the maximum driving force is obtained at α = 45 • with Biaxial tension compression with dominant compression A state of biaxial tension compression with dominant compression in n 2 direction yields a maximum driving force for phase-field evolution of An energetic branching point is obtained at α = ±45 • for m < m dc . In this case, the maximum driving force reads An energetic triple branching point is obtained for m = m dc , i.e. the driving forces are identical at α = 0 • and α = ±45 • and can be computed according to Eqs. (61) or (63), see Fig. 16 for a visualization.
Uniaxial compression A state of uniaxial compression in n 2 direction yields the maximum finite amount of driving force Volumetric compression A state of volumetric compression yields no driving force for phase-field crack evolution.
Biaxial Compression A state of biaxial compression with dominant compression in n 2 direction yields the maximum finite amount of driving force for phase-field evolution at an angle of α = 45 • , see Fig. 18.

Crack surface approximation
The approximation of the discrete crack surface by its regularized approximation l is governed by Eq.
(2). In the following, characteristic 2D crack patterns are studied with respect to the quality of this approximation. The focus is on the relation between the regularization length l and the characteristic size of the finite element discretization h, that is necessary for an adequate approximation.

Single through crack
A rectangular patch of dimensions 1.0 m × 1.0 m with a discrete, horizontal through crack is investigated, see Fig. 19. The crack's approximation by l is based on a finite element solution with phase-field boundary conditions at nodes located at the discrete crack. In principle, the phase-field approximation of can be represented in two ways, i.e. as a single row of nodes with p = 1 or as a row of elements, i.e. two rows of nodes with p = 1, see e.g. Steinke (2021) and Steinke et al. (2016) for a transient analysis of the realistic approximation of crack kinematics as well as the discussion in Sect. 4.5 for the static case. Furthermore, finite element discretizations with triangular and quadrilateral elements with characteristic size h = 5 mm are studied. The guidelines of the crack approximation based on a row of nodes and details of the discretization with triangles (SST) and quadrilaterals (SSQ) are shown in Fig. 20a. The guidelines for the crack approximation based on a row of elements and the details of the discretization with triangles (SDT) and quadrilaterals (SDQ) are given in Fig. 20b.
The phase-field distributions obtained by the finite element simulation for SST and SDT are visualized in Fig. 21 for different regularization length parameters l. The differences between the representation of the crack by a row of nodes or a row of elements is restricted to the A more detailed investigation is possible based on a plot of the phase-field profile in y-direction, provided in Figs. 22 and 23 for SST and SDT, respectively. First of all, it should be noted, that negative phase-field values are observed at the nodes directly adjacent to the node with the phase-field boundary p = 1 with l = h/5. This is due to the fact, that the discretization with 2D finite elements, which are based on linear shape functions, is not fine enough to resolve the phase-field gradient properly. In consequence, the nodes along the profile exhibit oscillating positive and negative phase-field values with fast decreasing amplitude. Furthermore, due to the AT-2 functional studied here, the profile of the phase-field is rather wide and does not fit into the patch of height 1 m for l > 20h(= 10 cm). The close up view on the profile for SDT in Fig. 23 also visualizes the plateau of the phase-field values at p = 1 for the row of elements. However, this does not influence the gradient of the profile, which is entirely governed by the value of the regularization length l.
The length of the discrete crack in this example is given by = 1.0 m. The amount of regularized crack surface l is computed based on Eq. (2) for an integration over the whole finite element domain. The result of Fig. 23 Phase-field profile for SDT with respect to the regularization length l the study of the regularization length l = h/5 . . . 40h is provided in Fig. 24. It should be noted, that a good approximation of / l ≈ 1.0 is obtained both for SST and SSQ at l > 2h. On the one hand, this indicates, that the use of triangular or quadrilateral elements is equivalent in terms of the regularized crack approximation. On the other hand, this result confirms the conclusions obtained in Hofacker (2013), i.e. the lower limit for the regularization length l > 2h. In contrast, the results obtained for the crack approximation based on a row of elements, i.e. SDT and SDQ, are converging much slower to the upper limit. Furthermore, the upper limit in this case cannot be / l = 1.0 in principle, as the plateau contribution is a fundamental violation of the theoretical framework of crack approximation by the phase-field. Nevertheless, the impact of the plateau on the value of l is decreasing with increasing l. In detail, the error decreases from / l = 79.5% at l = 2h to / l = 97.5% at l = 20h. A further increase of the regularization length is not possible in this study due to the excess width of the phase-field profile with respect to the size of the specimen, see Fig. 22.

Crack tip approximation
A rectangular patch of dimensions 1.0 m × 1.0 m with a discrete, horizontal crack of length = 0.5 m is investigated, see Fig. 25. Two alignments of the crack are studied, i.e. an alignment to the left edge of the specimen and an alignment at the center of the specimen, in order to investigate the influence of one and two crack tip approximations, respectively. The phasefield approximation of the discrete crack is obtained by phase-field boundary conditions p = 1 at the loca- tion of the crack. The left aligned crack is studied for a regularized crack approximated by a row of nodes with triangular elements (HLST) and quadrilateral elements (HLSQ) as well as for an approximation by a row of fully degraded triangular elements (HLDT) and quadrilateral elements (HLDQ). Accordingly, these cases are studied for the centered alignment of the crack as well, that are referred to by the abbreviations HCST, HCSQ, HCDT and HCDQ, respectively. The finite element discretization of the single through crack example discussed in Sect. 4.2 is used, see Fig. 20. The phase-field distribution for HCST and HLST is visualized in Fig. 26. It should be noted, that for the center aligned crack, the width available for the phasefield profile of the crack tip is limited and the study should be restricted to l ≤ 10h in this case.
Given the amount of discrete crack surface = 0.5 m to be approximated in both cases of this example, the quality of the crack surface approximation can be analyzed with respect to the regularization length parameter l, see Figs. 27 and 28. First of all, the discretization with triangular and quadrilateral elements is identical in all cases. Thus, the choice of the element type is of minor relevance for the approximation of the crack surface again. Furthermore, a limit value of / l cannot be observed. In contrast, the influence of the crack tip approximation increases with increasing regularization length. A conclusion of the example discussed in Sect. 4.2 is the fact, that the discrete length of the crack is already represented by the phasefield profile in vertical direction with a very good accuracy. In contrast, the crack tip does not have a finite length. However, it is approximated by a phase-field Fig. 29 Geometry of the branched crack example and guidelines of the finite element discretization for different branching angles α profile with semicircle shape, see Fig. 26, that increases in radius with increasing regularization length. Therefore, increasing the regularization length parameter increases the error due to the semicircle crack tip approximation. Interestingly, the approximation based on a single row of nodes yields the most accurate result at l = 2h with / l = 96.7% and / l = 98% for HCST and HLST, respectively. In contrast, if the crack is approximated by a row of fully degraded elements, the most accurate results are obtained at l = 8h, where / l = 86.2% for HCDT, and at l = 10h, where / l = 90.1% for HLDT. Nevertheless, it should be noted, that the profile of the crack tip is usually small compared to an evolved crack and the impact of its semicircle profile approximation decreases with increasing total length of the crack.

Branched crack approximation
A branched crack configuration is examined in the last example on the quality of the crack surface approximation by the phase-field method. A rectangular patch of dimensions 0.5 m × 1.0 m is joined to a semicircle patch with radius 0.5 m, see Fig. 29. The discrete crack to be approximated consists of a horizontal straight component of length 0.5 m and two branches of length 0.5 m each, i.e. the total length of the discrete crack is = 1.5 m. The angle between both branches is specified by 2α and studied for α = 5 . . . 85 • in steps of 5 • . The semicircle shape is chosen in order to minimize edge effects on the profile of the crack branches. As already observed in the previous examples, the profile perpendicular to the crack surface is essential for a correct approximation. Hence, the semicircle shape can decrease the impact of the edges on the profile for small regularization length parameters by ensuring a 90 • angle between the discrete crack and the specimen edges. However, especially for large l, the profile is trimmed, which results in an underestimation of the regularized crack surface in general. In this example, only the crack approximation by a single row of nodes and a discretization with triangular finite elements is investigated. The guidelines of the meshing are the outer boundary of the specimen as well as the location of the discrete crack. Furthermore, exemplary discretizations of the branching point are visualized in Fig. 29 for α = 10 • , α = 45 • and α = 80 • .
As already discussed, the phase-field profile perpendicular to the crack is crucial for a correct representation of the discrete crack by its regularized approximation l . Therefore, the optimal case for a phase-field crack approximation is given by the initial example of the straight through crack in Sect. 4.2. On the one hand, the through crack does not have any crack tip that yields a significant error for the regularized crack surface l , see Sect. 4.3. On the other hand, the phase-field profile is not disturbed, e.g. by edges or neighboring cracks, and can form a flawless transition from the broken to the sound region of the specimen, which results in a very accurate approximation of the discrete crack surface. In contrast, if a branched configuration is investigated, Fig. 31 Quality of the crack surface approximation of a half crack at center alignment with respect to the ratio between element size h and regularization length l the vicinity of the branching point is a location of multiple supports of phase-field profiles overlapping and influencing each other. Starting at the branching point along the crack path in direction of the edges, the influence of the neighboring cracks on the phase-field profile decreases. However, the decreased influence strongly depends on the angle between the branches, see Fig. 30 for a visualization of phase-field profiles at α = 5 • , α = 30 • and α = 50 • for exemplary regularization lengths. While the individual branches can be identified clearly for α ≥ 30 • even for large values of l, the branches are merged and represented by a single fat phase-field profile for small branching angles. The overlapping results in an underestimation of the crack surface in principle, see Fig. 31. The underestimation is strongly influenced by the branching angle. Furthermore, a pronounced change is observed for α ≥ 30 • , i.e. the influence of the branching angle is decreasing significantly for this range. This also has an impact on the choice of l/ h for the best approximation of . In the range of α ≥ 30 • , an exact approximation with / l = 100% is obtained for l = 2h. For smaller branching angles α < 30 • , the optimal choice of the regularization length is decreasing as well until l = h for α = 5 • , which yields / l = 99.5%.

Kinematics of ideal plane cracks without friction
The crack kinematics for the special case of an ideal plane crack surface without frictional effects is discussed in Sect. 2.2 and encompasses the three types of structural deformation crack closure, crack opening and crack face sliding. A realistic approximation of cracks by the phase-field method must be able to reproduce these features in a numerical setting. The quality of the approximation is evaluated based on the numerically obtained reaction forced at the structural level as well as the local values of the stress tensor. In the following, benchmark examples with initial phasefield cracks at static and dynamic loading are presented and analyzed in order to evaluate different phase-field splits with respect to a realistic post fracture behavior. The focus of the last two examples is on the specific challenges arising from the implementation of a directional split in the framework of FEM.

Static investigation of crack kinematics
In the first example, a 2D single through crack is investigated, that is similar in geometry to the example studied in Sect. 4.1, i.e. a rectangular patch of dimensions 1.0 m×1.0 m with a discrete, horizontal through crack halfway up in vertical direction is examined, see Fig. 19. The through crack is approximated by an initial phase-field crack represented by phase-field boundary conditions p = 1 at a row of nodes or all nodes of a row of elements, see the discretizations in Fig. 20 for SST and SDT, respectively, where the nodes with a phasefield boundary condition are indicated by a horizontal thick line. In addition to the nodal phase-field boundaries at the location of the initial crack, displacement boundariesû x (t) =û · f x (t) andû y (t) =û · f y (t) are applied at the top edge of the specimen. Furthermore, the displacements at the lower edge are fixed in xand y-direction. The magnitude of the displacement iŝ u = 1 mm and the pseudo time dependent scaling functions are visualized in Fig. 32. Initially, a horizontal displacement is applied to investigate crack face sliding.  Dir 29.2 × 10 −3 31.2 × 10 9 73.1 × 10 −3 SST Vol dev 25.2 × 10 9 30.6 × 10 9 6.9 × 10 9 Spec 25.2 × 10 9 31.4 × 10 9 7.9 × 10 9 Dir 25.2 × 10 9 31.4 × 10 9 7.1 × 10 9 Linear elastic 31.4 × 10 9 31.4 × 10 9 8.2 × 10 9 Without crack Discrete crack 0 31.4 × 10 9 0 Then, the subsequent displacement in vertical direction results in crack opening followed by crack closing. The structural response is evaluated by a homogenized stiffness that is computed based on the mean value of the displacement U and the sum of the numerically obtained reaction forces F at the upper edge in the specific direction. The material parameters λ = 7.15 GPa and μ = 12.71 GPa are used and lead to K s = 8.184 GN/m and K c/o = 31.432 GN/m, for shear deformation and a compressive/tensile (closing/opening) deformation in vertical direction, respectively. It should be noted, that these values are obtained with the SDT mesh, but without a phase-field crack, i.e. they describe the linear elastic behavior of a bulk rectangular patch with displacement boundaries as mentioned above. These homogenized stiffnesses provide a measure for the structural response of the rectangular patch with respect to the material parameters and serve as a comparative value for the evaluation of the phase-field approximation. According to the considerations on crack kinematics visualized in Fig. 3  should be noted, that the opening and sliding is possible without reaction forces for the volumetric deviatoric split in SDT only, i.e. a realistic post fracture behavior is obtained, if the phase-field crack is approximated by a row of elements with each node having p = 1. The explanation of this phenomenon is obtained by a detailed investigation of the finite element implementation of the phase-field method for a crack opening, see Fig. 35. Let A and B the two parts of the rectangular patch, that are separated by the horizontal through crack. In the FEM discretization, these parts are connected via finite elements with modified stiffness according to the local phase-field values, that may be simplified as stiffness k in vertical direction, which is actually a combination of components of the elemental stiffness matrices. For a crack opening deformation, the split affects the correct identification of the strain and stress components associated with the structural opening direction. However, it is the degradation function g( p), that defines the amount of degradation of this specific stiffness. As a fundamental aspect of finite element implementation, the elemental stiffness is a result of the numerical integration of the element volume, i.e. the summation of contributions at the Gauss points of the element according to basic finite element theory. If the crack is approximated by a single row of nodes with p = 1, then the Gauss point value of the phase-field for the elements connected to those nodes is always lower than one. Therefore, the element stiffness is only partially degraded by g( p < 1) > 0 and a realistic behavior cannot be simulated. While the explanation and visualization of this feature in Fig. 35 is provided for the crack opening deformation, an analogous argumentation applies to the stiffness component for shear deformation. Furthermore, it should be noted, that this is a fundamental aspect of the phase-field implementation. Although the Gauss point value of the phase-field is more close to p = 1 for an increased regularization length l, a Gauss point value of p = 1 can only be obtained, if all nodes associated with an element exhibit p = 1. The homogenized stiffness for crack closure is only 93.3% of the expected value for the volumetric deviatoric split. The reason is visualized in Fig. 33d and is directly related to the volumetric deviatoric decomposition of the strains. Lets assume the local strain for crack closure is an uniaxial compressive strain in vertical direction of magnitudeε yy , which would be the case if the horizontal displacement of the upper and lower edge of the specimen is free. According to the considerations on crack kinematics visualized in Fig. 3, no degradation should be applied in this case. However, the volumetric deviatoric decomposition of a uniaxial strain in vertical direction reads The volumetric component is compressive, i.e. it is not associated with degradation in the volumetric deviatoric split. In contrast, degradation is always applied to components related to the deviatoric part of the strain. Therefore, a finite component of the initial uniaxial strain is associated with degradation and a wrong homogenized stiffness for crack closure is obtained. While this is a fundamental feature of the volumetric deviatoric split and applies to every element, an additional aspect is observed at the edge of the specimen, see i.e. the amount of deviatoric strain is increased. This results in increased degradation of the horizontal stiffness, which eventually yields an increase in the horizontal strain up to the point ofε xx =ε yy , i.e. the non-degraded volumetric components vanish totally. While this is the reason for the defective horizontal displacement shown in Fig. 33d, it also results in nonconvergent results for large amounts of crack closure deformations for the volumetric deviatoric split. The structural deformations obtained by the spectral split are visualized with a deformation scaling factor of 50 in Figs. 36 for the SDT discretization. The results with SST discretization are omitted, as they suffer from the same fundamental deficiency discussed above for the volumetric deviatoric split. The spectral split exhibits realistic post fracture behavior in the case of crack closure and crack opening, where the homogenized stiffnesses are identical to the expected ones within the tolerances of a numerical approximation, see Tab. 1. However, significant reaction forces are obtained in the case of crack face sliding, that result in a homogenized stiffness for shear deformation that amounts to 91% of linear elastic behavior without a crack. Again, the reason to the unrealistic post fracture behavior is given by the decomposition of the strain tensor. Let the structural shear deformation be simplified by a strain tensor, that consists exclusively of shear components of magnitudeε xy . Its decomposition according to the spectral split reads i.e. it is considered to be a combination of a tensile principal strain of magnitudeε 1 and a compressive principal strain of magnitudeε 2 . Furthermore, pure shear results inε 1 =ε 2 . In the spectral split, tensile principal strains are associated with degradation and compressive principal strains are related to persistent components of stress and strain energy density. Therefore, the decomposition of the shear strain, that should be totally degraded for a realistic post fracture behavior, results in a significant amount of non-degraded contributions to the elemental stiffness and the global structural response. The structural deformations obtained with the directional split are visualized with a deformation scaling factor of 50 in Fig. 37 for the SDT discretization. Due to the explicit design of the split to approximate a realistic post fracture behavior, both the structural deformations are in good agreement with the expected behavior and the homogenized stiffnesses provide values within the tolerances of a numerical approximation, see Tab. 1. It should be noted, that the fundamental deficiency regarding a proper crack approximation within an SST discretization is observed for the directional split as well, as this is an aspect of the FEM implementation rather than the phase-field split. In consequence, a row of fully degraded elements is a basic requirement for the approximation of an initial phase-field crack.

Type and orientation of finite elements
The realistic post fracture behavior of a 2D rectangular patch with an initial phase-field crack and a directional split is investigated with a special focus on the type and the orientation of the finite elements applied for the discretization. An initial crack with orientation r = [0, 1] T divides the specimen into two separate bodies A and B. Three types of discretization are analyzed, see Fig. 38. Two discretizations are based on quadrilateral elements with quadratic shape. In ETQ90, the element edges are aligned with the surface of the initial crack. In ETQ45, the element edges are at an angle of 45 • to the surface of the initial crack. It should be noted, that a continuous sequence of fully degraded elements along the location of the initial crack requires a zig-zag pattern of totally degraded elements, which results in a significantly larger plateau of p = 1 in the phase-field profile. The discretization ETT is based on triangular elements. According to the results of previous examples in Sect. 4.5, a row of fully degraded elements is considered to obtain a realistic post fracture behavior. The elements investigated are indicated by thick black lines in Fig. 38, where some nodes belong to body A and some nodes to body B. The characteristic element size is h = 5 mm and material parameter λ = 8.89 GPa and μ = 13.33 GPa are assumed.
The numerical integration of field quantities in finite elements is based on Gaussean quadrature formulas and is discussed in more detail e.g. in Zienkiewicz (1977). A 2D quadrilateral element with linear shape functions requires four Gauss points at the local coordinates ξ ≈ ±0.577 and η ≈ ±0.577, see Figs. 39 and 40. The summation over the Gauss point contributions is the basis for the computation of elemental quantities in the FEM framework like the residual vector and the stiffness matrix, where the definition of the crack driv- ing component of the stress tensor σ + is crucial. The elements under consideration in this example are indicated by a thick line in Fig. 38. They are analyzed as representatives of all elements that link the two bodies A and B separated by the initial crack. The residual at the nodes of these elements describes the local reaction forces against a specific structural deformation at the position of the specific element. The reaction force at the structural level is the result of all the local reaction forces obtained for the fully degraded elements. A realistic post fracture behavior at crack opening and crack face sliding deformation is characterized by no reaction forces. Therefore, the contributions of the fully degraded elements to the nodal residual in these cases must be zero. This implies a similar constraint on all Gauss points, as the residual vector is a result of the summation over all Gauss point contributions of the element.
The orientation of quadrilateral elements in ETQ90 is characterized by the element edges being parallel to both the plane of the initial crack as well as the crack orientation. Therefore, the crack opening and crack face sliding deformation at the global level are represented by a vertical and horizontal deformation of the upper edge of the element, respectively, see Fig. 39. In consequence, a uniform state of strain is obtained, that is identical for all Gauss points and reads for a vertical displacement u y = 1 mm and a horizontal displacement u x = 1 mm, respectively. The directional decompositions of the associated ground stress tensors read In both cases, the total stress is associated with the crack driving component of the stress tensor σ + . Thus, at full degradation, i.e. g( p) = 0, the reaction forces at the element level are zero and a realistic structural behavior is obtained, see the example in Sect. 4.5.1. Due to the discretization in ETQ45, nodes 1, 2 and 4 belong to the body A and node 3 belongs to body B. Therefore, crack opening and crack face sliding deformations result in a displacement of node 3 only, see at the Gauss points G 1 to G 4 , respectively, which results in the directional decompositions of the associated ground stress tensors as For Gauss points 2 and 4, non-zero persistent stress tensor contributions σ − are obtained, which causes non-zero components in the residual vector of the element. In consequence, a global crack face sliding deformation, that is actually without resistance for an ideal crack, results in significant reaction forces in the finite element approximation. It should be noted, that this behavior can be observed in every discretization, that deviates from the perfect alignment of element edges and crack orientation discussed as ETQ90. It is a direct consequence of the bi-linear approximation of the field variables inside the quadrilateral elements, that is already pointed out in Sect. 2.4. Linear shape functions within finite triangular elements result in a constant field of strain all over the domain of the finite element. Thus, the numerical integration of the function with respect to the strain field can respectively. As they are similar to the strains obtained in ETQ90, again, the realistic post fracture behavior is obtained. Furthermore, the orientation of the triangular element does not affect the strain field, as the strain field is constant. Therefore, in contrast to discretizations with quadrilateral elements, a realistic post fracture behavior can be obtained regardless the discretization. It should be noted, that this issue is present for 3D simulations as well, where brick element discretizations exhibit a tri-linear field of strain and tetrahedral elements provide a constant field of strain. Nevertheless, a constant approximation of strain results in a poor accuracy of the numerical results and a comparable accuracy to the bi-and tri-linear approximation of quadrilaterals and bricks, respectively, can be obtained only by an increase of the discretization density. Fortunately, the phase-field crack approximation requires a very high discretization density as well and the application of triangular and tetrahedral elements is a suitable approach to obtain a realistic post fracture behavior in principle.

Distortion of the crack orientation vector
A rectangular patch of dimensions 1.0 m × 1.0 m with an initial horizontal through crack , dividing the specimen into bodies A and B, is investigated. The initial crack is approximated by a row of fully degraded elements and the whole specimen is discretized by a mapped meshing with triangular elements with characteristic size h ≈ 1 mm, see Fig. 42. The crack ori- During the evolution of phase-field cracks, a deviation of the crack orientation between neighboring elements is essential to approximate crack kinking or branching. Furthermore, the deviation is observed as a result of the stress field around the crack tip, see e.g. Steinke and Kaliske (2018).
The phase-field profile is computed once at the beginning of each simulation and the subsequent computations are restricted to the mechanical DOF with fixed crack orientation vectors for each element. A crack opening deformation is obtained by a displacement boundary of the upper edge in vertical direction with a magnitude of u y = 1 cm. Analogously, the crack face sliding deformation is studied for a displacement boundary of the upper edge in horizontal direction with a magnitude of u x = 1 cm. It should be noted, that both cases result in a rigid body movement of body B, regardless of the deviation angle φ, i.e. the impact on A summary of the results with respect to the deviation angle φ is visualized in Fig. 44. It is recapitulated from Sect. 4.5.1, that the homogenized stiffnesses of a linear elastic specimen without a crack are K o = 31.43 × 10 9 N/m and K s = 8.19 × 10 9 N/m for crack opening and crack face sliding deformation, respectively. These values are employed for an evaluation of the numerical results with deviated crack orientations. Furthermore, it should be considered, that the actual discrete crack kinematics discussed in Sect. 2.2 require zero stiffness in both cases. At the structural level, the homogenized stiffnesses obtained due to a deviation of the crack orientation of a single element are small for a crack opening deformation with a maximum of 2.67% of the linear elastic value at φ ≈ 70 • , which is a very large angle of deviation anyway. In contrast, the homogenized stiffnesses for a crack face sliding deformation amount to 13.43% of the linear elastic value at φ ≈ 40 • .
In the following, the stress fields for three representative angles of deviation φ are evaluated. Again, the linear elastic response of the specimen without a crack are considered for the evaluation of the stress field magnitudes. Due to the displacement boundary conditions in horizontal and vertical direction at the upper and lower edge of the specimen, the linear elastic results for crack opening and crack face sliding deformation with magnitude 1 cm each, respectively. A deviation angle of φ = 5 • results in homogenized stiffnesses K o = 0.47 × 10 6 N/m and K s = 62.51×10 6 N/m for crack opening and crack face sliding deformation, respectively. Both values are small in comparison to the stiffnesses obtained for a linear elastic specimen without a crack, i.e. 0.0015% and 0.76% for opening and sliding, respectively. Nevertheless, the peak stress magnitude forσ xx amounts to more than three times the linear elastic value for a crack opening deformation, see Fig. 45a. Interestingly, a tensile peak is observed at the left upper edge of the element with the deviation of the crack orientation vector and a peak compressive stress of identical magnitude is observed at the right lower edge below the element with the deviation of the crack orientation vector. At a first guess, this looks like an inclusion of stresses, where the positive and the negative magnitude neutralize each other. However, the non-zero amount of homogenized stiffness indicates a significant impact on the structural level, which is in contrast to a self-eliminating inclusion of stress. A similar combination of positive and negative peak values is observed forσ yy andσ xy at different positions, see Figs. 45b and c. While the shear stress componentσ xy is significantly higher than the linear elastic reference value, too, the stresses in vertical direction are negligible. It should be noted, that the combination of tension and compression gives rise to a significant amount of crack driving strain energy density at the position of the tensile peak, that eventually leads to phase-field evolution perpendicular to the initial crack, if the deformation is increased and phase-field evolution is allowed. The results of an additional simulation for crack opening deformation with non-fixed crack orientation and phase-field evolution is shown in Fig. 48a. The lateral growth of the phasefield is located at the position of the tensile peak. Furthermore, the modified crack orientation vector in the elements around the additional phase-field evolution is provided. Instead of an explicit additional branch, the phase-field crack widens and a diffusive pattern of irregular crack orientations is obtained, that leads to further stress inclusions. The crack face sliding deformation gives rise to a very high amount of normal stress in x-direction, see Fig. 45d, while the amount ofσ yy andσ xy is rather small. Nevertheless, the combinations of positive and negative peak stresses in adjacent elements seem to be a characteristic pattern for a deviation angle of φ = 5 • . The deviation angle of φ = 40 • results in homogenized stiffnesses K o = 0.4 × 10 9 N/m and K s = 1.1 × 10 9 N/m for crack opening and crack face sliding deformation, respectively. The homogenized stiffness for crack opening deformation is already at 1.27% of the comparative value. However, the homogenized stiffness for crack face sliding deformation exhibits its maximum value of 13.43% in this case. In contrast to the small deviation angle of φ = 5 • , the combination of tensile and compressive peaks in adjacent elements is not observed for a crack opening deformation. Instead, large peak stresses are obtained for the normal stress in x-directionσ xx , i.e. ≈ 28 times the comparative value, and for the shear stressσ xy , i.e. ≈ 10 4 times the comparative value, see Figs. 46a and c, respectively. Furthermore, the vertical stress in y-directionσ yy amounts to approximately the linear elastic comparative value, see Fig. 46b, which results in additional phase-field evolution for a very small crack opening deformation at the position of the maximum tensile peak, if phase-field evolution is considered, see Fig. 48b. For the crack face sliding deformation, a stress inclusion with a minor tensile stress is observed for the stress componentσ xx , see Fig. 46d, that exhibits a peak value of ≈ −10 4 time the linear elastic comparative value. Furthermore, the vertical stressesσ yy rise to ≈ 1.5 × 10 4 of the comparative value, see Fig. 46e. Additionally, the shear stressesσ xy amounts to ≈ 4.2 of the comparative value. It should be noted, that the deformation applied for crack opening and crack face sliding are similar for Figs. 46 and 45. Nevertheless, the peak stress values of crack driving stress components, i.e. tensile and shear stresses, are significantly higher in the case of an increased deviation angle of φ. Therefore, additional phase-field evolution is triggered at smaller deformation, when the deviation angle is higher.
The deviation angle of φ = 70 • results in homogenized stiffnesses K o = 0.83 × 10 9 N/m and K s = 0.34×10 9 N/m for crack opening and crack face sliding deformation, respectively. Thus, this deviation angle results in the highest homogenized stiffness for crack opening deformation with 2.67% of the linear elastic comparative value. Compared to the previous results, the homogenized stiffness for crack face sliding deformation reduces to 4.15% of the linear elastic comparative value. The horizontal stressesσ xx , visualized in Fig. 47a, and the shear stressesσ xy , shown in Fig. 47c, exhibit reduced magnitudes of ≈ 26 and ≈ −6 × 10 3 times of the linear elastic comparative values. However, Fig. 47 Stress field contour plot for crack distortion angle φ = 70 • with crack opening deformation: aσ xx , bσ yy and cσ xy and with crack face sliding deformation: dσ xx , eσ yy and fσ xy the peak tensile stress in vertical directionσ yy amounts to ≈ 1.9 times of comparative values, see Fig. 47b. While this results in a high amount of crack driving strain energy density in the case of crack opening deformations, it should be noted, that it also coincides with the maximum amount of homogenized stiffness. The result of an additional simulation with phase-field evolution is visualized in Fig. 48c in this case. It should be noted, that in this case the additional phase-field evolves above the initial crack because of the fact, that the tensile peak of the vertical stressesσ yy is above the initial crack, too, see Fig. 47b. The peak stress values of horizontal and vertical stressesσ xx andσ yy , respectively, are both compressive in the case of a crack face sliding deformation. Furthermore, the peak shear stresŝ σ xy is only a relatively small ≈ 1.4 times of the comparative linear elastic value. Thus, a high angle φ of crack orientation deviation is less significant in this case of deformation.

Wing crack evolution
The major difference between a directional split with a crack orientation aligned to the principal direction and an energetically obtained crack orientation occurs for states of strain, where the phase-field driving force is dominated by shear stresses on the crack surface. The experimental investigation of a pre-notched gypsum specimen at compressive loading published in Bobet and Einstein (1998) is an example, where the impact of mode dependent fracture toughnesses can be demonstrated, see e.g. Bryant and Sun (2018). The geometry, boundary conditions and the location of the initial notches are shown in Fig. 49. A horizontal deformation u x is applied on the left edge in negative x-direction in order to apply a compressive loading on the structure with initial cracks . In consequence, so-called wing cracks evolve perpendicular to the initial notches at first. The wing cracks are driven by tensile stresses and exhibit Mode I crack opening deformation. However, due to the change of the load bearing structure by the wing crack evolution, the shear stresses at the tip of the initial cracks increase and result in a secondary crack evolution as a straight extension of the initial cracks. By the application of a mode dependent fracture toughness, Bryant and Sun demonstrate in (2018), that the wing cracks only evolve for a Mode II fracture toughness significantly higher than the fracture toughness for Mode I.
The material parameters are λ = 1.111 GPa, μ = 2.591 GPa, ρ = 2300 kg/m 3 and G c = 50 J/m 2 . The geometry is discretized by 758982 2D triangular elements with characteristic size h = 0.2 mm. The initial cracks are represented by a discrete gap of width 0.1 mm at the location of the initial cracks and the regularization length is specified as l = 1 mm. A horizon-

Fig. 50
Stress contour plot at u x = 50 μm for aσ xx , bσ yy and cσ xy 0.5 mm and time dependency f x according to Fig. 32 is applied at the upper edge of the specimen. The pseudo time step for the static simulation is constant with t = 10 −2 . The crack propagation is investigated for the directional split with energetic crack orientation with n = 1−100. The linear elastic result at displacement of u x = 50 μm is provided in Fig. 50 in terms of contour plots for the stresses with a close up on the edge of an initial crack. Significant shear stresses are observed at the tip of the initial crack in Fig. 50c. In addition, large tensile stresses in x-and y-direction are obtained along the side of the initial crack surface, see Figs. 50a and b. Iso-mode fracture toughness results with n = 1 are provided in Fig. 51a. The larger shear stress at the tip of the initial cracks results in phase-field evolution as a straight extension of the initial crack. In between the two initial cracks, their extensions propagate towards each other and merge. From the structural point of view, the specimen loses its structural integrity due to the shear crack propagating through the whole specimen at ≈ 45 • and a subsequent crack surface sliding deformation. The relation between the reaction force F x and the displacement at the edge of the specimen is provided in Fig. 52. With the shear stress contributions at maximum with n = 1, the lowest peak reaction force is obtained. It should be noted, that the crack orienta-(a) (b) Fig. 51 Phase-field profile and close up on crack tip for a n = 1 and b n = 100 tion is highly distorted around the phase-field cracks and significant convergence problems are observed in the numerical simulation for further phase-field evolution and post-fracture deformation analysis. In contrast, decreasing the amount of phase-field driving force from shear stress contribution by a large Mode II fracture toughness with n = 100 prevents the shear crack evolution. Instead, the typical wing cracks perpendicular to the initial cracks, as observed in the experiments, are obtained, see Fig. 51b. It should be noted, that the Fig. 52 Relation between reaction force F x and displacement u x at x = 152.4 mm highly distorted crack orientation results in a significant widening of the wing cracks, see the close up for u x = 0.53 mm and u x = 0.57 mm in Fig. 51b. In consequence, no further propagation of the phase-field can be obtained due to non-convergence of the results beyond u x = 0.58 mm. Nevertheless, the aspect of wing crack evolution is demonstrated in principle as a consequence of considering a mode dependent fracture toughness in the phase-field model. Furthermore, increasing the Mode II fracture toughness yields a decrease in phasefield driving force for this specific example, where shear stresses are dominated. Therefore, also the peak reaction force increases with increasing n, see Fig. 52.

Conclusion
The phase-field method is a powerful approach with the potential to become the most comprehensive numerical technique for a realistic approximation of fracture evolution and post fracture behavior. The elegant combination of the spatial description of the crack surface with the governing physical principle of crack surface formation, i.e. strain energy dissipation, already provides an approach, that can model and predict the propagation of cracks along realistic and experimentally validated paths. Additional effort is necessary for the further development towards a general description of post fracture behavior close to reality. The analysis and discussion provided in this paper are intended to contribute to this development, both via the development of the directional split as well as providing benchmark simulations for the numerical evaluation of future phasefield approaches.
The following developments and findings are presented and discussed: -A modified irreversibility formulation combines the straight forward implementation of a history variable approach and the recovery of a cracklike phase-field profile after phase-field evolution similar to the application of additional phase-field boundary conditions during the simulation. -Three basic kinematics are formulated for the simplification of an ideal plane and frictionless crack. The compression of the crack faces against each other results in contact force transmission. In contrast, crack opening and crack face sliding deformations are possible without reaction forces. A realistic phase-field model for crack approximation needs to be able to model these features. A suitable approach is based on the identification of stress and strain components normal to the crack surface and shearing components on the crack surface. For components perpendicular to the crack surface, a further distinction between compressive and tensile components is necessary. -Standard phase-field splits, i.e. the volumetric deviatoric and the spectral split, exhibit deficiencies for a realistic post fracture behavior. The introduction of the directional split, that is based on a local vector describing the spatial orientation of the crack, provides a phase-field split, that can approximate the post fracture behavior of an idealized crack close to reality. -The alignment of the crack orientation vector along the direction of the largest tensile principal strain is a suitable approach for most of the strain states. -The amount of crack driving strain energy density depends on the local orientation of the crack orientation vector. The decomposition of the crack driving strain energy density into components related to Mode I and Mode II fracture and a relation with mode dependent fracture toughnesses provide a measure of the phase-field driving force. The phasefield driving force can be formulated with respect to the crack orientation for 2D and 3D strain states. A maximization of the phase-field driving force with respect to the crack orientation under consideration of a mode dependent fracture toughness is derived analytically for a 2D state of strain. -The maximum of the phase-field driving force exhibits a non-unique definition of the crack orientation vector for specific states of strain. A nonunique crack orientation vector is considered as energetically branching. The constitutive behavior is different for the energetically equal crack orientations. Thus, the new crack orientation is specified such, that the previous orientation is at the smallest angle with the new orientation. -The quality of the crack surface approximation in 2D is best, when the crack is straight and the regularization length l is small compared to the size of the domain. For straight cracks, a good approximation is obtained for l ≥ 2h, where h is the characteristic size of the finite elements. The phase-field profile of a crack tip yields an additional artificial amount of regularized crack surface. The impact of this overestimation increases with increasing l and is negligible, if l is small in comparison to the domain and the overall length of the crack. The regularized approximation of branching results in overlapping phase-field profiles and an underestimation of the regularized crack surface. Increasing l yields more overlapping and increased underestimation of the crack surface. Again, the impact of overlapping is negligible, if l is small in comparison to the domain and the overall length of the crack. -For an initial straight crack, the volumetric deviatoric split, the spectral split and the directional split are analyzed with respect to the approximation of the three basic crack kinematics of an idealized crack. The volumetric deviatoric split fails in modeling crack closing deformation in a realistic and numerically stable manner. The spectral split cannot model crack face sliding deformations realistically. The directional split exhibits perfect agreement with the expected behavior, if the crack orientation is constantly perpendicular to the crack surface over the whole phase-field profile. -The directional decomposition of stresses in a 2D element with bi-linear shape functions yields persistent components at the material point level in case of a crack face sliding deformation, if the element edges are not aligned to both the crack surface and the crack orientation vector. A similar defect is expected for 3D simulations with tri-linear shape functions. A solution can be the application of constant strain elements, i.e. triangles in 2D and tetrahedrals in 3D. -The distortion of the crack orientation vector yields artificial stress inclusions for crack opening and crack face sliding deformations, that increase in magnitude with increasing distortion angle. The stress inclusions of tensile and shear kind can trig-ger additional phase-field evolution, that results in a lateral widening of the phase-field profile. In consequence, artificial stiffnesses against the reaction free crack kinematics crack opening and crack face sliding are obtained and the convergence behavior of the solution is affected in a negative manner. -Experimentally observed wing cracks in rock-like material require mode dependent fracture toughnesses G c,I and G c,II , where G c,II must be significantly larger than G c,I . The directional split is suitable for the simulation of wing crack formation in principle. However, the distortion of the crack orientation vector and the evolution of lateral widening for the phase-field profile remain open issues in this case.
Several recent approaches of fracture mechanics make use of an explicit representation of the crack normal, for instance (Morin and Acharya 2021;Storm et al. 2020Storm et al. , 2022. Thus, there is a demand for realistic predictions of the crack orientation in ongoing research activities. The presented concept to determine the crack orientation based on the principle of maximum energy dissipation is supported by the experimental findings in Rozen-Levy et al. (2020).