The concept of representative crack elements for phase-field fracture: Anisotropic elasticity and thermo-elasticity

Summary The realistic representation of material degradation at a fully evolved crack is still one of the main challenges of the phase-field method for fracture. An approach with realistic degradation behavior is only available for isotropic elasticity in the small deformation framework. In this paper, a variational framework is presented for the standard phase-field formulation, which allows to derive the kinematically consistent material degradation. For this purpose, the concept of representative crack elements (RCE) is introduced to analyze the fully degraded material state. The realistic material degradation is further tested using the self-consistency condition , where the behavior of the phase-field model is compared to a discrete crack model. The framework is applied to isotropic elasticity, anisotropic elasticity and thermo-elasticity, but not restricted to these constitutive formulations.


INTRODUCTION
Realistic modeling of load and direction depending material degradation, crack opening, and closure are of fundamental importance for reliable predictions of crack kinematics and crack evolution. Thus, describing these features is still under development for the phase-field method. May et al, 1 Strobl and Seelig, 2 Schlüter, 3 and Steinke and Kaliske 4 have shown, that the well-known volumetric-deviatoric split (V-D) 5 and spectral split 6 approaches with tension/compression decomposition lead to misleading predictions for the force transfer through the crack. They have proposed and developed a model for the crack kinematics in case of isotropic, linear elastic material. The directional decomposition 4,7 overcomes the observed discrepancies of the V-D and spectral split by performing the following steps: • local crack orientation is determined from a chosen criterion, • crack coordinate system is introduced and applied to the stresses and strains, • stress tensor is decomposed with respect to the kinematic considerations of an equivalent discrete crack, and • corresponding material tangent is derived from the decomposed stress tensor.
Also in Teichtmeister et al, 8 Bryant and Sun, 9 and Levitas et al., 10 considerations on the crack kinematics are used with the aim to obtain a consistent material degradation for the phase-field method. In those approaches, a crack orientation is explicitly introduced into the phase-field method. This step requires a few reconsiderations of the interpretation of the phase-field concept in the context of crack irreversibility. 6,11,12 Adopting the idea of a kinematic crack analysis, we propose to introduce a representative model of a corresponding discrete crack. Then, the degradation and crack closure behavior can uniquely be determined for the numerical crack without artificial assumptions like for the above-mentioned split models. Thus, we introduce a kinematic coupling of the phase-field model (continuous crack representation) to a representative model of a discrete crack (discontinuous crack representation) adopting the concept of homogenization. 13 This representative crack element (RCE) allows for the investigations of complex crack behavior (eg, crack surface contact and friction) and to uniquely determine the effective behavior (stress and material tangent) for the continuous crack representation. In contrast to some previous works, the novel concept is formulated consistently in a variational framework, which leads to thermodynamically consistent material approaches and can be applied to inelastic materials and nonlinear geometries (finite deformations). The paper at hand focuses on elastic deformations in the small deformation framework but is by no means restricted to this class of materials. The overall aim of the RCE concept is to obtain phase-field formulations, which can equivalently be used instead of a discrete crack model in any loading situation.
The outline of the paper is as follows. In Section 2, models for the crack kinematics in phase-field fracture are reviewed and analyzed. The concept of RCE is introduced and the variational framework is presented in Section 3. The subsequent two sections address aspects on the crack orientation and irreversibility in the context of the phase-field concept, as well as the convergence behavior. The RCE method is applied to isotropic and anisotropic linear elasticity and thermo-elasticity in Section 6. The implementation of the phase-field material model into a thermomechanical standard element (part of most finite element codes) is sketched in Appendix A. The partial derivatives of the crack projectors are given in the Appendix B.

CRACK KINEMATICS OF PHASE-FIELD FORMULATIONS FOR FRACTURE
The variational formulation of phase-field fracture for the quasi-static case is given by the free-energy potential which consists of mechanical deformation energy mech and fracture energy Γ . In phase-field fracture, the domain of the discrete crack Γ is approximated by a smeared description over the finite volume Ω where ph is the volumetric approximation of the fracture energy via phase-field potential, G c is the critical fracture energy density, s is the phase-field variable and (s, ∇s) is crack density functional. The phase-field variable is limited to the set [0, 1], where s = 0 indicates intact and s = 1 belongs to fully broken material. In some publications, the interpretation of s is vice versa and the formulation of the energy potentials is adjusted accordingly. The material degradation is introduced into deformation energy via the degradation function g(s). The original approach of Francfort and Marigo 14 considers full degradation of the bulk material bulk as With this formulation, crack propagation is predicted also for compression loads, which is unrealistic for most materials. Furthermore, the different material behavior for an opened and a closed cracked cannot be distinguished. Thus, a formulation considering tension-compression asymmetry of the form has been proposed by Henry and Levine 15 , where the deformation energy of the bulk material is additively decomposed into an active + and a passive − part. The evolution equation for the phase-field variable follows from the application of HAMILTON's principle. Obviously, only the active part of the deformation potential contributes to crack propagation. The decomposition into active and passive contributions further depends on the crack state, that is, opened or closed crack. Some existing split models for phase-field fracture define this decomposition of the deformation energy. The two tasks for these models are the identification of the crack state and the decomposition of the deformation energy into active and passive portions. The contribution of the deformation to the phase-field evolution is indirectly defined via Equation (5). Some split models give evolution equations by definition alternatively to Equation (5). Those models are variationally not consistent approaches. The following review of the split models is sorted by model similarities instead of publication date.
The first split is introduced by Henry and Levine 15 as Opened and closed crack states are identified by the volumetric part of the strain tensor (tr E). Here, negative volumetric strain indicates a closed crack and positive counterparts indicate an opened crack. For the closed cracks, a volumetric penalty term is added to the deformation energy using the parameter and the bulk modulus K. The model prevents crack propagation at compression for > 1. All deviatoric deformations drive the crack.
Amor et al 5 and Freddi Royer-Carfagni 16 have introduced the decomposition (VD-split). The VD-split and the spectral decomposition are the two most common split models in phase-field fracture formulations. The crack state is also identified by the volumetric strain. In contrast to the approach of Henry and Levine, 15 the volumetric strain is related to the active and passive portion in dependence on its sign where ⟨⋅⟩± are the MACAULAY brackets and (⋅) dev gives the deviatoric part of a tensor. Thus, negative volumetric strain does not explicitly contribute to phase-field evolution. Lancioni and Royer-Carfagni, 17 and Lancioni and Royer-Carfagni 18 have proposed split models for ashlar masonry structures. The former one make use of the V-D split decomposition, but apply the degradation only to the deviatoric strains The later one is based on the assumption that the fully damaged material response with a negative semi-definite stress, which therefore belongs to the passive portion of deformation energy.
A nonvariational formulation of the V-D split, called the hybrid-split, is proposed in Ambati et al. 19 . The model preserves the linearity of the symmetric mechanical behavior.
but makes use of the volumetric strain to formulate a modified evolution equation for the phase-field variable.
Hence, the unrealistic crack propagation at compressive loads is avoided, but crack closure is not considered in the mechanical response in order to improve the numerical costs.
An alternative for the crack state indicator and for energy decomposition is proposed in Miehe et al. 12 . The spectral-split makes use of the sign of the strain eigenvalues where i and m i are the eigenvalues and eigendirections of the strain tensor, respectively. and are the Lamé constants. Crack opening and closing are assumed to be based on the sign of the eigenvalue. Hence, in this model three independent, orthogonal crack mechanisms are induced. A combination of VD-and spectral-split is proposed in Zhang et al. 20 The VD-split is applied first and afterward the spectral decomposition is applied only to the deviatoric strain.
This approach shows a very similar behavior to the original spectral-split. A cohesive crack surface energy is introduced and coupled to the strain state at the crack surface in Levitas et al. 10 Therefore, the crack normal n 1 is approximated by the phase-field gradient n 1 = ∇s∕||∇s|| and is used to construct the in-plane projector M of the crack surface.
The deformation energy reads where the energy of the bulk material gets fully degraded and the upper term induces an isotropic (bi-axial) plane stress state at the crack surface. This cohesive stress depends only on the critical fracture energy and the in-plane volumetric strain, not on elastic material parameters. In-plane deviatoric strains and stresses are neglected. The nonvariational approach of Zhang et al 21 proposes the modified evolution equation, in order to consider the observed differences in crack propagation at modes I and II crack deformation. The two deformation modes are identified from the spectral decomposition as: where G c,I and G c,II are the critical fracture energies, respectively. Bryant and Sun 9 make use of the crack normal n 1 in order to decompose the deformation energy into modes I and II portions to: with in a two-dimensional formulation. The vector n 2 is the tangential direction on the crack surface. In this model, crack opening is identified from the strain ⟨P 1 ∶ E⟩ + normal to the crack. The influence of deformations tangential to the crack due to Poisson's effect is neglected. The normal direction is determined by applying the principal of maximum energy dissipation, which yields an optimization problem to be solved. The VD-split is extended to finite deformations by Borden. 22 The identification of the crack state is proposed by means of the determinant of the elastic deformation gradient J e Thus, the deviatoric part, controlled by the deviatoric left Cauchy-Green tensor b e , gets fully degraded like in the VD-split. While the isochoric part of the deformation energy is degraded at expansion but not at compressive deformations.
In Hesch and Weinberg, 23 a spectral split for finite deformations is introduced. The decomposed deformation energy reads based on the spectral decomposition of the deformation gradient where m i and M i being the eigenvectors of the left and right Cauchy-Green deformation tensor. The detection of the crack state makes use of the principal stretches along the principal directions. In contrast, the spectral decomposition in Teichtmeister et al 8 is applied to the effective Cauchy stress which yields In all of the aforementioned split models, the crack state is approximately identified by the definition of an artificial condition, for example, sign of volumetric strain, principal strains, principal stresses, or crack normal stress. Moreover, the decomposition of strains and stresses into active and passive portions are artificially introduced into these models, which is shown to predict unrealistic crack propagation and mechanical behavior. [1][2][3][4] Lateral contractions and stresses are neglected in all these criteria.
In Strobl and Seelig 7 and Steinke and Kaliske, 4 an idealized discrete model for the crack with contact is considered ( Figure 5). Isotropic, linear elasticity and small deformations are assumed. Then, basic deformation modes and the corresponding kinematics of the crack are analyzed at this fictitious crack. In both publications, equivalently, closed form solutions for the crack state and the decomposition of the strain/stress tensor with realistic and consistent crack kinematics are presented as directional-split. The strain-based split in Strobl and Seelig 7 is completed by a modified evolution equation for the phase-field variable and belongs to the nonvariational formulations. In Steinke and Kaliske, 4 the decomposition is formulated in terms of stresses in a local coordinate system, where the first axis coincides with the crack normal.
The evolution equations for those splits are consistently derived from Equation (5). The crack orientation is obtained from the phase-field gradient in Strobl and Seelig 7 and the maximum principal stress direction in Steinke and Kaliske. 4 The concept to derive consistent crack behavior from a discrete crack model is adopted and generalized to the computational framework of phase-field fracture with RCE in the following section.

Principals of the crack representation in the phase-field method
Another interpretation of the decomposition in Equation (4) into active and passive mechanical energy can be carried out in such a manner, that the mechanical potential energy is supposed to be the interpolation of the potential energies of the intact material 0 and of the fully broken material c using the degradation function With this clearly defined material states, the additive decomposition of stress and material tangent reads which are composed of the material response without any crack (S 0 , C 0 ) and the material response with a fully developed crack (S c , C ). Once, the two unique material states are identified and described, crack evolution (driving force) is the result of using the novel concept. The intact material behavior can be determined by standard constitutive relations. The fully degraded behavior is determined by the RCE, resolving the entire kinematics of a discrete crack.
For the examples in this paper, the quadratic degradation function and the fracture energy are employed, where l is the length scale parameter and is a very small number to avoid numerical instabilities for fully degraded elements.

Representative crack element and its kinematic coupling
The discontinuous model of a crack (RCE) is kinematically coupled to a material point of the continuous model (phase-field formulation) via (linear) coupling operators. A common choice might be Quantities with bar indicate the continuous model, whereas the rest is associated with the RCE. The displacement field of the RCE u depends on the displacement u, the displacement gradient H of the material point in the continuous model, as well as on the position vector in the RCE , a constant vector 0 and an unknown displacement fluctuation fieldũ. The relations of the displacements and the displacement gradients toward the continuous model are given as their volume average on the RCE domain Ω in the reference configuration. The operators need to be compatible to each other, 13 which yields and the minimal constrained space of kinematically admissible displacement fluctuations where Ω is the boundary of the domain and n is the surface normal vector. The field of displacement fluctuations is the unknown in the RCE. Thus, the variational problem of the RCE is formulated in terms ofũ, where Equation (49) defines the space of admissible solutions. The solution of this problem provides not only the displacement fluctuations but also the stress and material tangent at each material point in the RCE. The relation of the stress field P (1. Piola-Kirchhoff tensor) of the RCE to the corresponding stress of the phase-field model P is given by considering quasi-static deformations. It can be derived from the balance of virtual power of the continuous and the discontinuous model (energy conservation, see Blanco 13 ) and the kinematic coupling operators.u is the rate of the virtual displacements,Ḣ is the rate of the virtual displacement gradient and b is the body force vector. Once the variational problem of the RCE is solved, the stress and the material tangent for the continuous crack model are calculated. In the following, this paper focuses on small deformations using the linearized strain tensor E, the CAUCHY stress tensor S and the material tangent C = S E .

Solution strategies for the RCE problem
The RCE approach is a variational problem regarding the unknown displacement fluctuation field. Several options are available in order to solve the RCE and to use the results for the continuous crack model: • model the RCE as a FEM problem and embed it into the material points of the continuous model (FE 2 approach), • model the RCE as an FEM problem and derive a material model for the continuous model using snapshots of the model behaviour (decoupled approach), • develop an analytical solution for the RCE behavior.
Generality but also computational effort decrease from the first to the last solution method. The development of an analytical or semi-analytical solution is preferred wherever possible, because the computational effort of these models is comparable to ordinary material formulations; whereas the FE 2 method is known to be very expensive.

FEM formulation of the phase-field method with the RCE concept
The potential energy of the continuous crack model is given as This starting point leads after application to the strong form and discretized weak formulation to the linearized FEM approach for an element with the element matrices The element tangent consists of the mechanical tangent K uu , the phase-field tangent K ss and the coupling tangents K us and K su . Also, the reaction vector is assembled by a mechanical contribution R u and a phase-field contribution R s . F u and F s are the external loads. The terms Δu and Δs are the updates of the displacements and the phase-field in the Newton-Raphson iteration, respectively. The vector of the shape functions is N, and its derivatives are given by matrix B.

Test for consistent material degradation
The aim of the RCE concept is to explicitly use the crack kinematics of a discrete crack in order to realistically degrade the continuous crack representation. Then, the correct degradation can be verified by testing the following self-consistency condition. In case, the continuous crack model (phase-field) is used to describe the RCE, the discontinuous and the continuous RCE model are required to give the same reactions (stresses and reaction force) for any applied deformation.

CRACK IRREVERSIBILITY AND CRACK ORIENTATION
In many cases, fracturing is considered to be an irreversible process. Hence, any kind of crack healing is not considered.
In this case, the domain of the crack has to be subset of every subsequent crack domain. 3 Therefore, the crack domain can expand but not transform. Crack irreversibility is often imposed, by Dirichlet constrains on the phase-field variable wherever the material is fully degraded (s → 1) 24 or by a constraint on the fracture energy, which is used for the phase-field evolution ( + ( ) ∶= max ∈[0, ] + ( )). 6 The latter one is shown by Linse 25 how it does not obtain Γ-convergence for the phase-field formulation and overestimates fracture energy. In order to preserve the simplicity of implementation, a combination of both approaches is proposed in Steinke and Kaliske 4 which is applied to the following examples (s c = 0.85) in the paper at hand. The introduced concept of RCE is based on an explicit crack orientation in order to determine the consistent degradation behavior. The introduction of a crack orientation is not unique and requires a concept, also for its evolution and irreversibility. The models in Section 2, which rely on the information about the crack orientation, use one of the following three approaches: • direction of phase-field gradient, 7,8,10 • maximum absolute principle stress direction, 4,9 • direction by principle of maximum dissipated energy. 9 The idea (maximum energy dissipation) of the approach by Bryant and Sun 9 is well known in plasticity models. Further investigations are necessary to develop a consistent variational formulation of this concept in the framework of phase-field fracture. Also nonunique solutions, that is, multiple/each crack orientation is energetically equivalent, and crack branching is of large interest.
The concept of crack irreversibility also affects crack orientation. Obviously, the orientation must not change for a fully evolved crack anymore. In the following applications, the evolution of crack orientation is associated to phase-field evolution. Consequently, the crack orientation is introduced as internal variable, which is updated whenever the phase-field energy evolves ( + ( ) > 0).
The determination and irreversibility of the crack orientation and its influence on crack propagation and computational costs require further investigations. These aspects are out of the scope of the current paper presenting a first introduction of the novel concept of RCE. Therefore, the direction based on maximum principal stress is used for the determination of the crack normal in the subsequent examples, which is in good agreement with experimental observations of brittle, isotropic materials.

CONVERGENCE CRITERIA FOR THE OPERATOR SPLIT METHOD
The presented variational phase-field formulation is separately convex with respect to the deformations and the phase-field variable, but not convex in the generalized degree of freedom, which consists of the displacements and the phase-field variable. This biconvex problem leads to a poor convergence behavior applying standard Newton-Raphson method to the fully coupled formulation (monolithic solution method). However, convergence is observed for problems with damping behaviour, for example, quasi-statics with a viscous regularizations 12 or dynamics with inertia effects. 4 Alternatively, the operator split method (staggered solution scheme) is proposed in Miehe et al 6 for phase-field fracture. This algorithmic decoupling yields two variational subproblem, that is, the mechanical (deformations) and the phase-field. The solution procedure in Miehe et al 6 proposes to solve the phase-field subproblem first and the mechanical problem afterward in each time step. Very small time steps (or an adaptive time scheme) are applied in order to keep the error, caused by operator split method, small. The error yields an underestimation of the crack propagation. A staggered cycle (loop) around the two subproblems is introduced in Ambati et al 19 with the structure • staggered loop: ⚬ phase-field problem, ⚬ mechanical problem, and ⚬ check staggered convergence.
Several energetic convergence criteria for the staggered cycle are proposed and compared by that authors. They are found to be not reliable. At the paper at hand, standard convergence criteria, that is, relative residual norm and relative energy change 26 are applied to the decoupled formulation with The reference values r (1) and e (1) are determined in each time step before the staggered cycle is entered. The results of convergence analysis applying this approach are presented the following section.

APPLICATIONS
The RCE framework for phase-field fracture is not restricted to a certain type of model or structure for the RCE. A simple model is proposed here for the following examples. The model consists of a cube where a discrete crack in the middle plane separates the model into two blocks of intact, homogeneous bulk material with plane surfaces (see Figure 1). A local coordinate system is introduced, where the normal vector of the crack surface coincides with the first axis. The RCE surfaces are divided into corresponding pairs ( Ω -, Ω + ), defined in Figure 2. All surfaces are considered to remain plane and all pairs to remain parallel during deformations. The constraint for the kinematic admissible displacement fluctuations (compare Equation (49)) is fulfilled by the piecewise periodic boundary conditioñ with Hence, both blocks deform homogeneously. Rigid body motions are prevented by applying a point symmetry with respect to RCE center. In comparison to the integral formulation of the coupling operators, periodic boundary conditions has no restriction to generality due to the periodic geometry of the RCE. Also for the contact model, periodicity of the contact surfaces is assumed to prevent that tangential relative displacements of the two blocks lead to loss of the contact partner. This kind of contact model is restricted to either straight cracks or moderate tangential displacements of the crack surfaces. The constraint of plane RCE surfaces induces linear displacement fields in both material blocks. Applying Equation (47) and considering periodic boundary conditions as well as point symmetry yields the relation of the displacement gradient at the RCE H and the displacement gradient of the phase-field formulation H. The components of the displacement gradient toward the crack normal (H 11 , H 21 , and H 31 ) depend on the relative displacements of the crack surfaces u Γ (see Figure 2). Thus, the deformation of the two blocks reads where F is the deformation gradient and Γ = Γ ∕ 1 are the crack deformations obtained from the relative displacements of the crack surfaces u Γ normalized by the RCE length l 1 (see Figure 2). The crack projectors P i are defined as using the basis vectors n i of the local RCE coordinate system. The orientation of the local RCE coordinate system with respect to the global coordinate system needs to be determined by a crack orientation criterion. The unknown crack deformations Γ i can be obtained by minimization of the total free energy of the RCE with respect to the displacement field. The vector n+1 contains internal (state) variables in case of inelastic material for the time point of interest t n+1 . The kinematic coupling with Equations (46) and (64) relates the displacement field of the RCE to the phase-field deformations and the crack deformations for quasi-static problems. Then, the extremal problem with respect to the unknown crack deformations reads, in case there are no constraints on the crack deformations (for an open crack). Subsequently, crack surface penetration has to be tested via Γ 1 < 0 and, possibly, a closed crack solution can be determined from the constrained minimization problem,

Linear, isotropic elasticity at small deformations
The linearized strain tensor E is the symmetric part of the displacement gradient. Then, the strains at the two blocks read where the symmetric part of a tensor is given by: The stresses are decomposed into a part of the intact material and a degradation due to the crack, using the material stiffness C of the intact material. Applying the chain rule for differentiation, also the material tangent for the continuous crack model consists of the intact material tangent and a degradation term Particularly, the energy potential of the RCE for isotropic elasticity of Hooke type readŝ where the material tangent is The solution of the open crack (unconstrained crack deformation) minimization problem yields Hence, the strains in the two blocks correspond to the strain of the continuous crack model modified by a compensation term for the crack normal tension and two compensation terms for the shear strains bridging the crack.
If the condition for a closed crack, holds, then the constrained minimization problem results in The open and closed crack solution only differ in Γ 1 . Thus, the analytical solution of the RCE for isotropic elasticity and small deformations can be obtained from the unconstrained minimization problem via which is fully equivalent to the model of the directional decomposition. 4,7 In case of a closed crack, the stress response is fully isotropically elastic but with compensations for the two shear stresses bridging the crack. For the open-crack case, the normal stress needs to be compensated and the lateral stresses (addressed by the hydrostatic term) are changing. By Equation (84), the mechanical behavior of the RCE can be identified to be bilinear: linear in the deformations leading to an open crack and linear in the relations leading to crack closure. Self-consistency is tested by means of a fully cracked plate (see Figure 4). A penalty formulation is used for the node-surface contact model at the discrete crack. Rigid body modes are suppressed using a point symmetry constraint with respect to the midpoint. The effective strains E are applied to all pairs of surface points ( -, + ) via periodic boundary conditions of the form within small deformation framework. For four representative deformation modes shown in Figure 5, von Mises stress is plotted in Figure 6.   The visualized strips of the plate (cf, Figure 4) are (i) a discrete crack model with contact and the phase-field approach using (ii) the RCE method, (iii) the VD-split, and (iv) the spectral split.
Obviously, only the results by the RCE method (Strip II) coincide with those of the discrete crack (Strip I) for all applied deformations. The predictions of the degradation by the VD-and spectral split fail for several loading conditions, which has been also shown in Strobl and Seelig, 2 Schlüter, 3 and Steinke and Kaliske. 4 The unrealistic degradations for the VD and spectral split yield different results, while modeling a crack with a discrete or a continuous approach. The canonical benchmark of a single-edge notch at simple shear loading (see Figure 7) is usually presented using a discrete notch (eg, Reference 6) instead of a numerical phase-field notch. For a continuous crack model, which passes the self-consistent test, the results using a discrete precrack coincides with those using a continuous crack model for the initial crack (cf, Figure 8). The crack orientation (and consequently the RCE orientation) is modeled using the largest principal stress axes for the crack normal. The load is applied along 30 steps up to the displacement of u y = 0.02 mm. The correct mechanical behavior is further demonstrated for Mode II load of the single-edge notch (cf, Figures 9 and 10). Both models consist of about 23 000 elements.
The presented phase-field formulation is solved by applying the operator split method (staggered solution scheme) introduced by Miehe et al. 6 . The mechanical and the phase-field part of the formulation are separated and solved alternately until convergence for this staggered cycle is achieved. 19 The convergence of the formulation is firstly investigated at the two nested Newton-Raphson loops (one loop for phase-field and one loop for mechanical solution) and then at the staggered loop. The phase-field subproblem (see Equation (45)) is linear and converges in one step. The self-consistency test (see Figures 5 and 6) allows to study the mechanical convergence in case of constant phase-field distribution and crack orientations (final crack behavior). Linear elasticity in combination with tension-compression asymmetry (due to contact) yields bilinear mechanical behavior, which is confirmed by the number of mechanical iterations in Table 1. The benchmark of the single-edge notch at simple shear allows to study the mechanical convergence with variable crack orientations during crack propagation. The mechanical behavior is nonlinear due to the principal axis and, thus, crack orientation can change in each mechanical iteration. The obtained convergence of the mechanical problem is shown in Figure 11. The very first mechanical solution (blue line) requires some more iterations than the subsequent solutions. Most of the mechanical solutions (98%) are achieved using three to five iterations (green lines) and four solutions (1%) have converged within six to seven mechanical iterations (orange lines). The convergence criteria for the staggered loop are based on standard convergence tests applied to the decoupled problem, that is, relative residual norm and relative change in energy (cf, Section 5). The staggered convergence behavior for the single-edge notch at simple shear is shown in Figure 12. Obviously, the number of staggered iterations depends on the crack propagation. Thus, the iteration number is small at the first period of the simulation before the phase-field variable changes significantly. The reduction in staggered iteration numbers correlates with the decreasing speed of crack propagation, which is obtained because the crack grows close to the border of the plate, where strains are constrained by the applied boundary condition.

Linear, anisotropic elasticity at small deformations
The energy potential of the RCE for linear anisotropic elasticity readŝ In order to solve the unconstrained minimization problem of Equation (69), has to hold, which yields the system of equations The inversion of the symmetric 3 x 3 matrix [C Γ ] results in crack deformations by For the case that Γ 1 < 0 holds, the constrained problem reads with The anisotropic solution is demonstrated for a transversely isotropic material defined by the engineering constants E 1 , E 2 , 12 , 23 , G 12 , and G 23 for the compliance tensor in Voigt's notation and the rotation angle around the n 2 -axis. The deformation modes are given in Figure 13 and the results of the self-consistency test (compare Figure 4) are shown in Figure 14. For this material, simple crack closure leads to shear deformations in the crack surfaces. The solution of the discrete model in Figure 14 agrees well with the one of the numerical crack model, but is slightly nonhomogeneous due to the contact formulation. The anisotropic approach is applied to an embedded crack in a thin plate (see Figure 15). The material parameters are E 1 = 209.8 GPa, E 2 = 629.5 GPa, 12 = 23 = 0.3, G c = 2.7 kN/m, and = [70 • , 90 • ]. Thus, the mechanical stiffness in longitudinal direction is three times larger than in transversal direction. As a note, only the anisotropy of mechanical stiffness is considered and not an anisotropic fracture toughness. The model dimensions are l x = l y = 1 m, the crack length is 0.125 m and the length-scale parameter is chosen as 6.64 mm. The discretization contains about 90 000 elements. The load is applied normal to the crack within 100 time increments up to the tensile displacement of u y = 1.0 mm, while the lateral displacements at the plate boundary are constrained.
In Figure 16, crack propagation for the isotropic and the transversely isotropic material is shown. Even though the critical fracture energy is still isotropic, the crack path is changing compared to the one of the isotropic material. This phenomenon is due to different load situations at the crack tip, which leads to different directions for the maximum principal stress. The transversely isotropic material causes crack branching, where the branching angle depends on the material orientation.

Linear, anisotropic thermo-elasticity at small deformations
A thermomechanical material is considered with partially decoupled features. Temperature induces thermal strains, but mechanical deformations do not influence heat conduction. The total strain decomposes additively into an elastic and In analogy to anisotropic elasticity, crack deformations yield where and i, j ∈ {1, 2, 3} for the unconstrained minimization (open crack) or i, j ∈ {2, 3} for the constrained minimization (closed crack). The thermal strains depend on the temperature, which are assumed to be constant at the RCE. In this model, also the temperature drives the crack (via thermal strains) without defining an explicit model for the contribution to the driving force. The temperature change yields two different deformation modes for the precracked plate with isotropic material (compare Figure 17). For both modes, the results of a discrete crack model with contact and those from the phase-field model are compared in Figure 18 for isotropic thermo-elasticity.
The anisotropic thermomechanical RCE approach is applied to the embedded crack as well (see Figure 19). The material parameters are E 1 = 209.8 GPa, E 1 ∕E 2 = [1.0, 3.0], 12 = 23 = 0.3, G c = 2.7 kN/m, th = 12 × 10 −6 1∕K, and = 90 • . Thus, the mechanical stiffness parallel to the crack is up to three times larger then normal to the crack. The model dimensions are 1 m, the crack length is 1 8 of the plate width and the length-scale parameter is 6.64 mm. The discretiszation contains about 90 000 elements. The plate is clamped at the edges and the temperature is proportionally reduced by −20 K along 100 time increments.
The comparison of the crack propagation using isotropic and transversely isotropic material in Figure 20 shows different crack patterns. The different crack paths are caused by different directions of the maximum principal stress axis and the respective crack driving forces. The maximum principal stress field is influenced not only by geometry (crack tip) but also by material anisotropy. Thus, also the crack paths of the two transversely isotropic materials with identical longitudinal directions but different longitudinal stiffnesses result in different crack paths (compare Figure 20C,D).   Figure 22 shows reasonable results. The amount of shear load P S in the first loading stage has a significant influence on the crack path in the experiment as well as in the simulation. The heterogeneous microstructure of the specimens may result in the observed crack branching with smaller side cracks in the experiment. The specimens in the experiment also undergo torsion and/or small Mode III deformation, which results in the difference of the front-and back-side crack path and the asymmetry of upper and F I G U R E 21 Model

Nacre-inspired notched plate at tension load
The nacre-inspired structure for a notched plate at tensile load is presented in Jeong et al 28  The length of the notch is 0.24 mm. The critical fracture energy for both components is G c = 5 × 10 −5 kN/mm and the length-scale parameter l = 0.0004 mm. The element size is 0.002 mm (about 770 000 elements) and the displacement load is Δu y = 5 × 10 −2 mm. The irreversibility approach of Miehe et al 6 is applied, which allows a better investigation of the crack branching behavior from the final crack pattern. The geometry chosen in the paper at hand differs slightly from that in Jeong et al, 28 where at the plate edges the structure is cutted in a nonperiodic manner. Moreover, instead of a two-dimensional simulation a three-dimensional analysis with a plate thickness of 0.006 mm is performed, where plain strain conditions are enforced by boundary conditions. In Figure 24, the crack paths obtained with the isotropic RCE model are compared for three different supports at the plate given in Figure 23. It becomes obvious from the chosen deformation states that in all of these loading situations, a secondary crack is initiated which merges with the crack starting at the notch. This secondary crack initiation is also found for the phase-field model with VD-split and with no split. In order to study the given microstructure, the dimensions of the plate need to be much larger then the dimension of the components to reduce such boundary effects. For the two cases with symmetric boundary conditions (cf, Figure 23A,B), crack branching points possibly occur with equivalently alternative crack paths. Thus, the final crack path in this bifurcation problem is the result of small numerical model disturbances. The crack patterns, obtained from the simulation, are not suitable for comparison of the results. The given model neglects a possible anisotropy of the critical fracture energy in the soft material layers and that the fracture energy possibly differs in the soft and hard components, as well as at the interface between both. Undergoing shear in the soft material layers, the crack is deflected toward the hard component for the crack orientation criterion based on the principal stress axis. The missing interface barrier and interface delamination yields in crack propagation through the hard aggregates in some cases. In order to ensure cracks in the soft component only, the crack orientation is prescribed in the soft layer by the layer normal. Crack orientations at the intersections and in the hard component are obtained by the maximum principle stress criterion.
The final crack patterns of the plate with the support at the right edge, using the phase-field formulations with no split, with VD-split, with spectral split and with the isotropic RCE, are shown in Figure 25. The results obtained with the spectral split differs significantly from those of the other three models. The corresponding reaction force depicted in Figure 26 also differs. This is because of the wrong prediction of the spectral split formulation for Modes II and III crack deformations, which are induced by the given microstructure. spectral split RCE no split VD-split

CONCLUSIONS AND OUTLOOK
Standard phase-field formulations for fracture with VD-split or spectral decomposition of the mechanical potential energy yield unrealistic predictions for the crack kinematics and, thus, for the material degradation. [2][3][4] A consistent degradation model is presented in Strobl and Seelig 2 and Steinke and Kaliske 4 under the restriction to isotropic elasticity and small deformations. The paper at hand presents a variational framework for phase-field fracture, which allows to derive realistic and kinematically consistent material degradation without such restrictions. For this purpose, the concept of RCE is introduced to analyze the fully degraded material state. The kinematic coupling between the phase-field model and the RCE is given in a variational framework adopting the homogenisation concept of Blanco. 13 Realistic degradation behavior can be tested with the self-consistency condition, in which the behavior of a discrete crack and the phase-field model are required to be equal for any applied deformation. The concept is applied to isotropic elasticity, anisotropic elasticity and thermo-elasticity. In order to avoid huge computational effort of a nested scheme (like in FE 2 ), analytical solutions are derived for the used RCE model. The assumptions for the analytical solutions are the neglect of softening for the bulk material (to prevent localisation) and moderate tangential displacements of the crack surfaces (to simplify the contact model).
The obtained degradation model for isotropic elasticity is equivalent to those developed by Strobl and Seelig 2 and Steinke and Kaliske. 4 The comparison to a discrete crack simulation shows the same results. This is in contrast to phase-field simulations using V-D or spectral decomposition. To demonstrate that the phase-field models using the novel RCE concept are equivalent to discrete crack models, the well-known single-edge notch at simple shear and Mode II load is shown with an initial phase-field crack instead of a discrete crack with equal results. The mechanical and the staggered convergence behavior are depicted as well.
Self-consistency is also analyzed for transversely isotropic material and thermo-elasticity. The results are in good agreement with those of discrete crack models. Also, tangential crack surface displacements for compression loads normal to the crack in the anisotropic model and crack opening/closure due to thermal strains are correctly predicted. The influence of the anisotropic elastic behaviour on the crack path is studied at a thin plate with an embedded crack at uniaxial deformation and thermal shrinkage.
The promising results and the simple implementation in case of analytical RCE models indicate the great potential of the phase-field method for fracture also for finite deformations and dissipative materials, which will be presented in a future publication.

APPENDIX A. IMPLEMENTATION OF THE PHASE-FIELD METHOD INTO A THERMO-MECHANICAL STANDARD ELEMENT
A short comparison of the thermomechanically coupled element formulation to the phase-field element above shows the same structure [ K uu K ut with the element matrices Hence, the phase-field method can easily be implemented into finite element codes with a thermomechanical standard element. In this case, the whole phase-field method can be implemented into a user-defined, thermomechanically coupled material routine.

B1. Derivatives of the projectors
The application of chain rule yields with the eigenvectors (principal axes) n i of the stress tensor and the symmetric dyadic product a⊗ s b = 1 2 (a ⊗ b + b ⊗ a). Using the relation for the eigenvector derivative 29