Phase-field modeling of fracture in heterogeneous materials -- jump conditions, convergence and crack propagation

In this contribution, a variational diffuse modeling framework for cracks in heterogeneous media is presented. A static order parameter smoothly bridges the discontinuity at material interfaces, while an evolving phase-field captures the regularized crack. The key novelty is the combination of a strain energy split with a partial rank-I relaxation in the vicinity of the diffuse interface. The former is necessary to account for physically meaningful crack kinematics like crack closure, the latter ensures the mechanical jump conditions throughout the diffuse region. The model is verified by a convergence study, where a circular bi-material disc with and without a crack is subjected to radial loads. For the uncracked case, analytical solutions are taken as reference. In a second step, the model is applied to crack propagation, where a meaningful influence on crack branching is observed, that underlines the necessity of a reasonable homogenization scheme. The presented model is particularly relevant for the combination of any variational strain energy split in the fracture phase-field model with a diffuse modeling approach for material heterogeneities.

the element edges do not have to be aligned with the crack path, i.e. remeshing is avoided when the crack changes its direction or branches.
The approach for a non-conforming description can be extended to a more general setting, where all discontinuities within a heterogeneous structure are captured in such a way. Cumbersome preprocessing such as manual identification of heterogeneities and meshing can be skipped or at least simplified if the structure is for example obtained from direct imaging [7,8,9]. Prominent representatives are the extended finite element method [10,11,12] or the finite cell method [13,14]. The latter is among the category of fictitious domain methods, where the geometry is fully embedded within a larger domain [15,16,17].
Another way to incorporate material heterogeneities is closely related to the phase-field model for cracks. A (multi)phase-field, which introduces a diffuse region Γ i , with a characterizing width length scale i , along the formerly sharp interface Γ i , marks the individual subdomains, and offers a broad range of applications in the field of thermodynamics, chemistry and mechanics [18,19,20]. The individual phase-field parameters are often referred to as order parameters. Far away from the interface, the properties of the constituents are recovered. Within the diffuse region, additional assumptions such as the treatment of strong discontinuities of the stress or strain are required, which will be highlighted in the following.
In classical homogenization theory, the Voigt and Reuss limits play a major role. Transferring the corresponding assumptions of a constant strain or a constant stress within multiple constituents to the diffuse interface region implies, that either the strain or stress does not exhibit a jump and the corresponding quantities are equal for each phase. It can be shown, that these assumptions converge towards the sharp interface limit [18]. However, better approaches exist, which account for the stress and strain discontinuity across the interface simultaneously [17,21,22,23]. The improvement manifests itself in better convergence rates and the fact, that the kinematics of an energy driven, moving interface are captured correctly [24]. The improved homogenization scheme is often referred to as partial rank-I relaxation, which will become apparent in the remainder of the paper, and yields a pointwise fulfillment of the equilibrium and compatibility within the diffuse region.
We presented a modeling framework [25], where a crack phase-field and a static order parameter were used to simulate crack branching and deflection in heterogeneous media. Besides a qualitatively and quantitatively good agreement with analytical predictions from linear elastic fracture mechanics [26], we observed inconsistencies of the modeling results when comparing the diffuse interface approach to a sharp interface simulation. We presumed, that the Voigt-Taylor homogenization approach, which was incorporated in the diffuse region, biased the crack phase-field driving force in a way, that a straight crack instead of a deflection at the interface was favored. As discussed above, the homogenization scheme significantly influences the energetic driving force of the phase transition.
This work extends our diffuse modeling framework [25] in a way, that a partial rank-I relaxation is incorporated within the diffuse interface region. The key novelty of the present paper is the variationally consistent combination of the aforementioned homogenization approach with an additive strain energy decomposition, which is necessary to account for physically meaningful crack behavior such as the tension compression asymmetry. It will be shown, that the novel scheme is superior to the classical Voigt-Taylor assumption.
The paper is structured as follows. Section 2 provides a variational formulation of the modeling approach. Subsequently in Section 3, a convergence study and crack propagation simulations provide an insight into the functionality and advantages of the model. Final conclusions summarize the paper in Section 4. Additional information concerning technical details are provided in Appendices A to D.  The investigations presented in this contribution build on a fully diffuse framework for crack propagation simulations in heterogeneous materials. The crack is described following the phase-field approach. Suppose, a one-dimensional rod x ∈ (−∞, ∞) of infinite length is cracked at x = 0. The sharp crack location can be fixed using the D distribution, cf. Figure 1a. In contrast to a sharp crack, the phase-field approach introduces an additional scalar field with a characteristic length c , cf. Figure 1b, where c = 1 and c = 0 resemble intact and fully broken material, respectively. Consequently, a three-dimensional generalization for the corresponding surface energy of the crack [3,27] can be derived, which is a measure for the energy required to form the corresponding crack surface. The fracture toughness G c stems from the energetic cracking criterion by Griffith [6]. In other words, the crack surface has been regularized [5]. In case of a homogeneous fracture toughness, Equation (1) is obtained by rewriting Equation (2) for the one-dimensional case and subsequent minimization of the functional, subject to the boundary conditions c(0) = 0 and c (x → ±∞) = 0. For a heterogeneous fracture toughness, phase-field profiles different from Equation (1) are obtained, cf. [25] for more detailed considerations. In Figure 2a on the left, a phase-field crack Γ c is depicted, which possibly interacts with the interface between the two subdomains Ω 1 and Ω 2 . Next, the diffuse representation of these subdomains is outlined.
inclusion of arbitrary shape in matrix The diffuse modeling approach is exemplarily depicted in (a), where a heterogeneity is embedded in a regular mesh. The dashed interface mid-surface is described by the level-set d(x) ≡ 0 of the signed distance function. The sharp interface Γ i is smoothed by a hyperbolic tangent transition (b), and becomes a diffuse region Γ i .

Diffuse modeling of elastic heterogeneities
In the present work, two linear elastic domains Ω i with i = 1, 2 are considered. A static order parameter p ∈ [0, 1] describes for every material point whether it belongs to the one Ω 1 or other Ω 2 subdomain, and introduces a smooth transition in the vicinity of the former sharp interface Γ i , cf. Figure 2a. The transition width is controlled by the interface length scale i . The order parameter does not change in time and is a priori known, e.g. derived from a direct imaging technique. Finally, the motion is described by the displacement u ∈ R 3 . The elastic energy density of each domain corresponds to where is a quadratic degradation function to account for previously introduced phase-field cracks. Various other degradation functions are proposed in [28,29], including multiple functions for isotropic and anisotropic contributions [30] and compared for example in [31,32]. The residual stiffness η 1 prevents numerical problems for fully degraded material c = 0.
Furthermore, only the positive part i ψ el + of the strain energy density is degraded to account for the tension compression asymmetry of cracks. Note, that positive does not refer to the sign of the strain energy density but to a one-dimensional analogy, where only positive stresses are degraded. A correct and physically meaningful, yet simple and efficient, splitting of the strain energy has been subject of discussion of several publications, cf. for example [27,33,34,35,28,36,37,38]. In case of large deformation kinematics, similar approaches can be found [39,40,41,42,43]. In this contribution, the so-called tensile split [27,Section 3 with the L constants i λ, i µ and M brackets • ± = 1 2 (• ± | • |) is adopted. For further information on the strains i ε ± , please refer to the cited section. It is noted, that the diffuse interface model is not restricted to the tensile split, but can be applied to any decomposition of type (3).
The diffuse interface region with p ∈ (0, 1) in the vicinity of the formerly sharp interface Γ i is referred to as Γ i . The order parameter p, which recovers the values p = 0 for Ω 1 and p = 1 for Ω 2 far away from the interface, cf. Figure 2a, varies according to which corresponds to a three-dimensional generalization of the one-dimensional solution of the M -M functional [44], cf. Figure 2b. The signed distance d(x) measures the shortest distance from any point in Ω 1 ∪ Ω 2 to the interface mid-surface in the direction of the interface normal n i = ∇p/|∇p|, cf. Figure 2a on the right. In this contribution, the signed distance d(x), and thus p(x) are a priori known. Within the diffuse region, the elastic energy density is defined as which also holds for the individual phases, because one or the other contribution vanishes for p = 0 or p = 1. The approach of interpolating the individual constituents is a widespread approach, cf. for example [18,19,20,21,24]. An alternative representation can be obtained by inserting Equation (3) into (7). The individual dependencies on i ε are omitted above for the sake of readability. Analogously to Equation (7), an interpolation of the individual strains is assumed in the diffuse region. Together with the strain jump ε = 2 ε − 1 ε, the independent variables for the strain energy density are now changed from which is more convenient in view of the presented model. The strain is obtained from the displacement u by whereas different homogenization assumptions for ε are discussed in Section 2.4. The next section combines the energy contributions from the phase-field crack and the linear elastic deformation, and derives the coupled field problem.

Global incremental potential and governing differential equations
Starting from the findings and definitions of the previous section, the internal energy density is presented, which is used in the following derivations of the governing differential equations. It contains bulk material and crack phase-field contributions, ψ el and ψ c respectively. The global incremental potential (14) is defined according to [45,Equation (90)]. The index n refers to the previous, converged time increment, whereas quantities without the index are related to the current (n + 1) th increment. The time increment is denoted by τ = t n+1 − t n . It is noted, that a first order E backward time integration scheme is employed for the rate c of the crack phase-field, which enters the Here, η f ≥ 0 is a kinetic parameter, often referred to as viscosity, which is of purely numerical nature in this contribution. Its choice is discussed in Section 3. The incorporation of c can be understood as a viscous regularization, which makes it possible to use a monolithic solution approach. Various other approaches and their implications can be found in literature, see e.g. [2,46]. The vectort denotes given external tractions on N boundaries ∂Ω t . The strain jump, which enters ψ, is determined within a local minimization procedure, cf. Section 2.4, and is therefore known which is indicated by the star ε . Hence, u and c are obtained by the global minimization principle with admissible sets The first constraint in C is only relevant for the initial crack along Γ c 0 . Afterwards, the second condition, also known as fracture-like irreversibility constraint [3,47], suffices, where c th 1 is a small, positive threshold value. Previous investigations [25] confirmed, that setting c th = 0.03 does not impact the crack path, while ensuring irreversibility. A discussion on different irreversibility constraints can be found in [47]. Subsequently, the E -L equations can be derived, subject to N boundary conditions where n b denotes the outward normal vector along the corresponding boundary. Equation (4) has been used in Equation (19). The thermodynamically consistent relation for the C stress tensor is defined by It is noted, that the chain rule has been used together with Equations (10) and (11). The next section is dedicated to the determination of the strain jump ε which enters the incremental potential (14). The validity of Equation (22) is demonstrated in Appendix A.

Homogenization assumptions within the diffuse interface region
In this section, two different approaches from the homogenization theory are examined to determine the unknown strain jump ε , which enters the incremental potential (14) as ε .
In particular, the Voigt-Taylor approach and a partial rank-I relaxation are introduced. In contrast to existing works [20], the combination of an additive strain energy decomposition and the above approaches is presented. This work incorporates the tensile split by Miehe [27] but is not limited to it.
Voigt-Taylor approach: In classical homogenization theory, the Voigt-Taylor approach assumes equal strains in every constituent, i.e.
From that, the strain jump ε = 0 follows, which is inserted in the incremental potential above. In other words, the Voigt-Taylor approach is automatically implemented in any diffuse modeling framework similar to this one, if the strain jump is not treated by some means or other. Different from the assumption of equal strains in the constituents, the Reuss-Sachs approach assumes equal stresses. It is, however, not investigated in this contribution, because there are various other publications doing so [21,22,24]. The two approaches yield upper and lower energetic barriers for the energy, respectively, cf. Partial rank-I relaxation: The previous approaches do not include any information about the interface but only interpolate between the bulk energies, cf. Equation (7). In contrast, the partial rank-I relaxation correctly incorporates all mechanical jump conditions by means of the strain jump, which can be expressed by the H condition [48, Equation (2.2.9)] without loss of generality1, and the stress jump σ · n i = 0. The former is stated here for a symmetric tensor ε. The vector a is the strain jump amplitude and n i = ∇p/|∇p| 1 The H condition holds for a continuous deformation field and is also referred to as kinematical compatibility condition. A proof is provided in [48,Sec. 2.1.6]. Because of the regularized description of the crack, the displacement field is always continuous. Thus, the condition also holds when a crack meets the interface.
the interface normal. The latter is fulfilled as follows: In order to calculate a, a pointwise minimization procedure [21,Equation (50)] with the associated necessary condition is pursued, which is a relaxation of ψ el with respect to a [24]. The stresses i σ ± are calculated according to Equation (22). For the cases p = 0 or p = 1, which indicate bulk material away from the diffuse interface region, the condition is fulfilled and the jump ε = 0 either way. For p ∈ (0, 1), the expression in the curly braces has to be zero, which is similar to the aforementioned stress jump condition at sharp interfaces. Here, a pointwise enforcement throughout the whole diffuse interface region is required. In absence of a crack, g(1) = 1, the analytical solution for linear elastic material according to [24,Equation (42)] is recovered because the additive strain energy decomposition is unnecessary. The same is valid for g(c) < 1 if there is no strain energy decomposition at all, cf. [20], because g(c) does not exhibit any jump across the interface. In this contribution, the key novelty lies in the incorporation of any (variational) additive decomposition of the strain energy into a degraded and persistent part. The nonlinearity of the spectral decomposition for the tensile split requires a local N -R scheme to solve for a as soon as a crack emerges within the diffuse interface region. In the context of an elasto-plastic multiphase-field model, Schneider et al. [19] had to solve a very similar equation to (26) because of the nonlinearity due to plasticity. The local N -R scheme is presented in Appendix B.

Numerical implementation
The weak form including consistent linearization for the global N -R scheme can be found in Appendix C. The spatial discretization is carried out using locally refined Truncated Hierarchical non-uniform rational B-splines [49] (NURBS) with quadratic shape functions. An adaptive refinement strategy is pursued, which allows for efficient computations with a high resolution of the steep gradient in regions where the crack develops and propagates. Additionally, the mesh is prerefined along the interface Γ i to resolve the diffuse transition for p sufficiently. Spatial convergence with respect to the crack length scale c is ensured by choosing the finest element level in a way, that the characteristic element size is at least three times smaller than c , cf. Figure 3. All simulations are carried out using a Matlab-based in-house finite element code.

Convergence study
For verification of the implementation, the convergence of the presented modeling approach with respect to the sharp interface solution is tested. For this purpose, a quarter of a circular bi-material disc is subjected to a radial displacement u R = r 2 /1000, cf. Figure 3. In the study, the simulated domain is reduced to a square. The geometry parameters are given as r 1 = 3 mm, r 2 = 15 mm and a = 8 mm. Plane strain is assumed. The signed distance for Equation (6) reads and total error in energy Investigation without a crack At first, the model is tested for the linear-elastic case in the absence of a crack. The material parameters are set to ν 1 = ν 2 = ν = 0.3, E 2 = 100 GPa, c = 50 µm and a very high value G c = 10 5 kN mm −1 to prevent any crack formation or propagation. The viscosity η f is set to zero. Two cases E 1 /E 2 = {2, 20} allow to study the influence of the intensity of the elastic dissimilarity. The coarsest element level of the mesh is depicted in Figure 3 on the bottom left. An h i -refinement study is undertaken, where the interface length scale and the mesh size are reduced simultaneously, i.e. more refinement levels are added along the arc with radius r 1 . The ratio i /h ≈ 2.7 is kept constant.
The results of the convergence study are depicted in Figure 4. In general, it can be seen, that the partial rank-I scheme exhibits the same or better convergence rates and error levels than the Voigt-Taylor approach with respect to the analytical sharp interface solution. Furthermore, the obtained results are consistent with literature, compare Figure 4a to [17, Figure 8b], or Figure 4b to [24, Figure 5b]. Additionally, the Voigt-Taylor approach reacts more sensitive to a more intense elastic dissimilarity which can be seen when comparing the left and right plots in Figures 4a and 4b. This is also reported in [18].
Investigation with a crack across the interface Next, a phase-field crack along Γ c 0 with e = 0.5 mm is introduced as shown in Figure 3, with the corresponding finite element mesh on the bottom right. The mesh is prerefined to the finest level along the crack to rule out any inaccuracies due to approximation errors of the phase-field. It serves as coarsest refinement level and is gradually refined along the interface for the h i -convergence study. All parameters and the boundary conditions are chosen similar to the previous investigation, except that only results for E 1 /E 2 = 2 are considered for the sake of brevity. It is noted, that the tractiont does not reproduce the analytical solution as before, nor does it lead to a circumferentially constant radial displacement u R , which is why an overkill solution is used as reference. For the overkill solution, the NURBS mesh has to account for the physical C 0 -continuity of the displacement along the arc r 1 . Additionally, the crack boundary condition along Γ c 0 has to be set. This is why the finite element mesh is different from the mesh of the previous convergence study. Convergence of the overkill solution has been verified.
The results for three different values of the residual stiffness η, cf. Equation (4), are depicted in Figure 5. All plots are identical, i.e. the residual stiffness does not influence the convergence, which is expected. However, for crack propagation simulations, a zero value most probably leads to numerical problems. Furthermore, the obtained results are more or less identical to the intact case, cf. Figure 5 to the left plot in Figure 4b: Convergence rate and error level are matching.
The two convergence studies suggest the successful verification of the implemented method and validation by means of literature results. The major question whether any impacts on propagating cracks can be observed is answered in the next section.

Crack branching at interfaces
Besides the improved convergence of the presented model with respect to the sharp interface solution it is of great interest whether any impacts on crack propagation at diffuse interfaces are recognized. Hansen-Dörr et al. [25] thoroughly investigated many different setups, where a crack approaches a possibly inclined interface and, depending on the ratio of the elastic moduli on both sides of the interface and the ratio of the bulk material and interface fracture toughnesses, decides to branch or to go straight into the adjacent bulk material. A comparison to analytical results from He and Hutchinson [26] and to a sharp interface reference revealed inaccuracies of the Voigt-Taylor approach, which was used therein. This section demonstrates, that the presented model extension to a partial rank-I relaxation yields the same critical deflection ratios as the sharp interface reference, while the Voigt-Taylor gives different results.
The setup for this study, which is similar to [25], is depicted in Figure 6. The corresponding NURBS mesh is shown in Figure 7a. A square domain with an initial phase-field crack Γ c 0 is subjected to a horizontal displacement load. For this purpose, the concept of the so-called surfing boundary condition [50] is exploited, which proved to be very suitable for quantification purposes [50,51]. The key idea of this approach is to introduce a point with the time dependent positionx, which is moving along the y-axis A displacement of hyperbolic tangent-like shape is applied on the side edges ∂Ω u,s , with respect to the moving point, assuming u ref = 7.5 µm, d = 0.5 mm, v = 0.3 mm s −1 and y 0 = −4 mm. Similar parameters proved to be suitable in previous studies [25]. Homogeneous N boundary conditions are considered along the top and bottom edges. A horizontal interface, implied by the grey hatched bar in Figure 6, divides the domain in two halves. The fracture toughness of the interface is different from the bulk material. This is accounted for  In total, seven hierarchical levels are used. The simulation results for α = −0.5 and G b c /Ĝ i c = 23 are depicted in (b) to (d). For these three plots, the two black lines mark the interface midline with d ≡ 0. Shortly after meeting the interface, branching is predicted by the sharp reference solution. In case of a diffuse interface, the partial rank-I relaxation predicts the same phenomenon, while the Voigt-Taylor approach fails. The visibility of elements with c < 0.1 is turned of which is why the contour legend stops at 0.1. by locally reducing the fracture toughness according to a G -like distribution cf. Figure 6 blue plot in the middle. Far away from the interface, the bulk material fracture toughness G b c = 2.7 N mm −1 is recovered. Along the interface, a potential crack is supposed to propagate for a critical energy release rate equal to G i c . However, the crack experiences a bulk material influence due to the interaction of both length scales i and c [25], and tends to propagate at critical energy release rates, which are higher compared to the minimum of the function in Equation (32). Hence, a compensated interface fracture toughnessĜ i c is specified in a way, that a crack propagates along the interface midline at the true interface fracture toughness G i c , i.e.Ĝ i c < G i c . For more details on the compensation procedure, the reader is referred to [25,52,53]. The ratio of the elastic moduli with E 1 = 210 GPa varies, while the P ratio ν = 0.3 is kept constant, i.e. the interpolation in Equation (7) directly translates into a hyperbolic tangent function for E 1 and E 2 , cf. orange line in Figure 6 on the right. For the sharp interface reference solution, the elastic moduli vary according to the purple step function. The elastic dissimilarity is described with the first D ' parameter [54] α = Following the study of Hansen-Dörr et al. [25], the length scale parameters are set to i = 18.75 µm and c = 15 µm, while the viscosity and the residual stiffness take the values η f = 10 −5 kN s mm −2 and η = 10 −5 , respectively. A parameter study with α = {−0.5, −0.25, 0.25, 0.5} and different ratios G b c /G i c is conducted, and the results for the sharp interface and the diffuse interface description using the Voigt-Taylor approach or partial rank-I relaxation are compared with respect to the crack behavior at the interface: Fixing α, the tendency of the crack to branch is higher for larger ratios The NURBS mesh is adaptively refined according to the crack path. All elements, where the phase-field fell below a value of c = 0.5, are marked for refinement [49]. For the sake of brevity, only the results where the two homogenization approaches yielded different crack phenomena at the interface are presented.
Figures 7b-7d present the crack patterns shortly after the crack met the interface, with For better interpretation, the elements with c < 0.1 are made invisible. It can be seen, that the sharp interface solution and the diffuse interface approach with a partial rank-I relaxation predict a branching of the crack. On the contrary, the Voigt-Taylor approach yields a straight crack with no such tendency. Even 3, did not change this fact. This is in line with the presumption of Hansen-Dörr et al. [25], who also observed such deviations. For further loading, all cracks started to propagate further in y-direction, while the two crack branches for the sharp interface and the diffuse partial rank-I relaxation arrested. For α = −0.25 and , a similar behavior is observed. For α = {0.25, 0.5}, the three approaches always yield the same outcome. On the one hand this could imply an increased sensitivity of the crack behavior towards the homogenization scheme for negative α values. On the other hand, only discrete values of G b c /G i c are investigated and it is possible, that only those were picked from the infinite amount of possible ratios, where the results are identical. This is, why future model investigations will include a deeper and also three-dimensional analysis of the phenomena reported above. To conclude with, it was shown, that the present model extension meaningfully influences the crack behavior and yields closer results to the sharp interface solution.

Conclusions
In this contribution, a variational diffuse modeling approach was presented. Heterogeneities were described by an order parameter which is obtained from the signed distance function and a hyperbolic tangent, and smoothly bridges material discontinuities. The mechanical boundary value problem is monolithically coupled to a phase-field model for cracks.
Due to the diffuse transition between two materials, the mechanical jump conditions are not necessarily fulfilled for a general case. Instead, without taking any measures, the Voigt-Taylor homogenization approach is implemented within the diffuse region, which states equal strains in both phases. As shown above, this fails to reproduce several crack patterns compared to the sharp interface. A remedy to this issue is the partial rank-I relaxation, which was extended to account for a strain energy density decomposition. The decomposition is in general necessary for physically meaningful crack behavior like closure. The partial rank-I relaxation yields a point-wise fulfillment of the mechanical jump conditions and was successfully validated by means of convergence investigations. The presented scheme includes a local N -R loop, which is the main drawback of the approach: Although, the local algorithm always converged without any problems, it is more time consuming than avoiding it by employing the Voigt-Taylor approach. A final study of a crack propagating towards an interface demonstrated, that the partial rank-I relaxation does not only improve the convergence of the global energy, but also influences the crack growth locally towards the sharp interface solution since it alters the crack driving force directly. The Voigt-Taylor approach did not achieve this for the given length scale of the interface, which is in line with observations from literature, and underlines the necessity to fulfill the mechanical jump conditions.
where 1 σ and 2 σ correspond to the stresses in the individual phases. Subsequently, the strains of phase 1 and 2 are differentiated with respect to the interpolated strain using Equations (10), (11) and (24), and one obtains where δ i j denotes the K symbol. Inserting Equations (35) and (36) into Equation (34) yields In case of the Voigt-Taylor approach, ε ≡ 0 and the derivatives of ε vanish. For the partial rank-I relaxation, ε depends on the strain state ε, cf. Appendix B, which is why the chain rule is applied. The corresponding derivative is obtained by exploiting Equation (24). Hence, the last summand in Equation (37) vanishes, because Equation (26) holds. This has been confirmed in combination with the local N -R scheme B, which yielded deviations from zero in the order of 10 −20 . Finally, is obtained, which is equal to Equation (22).

B -Local N -R scheme for the strain jump
A local N -R scheme is employed in order to determine the strain jump amplitude a, for which holds, cf. Equation (26). For this purpose, the local residual of the ith local iteration is defined. Index notation is used for better readability of the tensor products. It is noted, that the stresses of the individual phases comprise the degraded and persistent stress contributions. A first order Taylor series expansion of the local residual with respect to a variation of the strain jump amplitude yields The incremental increase of the strain jump amplitude can be obtained by rearrangement. The local tangent reads The specific choice of the strain energy split influences the stress tangents i C l moh of the individual phases. Hence, the scheme can be applied to any split, which fits in the present model.
For the consistent linearization in Appendix C, the derivative ∂a o /∂ε l m is needed. It is obtained, by implicitly differentiating the converged local residual with respect to the strain. Together with Equations (35) and (36) and some rearrangement, is obtained, which finally gives

C -Weak form and consistent linearization
The weak form of the coupled differential equations (18) and (19) reads where δε = 1 2 (∇δu + (∇δu) ) = ∇ s δu, δu and δc denote the variations of the strain, the displacement and the phase-field. The total domain is subsequently divided into n el finite elements for which, the primary field variables and corresponding variations are approximated by with quadratic shape functions N α , where α is the global node number, and the corresponding nodal value (•) α . The weak form (48) has to hold for arbitrary variations δu and δc, which is why the residuals of the ith iteration and (n + 1)th time step are defined. All variables with no time index refer to time step n + 1.
A first order Taylor series expansion of the residuals with respect to a variation of the primary field variables gives In order to determine the incremental increase of the displacement and the phase-field for the (n + 1)th time step, the system of equations is solved for the global node number n nd in every iteration i. It is noted, that the global stiffness matrix is symmetric due to the variational structure of the problem. The four submatrices are given as  Please refer to Equations (37), (38), (44) and (47) for the evaluation of the two above equations. For better readability of the tensor products, the index notation is used. It is noted, that the evaluation of ∂ 2 a 0 /∂ε l m ∂ε i j is not necessary, because the stress jump condition is pointwise fulfilled for p ∈ (0, 1).
The radial and tangential stress components are obtained from σ r r (r) =Ẽ for the square domain of unit thickness, and is used as reference for the convergence investigation of the total energy for the case without a crack. The stress components in Cartesian coordinates read σ x x (x, y) = σ r r (r) cos 2 ϕ + σ ϕ ϕ (r) sin 2 ϕ , σ y y (x, y) = σ r r (r) sin 2 ϕ + σ ϕ ϕ (r) cos 2 ϕ , σ x y (x, y) = σ r r (r) sin ϕ cos ϕ − σ ϕ ϕ (r) sin ϕ cos ϕ , with the angle ϕ = arctan(y/x) and radius r = x 2 + y 2 . The tractiont along ∂Ω top and ∂Ω right for the convergence investigation with and without a crack is obtained from σ x x (x,ȳ) σ x y (x,ȳ) σ y x (x,ȳ) σ y y (x,ȳ) · n b =t with {x,ȳ } = {x, y | x, y ∈ ∂Ω top ∪ ∂Ω right } .