Analysis and simulations for a phase-field fracture model at finite strains based on modified invariants

Phase-field models have already been proven to predict complex fracture patterns for brittle fracture at small strains. In this paper we discuss a model for phase-field fracture at finite deformations in more detail. Among the identification of crack location and projection of crack growth the numerical stability is one of the main challenges in solid mechanics. Here we present a phase-field model at finite strains, which takes into account the anisotropy of damage by applying an anisotropic split of the modified invariants of the right Cauchy-Green strain tensor. We introduce a suitable weak notion of solution that also allows for a spatial and temporal discretization of the model. In this framework we study the existence of solutions and we show that the time-discrete solutions converge in a weak sense to a solution of the time-continuous formulation of the model. Numerical examples in two and three space dimensions illustrate the range of validity of the analytical results.


INTRODUCTION
In solid mechanics, one of the main challenges is the prediction of crack growth and fragmentation patterns. Regarding the modeling side, complicated structures and a non-regular behavior of cracks turn numerical simulations into a difficult task. The classical brittle fracture approach of Griffith and Irwin [1,2] is based on an energy minimization setting for the whole structure. Let us consider a solid with domain  0 ⊂ ℝ 3 and boundary  0 ≡ Γ ⊂ ℝ 2 deforming within a time interval ∈ [0, ]. Each crack, that is located in a solid, forms a new surface Γ( ) of a priori unknown position which needs to be identified. Therefore, the total potential energy of a homogeneous but cracking solid is composed of its bulk energy with a Helmholtz free energy density Ψ and of surface energy contributions due to growing cracks: For rate-independent brittle fracture the material's resistance to cracking, the fracture toughness  , corresponds to Griffith energy release rate with unit force per length. Rate-independent crack growth sets on as soon as the energy release rate  reaches the critical value  . The evolution law encoded therein at time amounts to a minimization of the energy (1) under the constraint that Γ( ) ⊂ Γ( ) for all < ∈ [0, ]. However, the minimization of the energy functional (1) is a challenging task because of the moving boundary Γ( ). Several sophisticated discretization techniques exist, e.g., cohesive zone models [3][4][5], eroded finite elements [6,7], or eigenfracture strategies [8,9] to name some of them. Another approach to such moving boundary problems is a diffuse-interface approach, which approximates the two dimensional crack surface by a three dimensional damaged volume. These types of phase-field models for fracture have gained much attention in recent literature, cf. e.g. [10][11][12][13][14]. The main idea of this ansatz is to mark the damage state of the body by a continuous order parameter ∶ [0, ] × ℝ 3 → ℝ, which evolves in space and time. From a mathematical point of view the phase-field fracture model we shall investigate in this work, is a modification of the Ambrosio-Tortorelli functional [15], which can be used to model rate-independent volume damage, and which we here augment by a viscous dissipation potential for the variable . The modifications made allow it to consider the evolution problem at finite strains using polyconvex energy densities. They take into account the anisotropy of damage, meaning that damage only increases under tensile loadings but not under compression. For this, it will be important to consider the stored energy density as a function of the modified principal invariants of the right Cauchy-Green strain tensor, and based on this, to introduce an anisotropic split of the energy density, which we shall explain in Section 2.
This model is analyzed in detail in Section 3. We introduce a suitable notion of solution and in this setting we prove the existence of solutions using a time-discrete scheme via a minimizing movement approach. In particular, we show that solutions constructed with a staggered time-discrete scheme converge in a weak sense to a solution of the time-continuous formulation of the problem. This convergence result is confirmed within a series of numerical examples in Section 4, where we also provide further details on the spatial and temporal discretization. We demonstrate a simple but typical problem of a mode-I-tension test in two and three dimensions to study different influencing factors.

THE PHASE-FIELD FRACTURE APPROACH
In this section the focus is set on the phase-field approach for fracture to overcome the difficulty of moving boundaries. The phase-field is introduced as an additional parameter which is by definition a continuous field. Thus, the moving crack boundaries are "smeared" over a small but finite length. The order parameter ∶  0 × [0, ] → ℝ with ∈ [0, 1] characterizes the state of the material, whereby = 1 indicates the unbroken state and = 0 the broken state. The surface integral in (1) is approximated by a regularization using a crack density function ∶ ℝ × ℝ 3 → [0, ∞) This approximation (2) is inserted into (1) so that the total potential energy reads The crack density function is by definition zero along the cracks only. In general different approaches, e.g., a second order or fourth order crack density function can be applied, cf. Figure 1. Typically, a second-order phase-field approach is defined as with the fixed parameter ∈ (0, 1) which is a measure for the width of the diffuse interface zone, see Figure 1. The lengthscale parameter weights the influence of the linear and the gradient term whereby the gradient enforces the regularization of the sharp interface. The insertion of (4) in (3) leads to a potential which corresponds to the Ambrosio-Tortorelli functional, cf. [15]. F I G U R E 1 Uniaxial model with a crack at = 0 and with a continuous phase-field ∈ [0, 1]; phase-field approximation for a second order and a fourth order approximation of the crack density function

Governing equations
The elastic boundary value problem is based on the balance of linear momentum div( ) + = in  0 .
Here and in the following denotes the first Piola-Kirchhoff stress tensor, is the prescribed body force and will be the traction. In this contribution the focus is set on a simple elastic boundary problem such that the body force will be neglected. The solid's boundary is divided into displacement and traction boundaries  D 0 and  N 0 with  D 0 ∪  N 0 =  0 ,  D 0 ∩  N 0 = ∅, and = ( ) on  N 0 .
with the outward normal and fixed prescribed displacements ( ). All fields refer to the reference configuration. Furthermore, the phase-field evolution equation is stated in general form aṡ where the parameter denotes the kinematic mobility and summarizes all driving forces which typically represent a competition of bulk material and surface forces, cf. [16]: In particular, the phase-field model is based on the crack density function (4) and the driving force for crack growth that consists of components of the free energy.
Here it has to be mentioned that our analytical results also account for the unidirectionality of the damage evolution, i.e., that damage of the material can only increase during the evolution, but not heal. This can be done by reformulating Equation (8) for the damage evolution as follows: Under the assumption that the kinematic mobility is strictly positive we can equivalently rewrite (8) as −1̇= − . We then introduce (the density of) the quadratic dissipation potential  . According to the definition = 1 for the unbroken and = 0 for the completely broken state of the material, it thus ensures that the damage of the material can only increase in time, which means the unidirectionality of the damage evolution. Since the unidirectionality constraint is not incorporated in the numerical simulations presented in Section 4, we use the prefactor = const ≥ 0 to indicate that we switch this constraint on or off, so that we can consider two different types of models: A model with = const > 0, where the unidirectionality constraint is active, and a model with = 0, where unidirectionality is not incorporated (i.e. 0 ⋅ ∞ = 0) and which is used in the simulations. Our analytical results cover both cases = const > 0 and = 0. Using (10) we see that the evolution equation (8) for = 0 is given by Ḋ(̇( )) = − ( ) in  0 . (11a) In the case = const > 0 the dissipation potential is non-smooth foṙ= 0. Then equation (11a) is replaced by the subdifferential inclusioṅ(̇( featuring the multivalued subdifferential of the convex potential . More precisely, In Section 3 we will introduce a suitable weak formulation for the evolution problems given by the two cases (11a) and (11b).

Finite elasticity and the anisotropy of damage
In this work we set the focus on finite strains and we will subsequently introduce a nonlinear material model and its split into compressive and tensile parts, which makes sure that only the latter are responsible for crack growth.
In the finite deformation regime a deformation mapping ∶  0 × [0, ] → ℝ 3 is regarded and the deformation gradient Concerning the following notation the fields in capitals refer to the reference configuration. Furthermore, the volume map is given by the determinant of the deformation gradient, det . Further, the assumptions of hyperelasticity, objectivity, and isotropy of the constitutive law imply that the free energy density can equivalently be written as a function of the principal invariants of the right Cauchy-Green strain tensor = ⊤ . In three space dimensions, they are given by The cofactor cof maps the element area vector from the reference to the current configuration. Since many materials behave quite differently under bulk and shear loads, it is often convenient for numerical simulations to introduce a multiplicative decomposition of into a volume-changing and a volume-preserving part: with (det ) 1∕3 being the volume-changing deformation,̄is the volume-preserving deformation and ∈ ℝ 3×3 is the identity matrix. This split goes back to Richter [17] and was also successfully used by a series of other authors, cf. [18][19][20][21]. This leads to a set of modified principal invariants 1 1 Please note that the notations and do not stand for the left and right polar decomposition here. Throughout this contribution and will always denote the modified invariants defined in (15).
which in fact guarantee a stress-free reference configuration since cf. (88). The free energy density Ψ may now be expressed either in terms of the invariants ( 1 ( ), 2 ( ), 3 ( )) or in terms of the modified invariants ( ( ), ( ), 3 ( )). As we shall see, the use of the modified invariants is more convenient. The form of the free energy density is further specified by taking into account the anisotropy of fracture. Cracks grow and damage increases only under tensile deformations but not under compression. To incorporate this feature, we first renormalize the stored energy density of the reference configuration. Here we assume that the reference configuration is associated with the state ( , ) = ( , 1), = 1 marks the undamaged state of the body. We now impose that ( , 1) = 0. For this, we note that 1 ( ) = ( ) = 2 ( ) = ( ). Thus, the condition ( , 1) = 0 can be achieved if (⋅, 1) is composed e.g. of renormalized power terms of the form ( ( ⊤ ) − ( )) or ( ( ⊤ ) − ( ) ) with ≥ 1 and being a placeholder for the invariant functions 1 , 2 , 3 , , . Moreover, compressive deformations are characterized by ( ⊤ ) < ( ), while tensile deformations are associated with ( ⊤ ) ≥ ( ). Hence, the anisotropy of fracture can be incorporated into the model by making the specific ansatz or, respectively, where we introduced ± ( ( ⊤ )) = ± max{±( ( ⊤ ) − ( ) ), 0} for ≥ 1, = 1, 2, 3, and for the two sets of invariants ( 1 , 2 , 3 ) ∈ {( 1 , 2 , 3 ), ( , , 3 )}. In the discussion below we will also use … as a placeholder for˜and .
The degradation function is defined with the specific ansatz where the parameter > 0 is a very small value ≪ 1 to catch numerical instabilities in the cases of = 0. In fact, from a mathematical point of view the introduction of > 0 ensures the coercivity of (⋅, ) for any ∈ [0, 1]. Thus, at level > 0, only incomplete damage is modeled, because the material is still able to carry stresses even if = 0.
The split with ± is explicitely tailored for energy densities of power-law type alike with ≥ 1, and ≥ 1. When using a given density … to introduce˜±( 1 , 2 , 3 ) and ± ( , , 3 ) as in (16), there must not hold˜±( 1 , 2 , 3 ) = ± ( , , 3 ), since the densities˜± and ± use different sets of invariants to which the functions ± are applied. Because of premultiplication with ( ) the two definitions of Ψ in terms of either˜or are different and lead to different results. In fact, for numerical simulations the use of has turned out to be more stable [22].
In order to fulfill the assumption of hyperelasticity we also want to make sure that the density … ± ∶ ℝ 3 × ℝ 3 × ℝ 3 → ℝ is continuously differentiable with respect to the invariants ( 1 , 2 , 3 ). It can be checked that this holds true for energy densities … of the form (18) with > 1, and ≥ 1. In this way, when focussing on the case = 1, the densities … ± defined in (16), are given by … ± ( 1 ( ⊤ ), 2 ( ⊤ ), 3 where summation now runs over the 3 invariants. Then it is and lim ℎ→0 (±( ( ) ± ℎ) ∓ ( ) ) −1 = 0 thanks to > 1. Hence, each of the derivatives … ± ∶ ℝ + × ℝ + × ℝ + → ℝ is continuous with … ± ( 1 ( ), 2 ( ), 3 ( )) = 0. This ensures that the reference configuration (associated with the deformation gradient = ) is stress-free. In particular, for continuously differentiable densities … ± ∶ ℝ 3 × ℝ 3 × ℝ 3 → ℝ, due to the assumption of hyperelasticity and relation (16), the first Piola-Kirchhoff stress can be calculated with the chain rule as and ( , ) = 0 for any ∈ ℝ because … ± ( 1 ( ), 2 ( ), 3 ( )) = 0. In contrast, for = 1 in (19) one can see that ± is no longer continuously differentiable in ( ) since the left and the right differential quotient in ( ) do not coincide. Then, the choice of the set of invariants has an influence on the continuity properties of the first Piola-Kirchhoff stress. Exemplarily, let us consider the density of a Neo-Hooke material ( 1 ( ⊤ ) − 3) in the split (16), i.e., 1 , 1 = 1, 1 = and 2 = 3 = 0 in (19). Using (16) we have˜( , ) = ( ) ( + 1 ( 1 ( ⊤ ) − 3)) + ( − 1 ( 1 ( ⊤ ) − 3)). For ≠ , the densitỹ (⋅, ) is continuously differentiable so that the first Piola-Kirchhoff tensor takes the form For = the energy density is non-differentiable, since max{ 1 (⋅), 3} has a kink in , but the left and right limits of ( , ) do exist for = , and for ( ) = 1 these two limits coincide. However, a discontinuity of ( , ) in = cannot be prevented if ( ) ≠ 0, i.e., when damage occurs. Instead, when using the Neo-Hooke law with the modified first invariant in the split (16), i.e., ( , ) = ( ) ( + 1 ( ( ⊤ ) − 3)) + ( − 1 ( ( ⊤ ) − 3)), we find that with lim | |→0 ( ± ) = 0 and with ± ∈ material with modified invariant ( ⊤ ) is continuous with ( , ) = 0 for any ∈ ℝ, which is not the case when the principal invariant 1 is used. For numerical simulations it may be more convenient to introduce the anisotropic split for energy densities of power law type in a slightly different way, as proposed in [22]. While the split defined in (16) and (19) is incorporated into the terms of the free energy density, the anisotropic split proposed in [22] is directly imposed on the modified invariants; more precisely, the anisotropically splitted invariants are defined by with being a placeholder for the invariants , , 3 . In Figure 2 the split of the invariants is visualized. We note that This is due to the fact that the split (24) is again tailored to densities of polynomial type as in (18). In particular, for of the form (18) we can check that ± ( , , 3 ) = ( ± , ± , ± 3 ). Therefore, the above discussion about the continuity of the first Piola-Kirchhoff tensor remains valid also when the anisotropic split is formulated in terms of (24): Continuity of (⋅, ) in holds true for exponents > 1 in (18) for both sets ( 1 , 2 , 3 ) and ( , , 3 ). For = 1 the densi-ties˜±( 1 ( ⊤ ), 2 ( ⊤ ), 3 ( ⊤ )) = … ( ± 1 ( ⊤ ), ± 2 ( ⊤ ), ± 3 ( ⊤ )) are non-differentiable in ( ⊤ ) = , so that stress (⋅, ) has a discontinuity in , which cannot be prevented for ( ) ≠ 1. Instead, when the modified invariants are used, (⋅, ) is continuous with ( , ) = 0. To manifest this statement we will consider a simple mode-I-tension test in a loading and unloading regime. We regard a two dimensional plate with a horizontal notch. The geometry and the related boundary conditions are depicted in Figure 3. On the lower boundary of the plate the displacements are constrained in horizontal and vertical direction and on the upper side prescribed displacements are applied incrementally. The displacements are increased until the time = 0.5 sec is reached -after that the plate is relieved to the reference configuration. Furthermore, the mesh presented in Figure 3 is on the basis of the hierarchical refinement strategy, see [23], and consists of 20 × 20 quadratic B-spline elements before making use of the refinement. After three local refinement levels in total 2656 elements with overall 12288 degrees of freedom are employed for the tension test.
The simulation is based on the non-linear Neo-Hookean material model with the proposed anisotropic split (24), here considered for two dimensions, The material parameter are chosen as follows: shear modulus = 80.7692 × 10 9 N m 2 , bulk modulus = 121.1538 × 10 9 N m 2 and specific fracture energy as  = 2.7 × 10 3 J m 2 . The process of loading and unloading is illustrated by the load-deflection  Figure 4. The focus is set on the stresses = (det ) −1 during the simulation after different time steps, cf. Figure 5. The snapshots of stresses demonstrate that the configuration is stressfree after unloading again. That means, the stresses are continuous with ( , ) = 0.
Our analysis uses the formulation (16) to take into account the anisotropy of damage. It is carried out directly for the density ( , ), hence it is independent of the set of invariants, but it relies on the assumption that the energy density is continuously differentiable.

ANALYTICAL SETUP, DISCRETIZATION, AND CONVERGENCE RESULT
For the mathematical analysis of the regularized crack model given by Equations (5) and (8) maps a 3 × 3-matrix onto the vector of its minors. We remark that the ansatz (∇ , ) ∶= ( ) 1 ( ∇ ) + 2 ( ∇ ) used in (27) is in accordance with the setting proposed in (16). In particular, one may choose 1 = + and 2 = − from (16). In this way, the anisotropy of damage can be reflected in the model, given that the densities are sufficiently smooth. The assumptions on the densities ∶ ℝ 3×3 × ℝ 3×3 × ℝ → ℝ are specified more precisely in (36) in Sec. 3.1. For the existence analysis, carried out in Sections 3.1-3.3, it will be more convenient to regard the densities directly as functions of the minors of gradients, and not as functions of the modified invariants of the right Cauchy-Green tensor as introduced in (15). Instead, in Section 3.4 we will translate the assumptions (36) imposed on the densities into assumptions for densities being functions of the modified invariants. In accordance with (10) we additionally introduce the viscous dissipation potential  ∶ → [0, ∞) with −1 the inverse of the kinematic mobility from . According to the definition = 1 for the unbroken and = 0 for the completely broken state of the material, it thus ensures that the damage of the material can only increase in time, which means the unidirectionality of the damage evolution. With the prefactor = const ≥ 0 we indicate that we switch this constraint on or off, so that we can consider two different types of models: A model with = const > 0, where the unidirectionality constraint is active, and a model with = 0, where unidirectionality is not incorporated (i.e., 0 ⋅ ∞ = 0). The latter case is considered in the phase-field flow rule (8) and used for the numerical simulations presented in Section 4.

Notion of solution for the body with damage
The elastic body undergoing damage is thus characterized by a suitable state space × , the energy functional  from (27) and the dissipation potential  from (29) and we refer to it as the (evolution) system ( × , ,  ).
In fact, we will find in Lemma 3.10 that D ( , ( ), ( )) ∶= ′ ( ( )) ( ∇ ( )) + ′ ( ( )) − Δ ( ) is bounded only in * the dual of the space ∶= 1 (Ω) ∩ ∞ (Ω). For = 0 we thus find condition (30b), which is the weak formulation of the phase-field equation (103). Formulation (30c) for > 0, i.e., when the unidirectionality constraint is active, is given in terms of a one-sided variational inequality. Together with the energy-dissipation estimate (30d) and (30a) it provides a characterization of the solutions for > 0. Such a formulation of the non-smooth damage evolution in terms of a one-sided variational inequality combined with an energy dissipation estimate has been applied in [26] in the small-strain setting, and later also e.g. in [27,28] in the case of small-strain (visco)plasticity.

Main result and analytical strategy
Our main result, Thm. 3.11, provides the existence of a solution ( , ) ∶ [0, ] → × of ( × , ,  ), which satisfies the governing Equations (30). Its proof will be carried out in Section 3.3 via a time-discretization. For this, we will consider an equidistant partition of the time interval In addition, for the analysis, we will also regularize the unidirectionality constraint (−∞,0] by its corresponding Yosida approximation, i.e., for each ∈ ℕ we introduce where (̇) + ∶= max{0,̇} is the positive part oḟ. Starting out from an admissible inital datum ( 0 , 0 ) ∈ × , at each time-step ∈ Π we then alternatingly solve which is also called a staggered time-discrete scheme and used for the simulations in Section 4. Existence of solutions ( , ) of (33) at each time-step will be shown in Prop. 3.9. With the discrete solutions ( , ) =1 we will construct suitable interpolants with respect to time and show in our main result, Thm. 3.11, that these interpolants approximate a solution of the continuous problem (30).

Comparison with other results in literature
Our evolution system ( × , ,  ) combines the energy functional  from (27), which is a modification of the Ambrosio-Tortorelli functional [15], with a quadratic dissipation potential  , cf. (29), which thus causes a viscous evolution of the phase-field variable. Without this viscous contribution, i.e., −1 = 0 in (29), the (standard) Ambrosio-Tortorelli functional, combined with the unidirectionality property ensured by > 0, models the rate-independent, unidirectional evolution of the phase-field variable. In this setting, it was shown in [29] that the standard Ambrosio-Tortorelli functional Γ-converges to the Francfort-Marigo model for brittle, Griffith-type fracture, cf. e.g. [30] as the diffuse-interface parameter → 0 in the definition of , cf. (4). Later, similar approximation results have been obtained in the rate-independent setting at small strains, allowing for the use of the linearized strain tensor and for modifications of the energy functional leading to cohesive or elasto-plastic fracture models, cf. [31][32][33][34][35][36]. Instead, the limit passage → 0 of an energy functional of Ambrosio-Tortorelli-type in combination with a viscous evolution of the displacements was investigated in [37].
In this work we do not consider the limit → 0. We rather study for > 0 fixed the existence of solutions in the sense of (30) for the system ( × , ,  ) at finite strains for a modified energy functional that takes into account the different evolution behavior of a viscous-type phase-field variable with regard to tensile or compressive loads. Due to this modification, since the energy contribution 2 , which accounts for compression, is not premultiplied by the function , we do not expect our model to approximate the Francfort-Marigo fracture model. Instead, we understand ( × , ,  ) as a very specific model for partial, isotropic damage, which has the property to localize damage in zones of width . This viewpoint on the Ambrosio-Tortorelli model has also been taken in the rate-independent setting e.g. in [38], where the alternate minimization scheme (30) has been further iterated in each time-step leading to parametrized -evolutions af the rate-independent problem, and in [39], where also a viscous regularization has been taken into account.
Other models for partial, isotropic damage allow for more general forms of the function and the regularization . While the very specific properties of the functions , in (27) make it possible to show that a solution satisfies ( ) ∈ [0, 1] a.e. in  0 for a.e. ∈ (0, ), cf. Prop. 3.9, this bound is in general lost for other physically reasonable choices of and . Then, additional indicator terms have to be incorporated into the energy in order in order to ensure that ∈ [0, 1] for physical and analytical reasons. For the analysis of general models for partial, isotropic damage with a rate-independent damage evolution we refer to the works, e.g., [40][41][42][43] and to the monography [44] for an overview on rate-independent processes. The viscous, rate-dependent counterpart is studied, e.g., in the works [26,[45][46][47], also in combination with dynamics, heat transport, and phase separation, and vanishing-viscosity limits from viscous damage models at small strains to rateindependent ones are investigated in the series of works [48][49][50].
Our approach to the analysis of system ( × , ,  ) extends the ideas used in [51], based on [52], for a rate-independent damage model at finite strains, to the present viscous setting by making use of the notion of solution studied in e.g. [26] at small strains. The study of the properties of energy densities as functions of the modifed invariants, cf. (15), builds on results drawn from the works [25,[53][54][55][56][57].

Analytical setup: Assumptions and direct implications
A physically reasonable deformation preserves orientation, which is ensured by Further natural requirements on the constitutive relations of particular importance are material frame indifference (34a) and the non-interpenetration condition (34b): However, they are not compatible with the convexity in the deformation gradient , which is a convenient claim in the setting of small strains. Instead, they can be compatible with the convexity in , cf. [58]. To see the incompatibility with convexity in the deformation gradient consider , ∈ SO(3), ∈ (0, 1), such that ( + (1 − ) ) ∉ SO (3), which conforms to a strain. Then convexity together with material frame indifference yields the following contradiction: The class of energy densities which fit to these natural requirements and which admit to prove existence are polyconvex energy densities. They were introduced by J.M. Ball in [53].
is the function, which maps a matrix to all its minors.
In [53], p. 362] it was established that the polyconvexity of̂∶ ℝ 3×3 → ℝ implies its quasiconvexity. By C.B. Morrey in [59 it was proven that quasiconvexity is the notion of convexity which is necessary and sufficient for the lower semicontinuity of the corresponding integral functionals, so that quasiconvexity together with other technical assumptions ensures the existence of minimizers. But quasiconvexity does not admit infinitely valued functions, i.e.̂∶ ℝ 3×3 → ℝ ∞ . However in [53, Th. 7.3, p. 376] it was shown that the polyconvexity of the densitŷ∶ ℝ 3×3 → ℝ ∞ together with other technical assumptions is sufficient for the existence of minimizers of infinitely valued functionals.

• Stress control:
For all ∈ℝ we have (⋅, ) ∈ C 1 (GL + (3), ℝ) and there are constants >0,̃≥ 0 such that for all ( , )∈ℝ 3×3 ×ℝ it holds • Uniform continuity of the stresses: There is a modulus of continuity ∶ [0, ∞] → [0, ∞], > 0 so that for all Herein, assumptions (36b)-(36d) ensure the existence of minimizers, see [60], p. 182, Thm. 2.10] and the discussion in Remark 3.6 below. In fact, the cases a), b), and c) provide three alternative sets of relations for the exponents , 2 , 3 building on compactness results for minors of gradients given in [53,54,60], see Remark 3.6. In analytical works on evolution problems for generalized standard materials often (the 2D-or D-version of) assumption (36d) a) is used, cf. e.g. [24,40,52,[61][62][63], since this ensures the continuity of the deformation and avoids ambiguities in the definition of its minors, see Remark 3.6. Within cases b) and c) we slightly weaken these assumptions in accordance with the results of [53,54,60] in order to include energy densities into our analysis, which are popular in engineering. We also refer e.g. to the works [64,65], respectively [66][67][68], where (the 2D-or D-versions of) assumption (36d) b) has been applied in the context of rate-independent fracture, resp. for the Γ-limit passage from finite-to small-strain linear elasticity for rate-independent processes in generalized standard materials. Furthermore, assumptions (36e) and (36f) ensure that a minimizer of (30a) satisfies the (weak form of the) corresponding Euler-Lagrange Equation (102) for all ∈ 1 (ℝ 3 ; ℝ 3 ) such that and D are uniformly bounded and satisfy ( ) = 0 on the Dirichlet boundary in trace sense. Moreover, there is > 0 so that for all

Assumptions on the domain, state spaces & given data
As in (27) we consider a body with reference configuration  0 ⊂ ℝ 3 consisting of a nonlinearly elastic material, such that This body undergoes a damage process driven by time-dependent exterior forces ( ) located on the Neumann part  N 0 ⊂  0 of the boundary. Moreover, the body is assumed to be clamped at the remaining part  D 0 of its boundary, so that the deformation is prescribed there: ( ) = ( ) on  D 0 . Thus, the set of admissible deformations at time ∈ [0, ] is given by with the weak 1, -topology. By the assumption on in (36d) we have in particular that > 3∕2 and thus, in three space dimensions, Adapting the ideas of [61] from the setting where > 3 to the present setting ∈ [2, ∞), we assume that the Dirichlet datum can be extended to ℝ 3 in the following way: For the time-dependent Neumann datum we impose that To handle the time-dependent Dirichlet conditions one assumes that the deformation is of the form with the weak 1, -topology. Regarding (6) this means that ( , ) = ( , ( )) = ( , ) on  D 0 . By the chain rule, the composition (45) leads to a multiplicative split of the deformation gradient: Under consideration of (27) and the explanations along with (30) we choose the set of admissible damage variables in (27) and the set of admissible test functions 0 in (30) as equipped with the respective weak topologies. The sets and form the state space × , which is endowed with the weak topology of the product space. For the closed subspace 0 ⊂ 1, ( 0 , ℝ 3 ) Friedrich's inequality is available: There is a constant = ( 0 , ) such that the following estimate holds for every 0 ∈ 0 ∶

Bounds and convergence properties for time-dependent Dirichlet data
We now provide bounds and convergence properties in the spaces ( ) and , cf. (41) and (46), which follow from the relations for the Dirichlet datum (43) and (45). First of all, the lemma below is a consequence of the growth restriction (43c). (40), (43) as well as (45) hold. For every ∈ and ( ) = ( , ) it holds

Lemma 3.4. Let
Proof. By the growth restriction (43c) one directly obtains Due to the realization of the Dirichlet condition in terms of a superposition, cf. (45), thanks to the regularity properties (43) of the Dirichlet datum the convergence of a sequence in translates to the convergence of a sequence respecting the Dirichlet condition in ( ) as follows:

Compactness of minors of gradients according to (36d)
We first motivate our sets of coercivity assumptions (36d) a), b), and c) in a remark.
Remark 3.6 (Results on compactness of minors). In three space dimensions the minors of a matrix ∇ are given by ∇ itself, its cofactor matrix cof ∇ , and its determinant det ∇ . Hereby, the entries of cof ∇ contain products of two derivatives , whereas det ∇ is composed of terms given by products of three derivatives . In this spirit, it is shown in [53, Cor. 6.2.2] that the map ↦ cof ∇ ∶ 1, ( 0 ) → ∕2 ( 0 ) is weakly sequentially lower semicontinuous if > 2, whilst the map ↦ det ∇ ∶ 1, ( 0 ) → ∕3 ( 0 ) is weakly sequentially lower semicontinuous if > 3. As done in [53, p. 370] the functions cof ∇ , det ∇ can given meaning in the sense of distributions by defining suitable distributions denoted as Cof ∇ , Det ∇ , and in general cof ∇ and Cof ∇ , resp. det ∇ and Det ∇ must not coincide. We point out the result [53, Thm. 6.2], which relaxes > 3 to the condition > 3∕2, 1 + 1 < 4∕3 and ensures: If ⇀ in 1, ( 0 , ℝ 3 ), then Cof ∇ ⇀ Cof ∇ and Det ∇ → Det ∇ in  ′ ( 0 ), i.e., the price is that the defined distributional minors cannot be identified with the minors as a function. Yet, as previously discussed, this ambiguity can be ruled out if > 3. Hence, taking a look at the coercivity assumptions (36d), it is therefore sufficient for > 3, cf. (36d) a), to claim boundedness of the energy density from below only with regard to |∇ | ; compactness of the other minors then follows by [53,Cor. 6.2.2]. In contrast, in the case of (36d) b), ∈ [2, 3) is possible, so that ∕3 < 1, and the weak sequential lower semicontinuity of the determinant cannot be concluded directly. Similarly, if = 2 compactness cannot be concluded for bounded sequences of cofactors. This is why coercivity assumption (36d) b), which is taken from [60, Thm. 2.10, p. 182], requires bounds with exponents larger than one for all of the three minors. In particular, the set of assumptions is designed exactly in such a way that, thanks to [60, Thm. 2.6, Part 5, p. 173], ambiguities between cof ∇ and Cof ∇ , resp. det ∇ and Det ∇ can be ruled out, even though < 3 is admissible. Moreover, assumptions (36d) b1) and b2) are tailored to the case that the energy density does not depend on det ∇ , resp. neither on cof ∇ nor on det ∇ , so that compactness in these terms is not needed. In this way our analysis also includes e.g. Neo-Hooke and Mooney-Rivlin materials, cf.

Analytical results on the convergence of discrete solutions
In this section we gather and explain all our analytical results.
For the interpolants ( , , , , ) we then verify a discrete version of the governing equations (30) and uniform apriori estimates. In addition, the following energy-dissipation estimate holds true for all ∈ ℕ, ∈ [0, ]: where the partial time derivative ( , ( , ( )), ( )) is given by (52). Furthermore, there is a constant > 0 such that the following apriori estimates are satisfied uniformly for all ∈ ℕ and all ∈ [0, ]: Thanks to the apriori estimates (61) we are now in the position to extract a (not relabeled) subsequence ( , , , , ) of the interpolants, which converge to a limit pair ( , ) that satisfies (30): for any sequence → in ( 0 ) for some ≥ 1 with 0 ≤ ≤ 1 a.e. in  0 .

Proof of Prop. 3.9
In order to establish the proof of Item 1, we will employ the direct method of the calculus of variations. For this, we will verify the coercivity and the weak sequential lower semicontinuity of the functional ( , ⋅, ⋅). To deduce the latter for the polyconvex functional ( , ⋅, ⋅) we use the following result on the convergence of minors of gradients, which goes back on [53,69], cf. also [52]. With this at hand we now establish weak sequential lower semicontinuity and coercivity. , and Lemma 3.4 it is for < 1: which states (64).

Proof of 2.:
To establish the weak sequential compactness of the energy sublevels we now consider a sequence ( , ) ∈ℕ ⊂ × with ( , , ) ≤ uniformly for all ∈ ℕ. Coercivity estimate (36d) thus allows us to employ Prop. 3.7, which, by Prop. 3.5, implies the existence of a subsequence ⇀ in , such that ⇀ in , cof ∇ ⇀ cof ∇ in 2 ( 0 , ℝ 3 × 3), and det∇ ⇀ det∇ in 3 ( 0 , ℝ). Moreover, thanks to (64) we also find a subsequence ⇀ in . It thus remains to deduce the weak sequential lower semicontinuity of each of the contributions of .

Proof of Prop. 3.10
We start with the proof of the discrete notion of solution (59). Proof of (59a): We observe that a minimizer of problem (33a) equivalently satisfies for all̃∈ Applying the definition of the interpolants (58) we find (59a). Proof of (59b): We use that a minimizer of problem (33b) satisfies the corresponding Euler-Lagrange equations for all̃∈ 0 : Again using the definition of the interpolants (58) we find (59b). Proof of (60): In order to find the energy dissipation estimate we test the minimality of in (33b) by −1 , exploit the minimality of , and add and subtract ( −1 , −1 , −1 ), Consider now ∈ ( −1 , ]. Then summing up the above relation over ∈ {1, … , }, results in Again by the definition of the interpolants and using that ∈ ( −1 , ] we see that this relation is equivalent to (60).

Proof of (61):
We now want to exploit the previously obtained discrete estimate (66) to deduce the apriori estimates (61). To do so, we apply (53) under the integral of (66). This allows us to apply the classical Gronwall inequality and, following the arguments of e.g. [61], one finds for every ∈ {1, … , } This translates into (61a) and, thanks to (64), also yields the estimates (61b)-(61d ) .
Observe that each of the three contributions 2 1 ( ∇ ), 1 ( ∇ ), and 2 ( ∇ ) generates a lower semicontinuous energy term. Therefore the above arguments can be carried out separately for each of the three terms (keep only one of them on the left, the other two then occur with negative sign on the right, alike the Neumann term in (68), and then pass to the limit by lower semicontinuity). In total, this implies that Further note that each of the three contributions 2 1 ( ∇ ), 1 ( ∇ ), and 2 ( ∇ ) also generates a lower semicontinuous energy term if integrated over any measurable subset ⊂  0 . Let now ( ) ⊂ ∞ ( 0 ) such that → in ( 0 ) and ‖ ‖ ∞ ( 0 ) ≤ 1 for all ∈ ℕ. Then, we may argue for = 1, 2 The first inequality of (71) is due to lower semicontinuity of the functional on measurable. To pass from the third to the fourth inequality in (71) we have exploited convergence (70b), resp. (70c), for the first term and the lower semicontinuity of the other two terms on  0 ∖ , resp. measurable; this yields an estimate from above as they here occur with negative sign. Hence we have shown that ∫ ( ∇ ( )) d → ∫ ( ∇ ( )) d for any ⊂  0 measurable. This, together with the uniform bound ‖ ( ∇ ( ))‖ ≤ is equivalent to weak 1 -convergence, cf. e.g. [71,p. 181,Cor. 2.58], and thus finishes the proof of (62f).
Proof of the evolutionary variational inequality (63c) for > 0: Let now > 0. To show the validity of (63c) we observe that (̇) +̃≤ 0 for everỹ∈ 0 with̃≤ 0 a.e. in  0 . Hence, when moving this term to the other side of the equation, we find that (59b) can be reformulated as an inequality, i.e., for all̃∈ 0 with̃≤ 0 a.e. in  0 : We can then pass to the limit on the right-hand side of (73) using convergences (62) and weak-strong convergence arguments and by arguing for the first term as in the case = 0. Proof of nonpositivity (63d) if > 0: From the third bound in (61a), we gather that By weak lower semicontinuity and convergence (62a) we conclude that which implies thaṫ≤ 0 a.e. in (0, ), a.e. in  0 . Proof of the energy-dissipation estimate (63f): Thanks to convergences (62) we can pass to the limit on the lefthand side of the dicrete energy-dissipation estimate (60) by lower semicontinuity arguments, also using in the case > 0 that  (̇( )) ≥  0 (̇( )) and  (̇( )) =  0 (̇( )) sincė( ) ≤ 0 for a.e. ∈ (0, ) by (63d). On the right-hand side, the energy at initial time is constant wrt. ∈ ℕ and we only have to take care about the limit passage in the powers of the external loadings. For this, we want to show that where ( , ( )) ∶= sup{ ( ,̂, ( )),̂∈ argmin( , ⋅, ( ))} is introduced as a surrogate for the partial time derivative from (52). We can conclude (74) if we first show that ( ) is a minimizer of ( , ⋅, ( )) (75) and secondly verify that Clearly, these two properties imply (74) due to the definition of . In addition, by the power control estimate (53), we see that ∫ 0 ( , ( )) d is well-defined and finite. We now prove statement (75). For this, we introduce a further interpolant, i.e.
Using convergences (62e) & (78) and by repeating the arguments of the proof of minimality condition (63a) we conclude (75). We now turn to the proof of the convergence of the powers of the energy (74). For this, we will adapt the arguments of [51,Sec. 3] and [61,Prop. 3.3] to the present, rate-dependent situation. More precisely, for  ∶ [0, ] × × ,( , , ) ∶= is uniformly continuous.

Examples of energy densities satisfying assumptions (36)
In this section we discuss well-known constitutive laws in nonlinear elasticity with regard to their admissibility for assumptions (36), which are at the core of our existence result. The related proofs of the main statements can be found in the Appendix. Note that assumptions (36) are formulated for an energy density in dependence of a matrix = ∇ ∈ ℝ 3×3 and its minors. However, by making use of the assumption of material frame indifference and isotropy, many constitutive laws used in engineering are equivalently formulated with respect to invariants of the right Cauchy-Green strain tensor ∶= ⊤ or with respect to the modified invariants introduced in Section 2, as it will also be the case in the numerical examples shown in Section 4. Then it is not obvious that constitutive laws given in this way also match with the assumptions (36) of our existence theorem. This is why we will now take a closer look at densities given as functions of the invariants 1 ( ), 2 ( ), 3 ( ) and at densities given as functions of the modified invariants ( ), ( ), 3 ( ). In the subsequent discussion we will neglect the dependence of the densities on the phase-field parameter; in other words, the densities studied here play the role of 1 and 2 in (36a). For our further investigations we recall the notation 1 ( ) ∶= tr , 2 ( ) ∶= tr cof , 3 ( ) ∶= det , for ∶= ⊤ , (85a) Here, the expressions in (85a) are the invariants of the right Cauchy-Green strain tensor ∶= ⊤ and , in (85b) denote the modified invariants of . Using the relations the modified invariants ( ⊤ ), ( ⊤ ) can also be reformulated directly in terms of , i.e., As in Section 2 we may also set ∶= (det ) −1∕3 and ∶= (det ) −2∕3 cof and find that | | 2 = ( ) = ( ) and | | 2 = ( ) = ( ).
In accordance with (85), we subsequently assume that we are given densities , • ,˜, and , which satisfy the relation for all ∈ ℝ 3×3 and = ( , cof , det ). In (86), the density˜∶ ℝ × ℝ × ℝ → ℝ is a function of the invariants (85a) of the matrix ⊤ and the density ∶ ℝ × ℝ × ℝ → ℝ is a function of the modified invariants (85b). The first Piola-Kirchhoff stress tensor is determined as with the expressions for the derivatives of the invariants gathered in the next lemma. We point out that (89) provides a stress control for the invariant functions alike (36e), which, in view of (87), will be used lateron to formulate sufficient conditions for the densities˜, in order to guarantee the stress control (36e) for .

Lemma 3.14 (Derivatives of the invariants).
Let the relations (85) hold true. For a matrix ∈ ℝ 3×3 with det > 0 it is Moreover, let be a placeholder for the invariant functions 1 , 2 , 3 , , ∶ GL + (3) → ℝ. The invariant function satisfies the following stress control estimate: The proof consists of a straight forward but lengthy application of the product and chain rule and we carry it out in detail in A.1. There we also give the proof of the next lemma, which provides continuity estimates for the invariants and their stresses in terms of moduli of continuity multiplied by the invariant, as required in the assumption (36f). Lateron, they will be used to deduce similar continuity estimates for the first Piola-Kirchhoff stress.
We now discuss the properties of the constitutive laws (91). We refer to the works [53,56,80], where the polyconvexity of some of (91) and many other constitutive laws, such as, e.g., Ogden's materials, has been discussed. We here collect statements on the polyconvexity of the constitutive laws (91): In view of (85c) we also have the following immediate results regarding the coercivity of the polyconvex constitutive laws (91a)-(91c).
Then the uniform continuity of the stresses (36f) is true. The proof exploits the form of the first Piola-Kirchhoff stress (87) and uses the stress control estimate (92) and the continuity relation (93) for the invariant functions , = 1, 2, 3. For more details we point ahead to Appendix C.5, where the calculations are carried out for depending on the modified invariants.
As a result of Proposition 3. 16

Assumptions (36) and the modified invariants
In the following we discuss the Assumptions (36) for a stored energy density that is a function of the modified invariants introduced in Section 2. Quite often in literature, the density˜is used and assumed to be a function of the modified invariants ( ), ( ), 3 ( ). In this spirit we will now consider the constitutive laws from (91) as functions of ( ), ( ), 3 ( ), i.e.
In [55] the properties of the modified Neo-Hooke and Mooney-Rivlin material (94a) and (94b) have been analyzed. In particular, it is shown that the term ( ⊤ ) = | det | −4∕3 | cof | 2 is not polyconvex itself, so that the modified Mooney-Rivlin material (94b) cannot be polyconvex either. Moreover, in [55] optimal coercivity properties are derived for functions of modified invariants, which guarantee the validity of Ball's existence result [53,Thm. 6.2], see also the discussion in Remark 3.6.
We now discuss the properties of the constitutive laws (94). Firstly, we refer to the works [56,57], which investigate the polyconvexity properties of the modified Arruda-Boyce model (94c) and of a modified Rivlin-Saunders-type model. Amongst others, they also give the following polyconvexity result: Table 3]. Let ∈ ℝ 3×3 . Then the following terms are polyconvex:
In addition, we now gather the following statement on the polyconvexity of further energy terms depending on modified invariants:  (96a) 6. The energy density 3 is polyconvex, where that the condition ≥ + 1 is not only sufficient but even necessary for polyconvexity. More precisely, they show that ≥ 3∕2 in (97a) and ≥ 3 in (97b) are necessary conditions.
Coercivity condition (36d) is formulated for the density as a function of . We now transfer (36d) into an analogous condition for as a function of , , 3 , resp. , , det , cf. (85).
Remark 3.25. If the modified invariants are used, the density depends on ( ⊤ ), hence at least on and on det . Therefore, the sub-cases (36d) b1) & b2) are irrelevant and the coercivity estimate has to feature the term | det | 3 with a suitable power 3 satisfying (98a). An analogous observation holds true if the energy density also depends on ( ⊤ ).
We further observe that, even if enriched by an additional summand involving the determinant, neither the modified Neo-Hooke material (94a) nor the modified Mooney-Rivlin material (94b) complies with the modified coercivity condition (98), since in both cases = 2, which does not allow for > ≥ 2. Yet, according to (98a), it is possible to find > 3∕2 for  (36). Remark 3.30 (Assumptions (36) and the anisotropic split (16), resp. (24)). As explained in Section 2.2 we may apply the anisotropic split (16) to the modified invariants in order to account for the anisotropy of damage. This neither affects the polyconvexity nor the coercivity properties of the constitutive law. As we have seen in (23), the modified invariants ensure the differentiability of the energy density in with ( , ) = 0. Thus, also the results on the stress control and the continuity of the stresses remain valid. If the anisotropic split is applied only to energy contributions in (99) with positive powers (but not to

NUMERICAL EXAMPLES
In this section we explain shortly the main equations for the discretization within the finite element framework and we demonstrate the robustness of the proposed model and the analytical results by a series of numerical examples in a next step.

Discretization
At first, the weak forms of the governing equations and the discretization are summarized. The elastic boundary value problem is based on the balance of linear momentum (5) and the crack phase-field evolution (8). For fixed time, the weak form of the coupled problem reads: Find ∶  0 → ℝ 3 and ∶  0 → [0, 1] such that and The term Ψ in (103) basically serves as a driving force for the phase-field. Moreover, the spaces of admissible test functions 0 and 0 are defined as 0 = { ∈ 1 ( 0 ; ℝ 3 ) | = 0 on  D 0 }, where 1 ( 0 ; ℝ 3 ) denotes the Sobolev space of square integrable functions with values in ℝ 3 and with square integrable weak first derivatives. Correspondingly, the space of admissible test functions for the phase-field equation can be formulated as 0 = 1 0 ( 0 ) ∩ ∞ ( 0 ) = { ∈ 1 ( 0 ) ∩ ∞ ( 0 )| = 0 on  0 }.
To apply the finite element method, the domain  0 is subdivided into a finite set of non-overlapping elements For discretization we use Lagrangian polynomials for both fields. In particular, the ansatz functions for the mechanical field are denoted by and the shape functions for the phase-field bỹ. The valueŝ( ) and̂( ) are the nodal displace-F I G U R E 6 Boundary conditions (left) of a mode-I-tension test and the related mesh based on a hierarchical refinement strategy (right) ments and the nodal values of the phase-field.
Inserting the proposed approximations (105) and (106) into the weak formulations (102) and (103), gives, after a straightforward calculation, the final finite element system.
The time integration is based on an implicit Euler-backward scheme regarding the phase-field parameter , whereby the time interval [0, ] is divided into pairwise disjoint equidistant subintervals with the time step △ ∶= +1 − . At last, the system of equations is solved by making use of a direct solver. In general there exist two solution strategies for the non-linear system (102) and (103), the monolithic and the staggered scheme. In the first the fully-coupled system is solved in each time step. In the staggered scheme the solution is split for the phase-field and for the mechanical field , which means that in each time step both quantities are solved successively. Our analysis in Sect. 3 is based on the latter approach and also for the numerical examples presented in the subsequent exposition we rely on the staggered scheme. In [81] both solution strategies are discussed in more detail.

Mode-II-shear test in two dimensions
As a first numerical example we choose a simple mode-II-shear test in two dimensions and consider a squared plate with horizontal notch. At the lower boundary of the plate the displacements are constrained in horizontal and vertical direction and at the upper side prescribed displacements are applied incrementally in -direction, see Figure 6. Furthermore, the mesh presented in Figure 6 consists of 128 × 128 quadrilateral elements.
The following simulations use the non-linear Yeoh material model (99a) and the proposed anisotropic split (16) so that the strain energy function can be formulated as   (17). The dimensionless material parameter are chosen as 1 = 2.6923 × 10 10 , 2 = 1.3462 × 10 10 and = 2.01923 × 10 11 . Refering to the SI unit system this corresponds to a Young's modulus of = 2.1 × 10 1 1 N m 2 , a Possion's ratio of = 0.3 and the critical energy release rate of  = 2.7 × 10 3 J m 2 . The mobility parameter is = 1 (Pa s) −1 . The related length-scale parameter depends on the element size ℎ and has to fullfill the inequality ≥ 2ℎ in general, cf. [82], which enables the approximation of a diffuse interface zone. In this case using three unniform levels of refinement, the length scale parameter is set to = 2ℎ min = 7.8125 × 10 −6 m. The snapshots of the phase-field and the related crack propagation related to the shear test are shown in Figure 7.
In a next step the influence of the step-size of time is investigated in more detail. Therefore, the time step size is varied such that a bigger and a smaller time step is applied. Within this assumption different time steps are examined and the related load-deflection curves are shown in Figure 8. The results show the convergence to a solution of the time-continuous formulation within decreasing stepsize.

Mode-I-tension test in three dimensions
In a next step we introduce an example in three dimensions. We consider a block with a horizontal notch which consists of 10 × 4 × 10 elements before refinement. The geometry and the related boundary conditions can be found in Figure 9.
All the boundary conditions are realized by using Dirichlet boundary conditions. In Figure 9 also the mesh of the block is shown after applying the hierarchical refinement strategy. Also in this example we use the Yeoh material model (107) with positive coefficients , > 0 for any ∈ {1, 2} and ( ) defined in (17). The material parameter are chosen as 1 = 2.6923 × 10 10 , 2 = 1.3462 × 10 10 and = 2.01923 × 10 11 which are based on the relations given in [79] and correspond to the same material values as above. The length-scale parameter depends on the mesh size and is chosen as = 1.25 × 10 −5 m. Moreover, a time step of Δ = 0.01 sec is applied.  Figure 10 the crack growth during the simulation can be observed. The crack propagates within this loading along the expected crack path.
The load deflection curves for different time step sizes are shown in Figure 11, the block is cracked completely at a vertical displacement of ≈ 0.28 × 10 −3 m.

A C K N O W L E D G E M E N T S
M.T. has been partially supported by the Deutsche Forschungsgemeinschaft in the Priority Program SPP 1748 "Reliable simulation techniques in solid mechanics. Development of non-standard discretization methods, mechanical and mathematical analysis" within project "Finite element approximation of functions of bounded variation and application to models of damage, fracture, and plasticity" -Project Number 255461777 (TH 1935/1-1, TH 1935/1-2). The authors K.W. and C.B. gratefully acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG) under the project "Large-scale simulation of pneumatic and hydraulic fracture with a phase-field approach", also as part of the Priority Program SPP 1748 with the Project Number 255801726. Open access funding enabled and organized by Projekt DEAL.

R E F E R E N C E S
< 0. Thus, is not rank-one convex, hence not polyconvex.