A unified phase-field model of fracture in viscoelastic materials

The phase-field approach has proven to be a powerful tool for the prediction of crack phenomena. When it is applied to inelastic materials, it is crucial to adequately account for the coupling between dissipative mechanisms present in the bulk and fracture. In this contribution, we propose a unified phase-field model for fracture of viscoelastic materials. The formulation is characterized by the pseudo-energy functional which consists of free energy and dissipation due to fracture. The free energy includes a contribution which is related to viscous dissipation that plays an essential role in coupling the phase-field and the viscous internal variables. The governing equations for the phase-field and the viscous internal variables are deduced in a consistent thermodynamic manner from the pseudo-energy functional. The resulting model establishes a two-way coupling between crack phase-field and relaxation mechanisms, i.e. viscous internal variables explicitly enter the evolution of phase-field and vice versa. Depending on the specific choice of the model parameters, it has flexibility in capturing the possible coupled responses, and the approaches of recently published formulations are obtained as limiting cases. By means of a numerical study of monotonically increasing load, creep and relaxation phenomena, rate-dependency of failure in viscoelastic materials is analysed and modelling assumptions of the present formulation are discussed.

for crack propagation into an energy minimisation problem. Bourdin et al. [5,6] introduced a regularisation of the underlying energy functional. The key idea of this approach is a diffuse crack representation, i.e. the crack surface is described by an auxiliary field variable which continuously varies from the intact to the fully broken material state. In general, this auxiliary variable is referred to as phase field. In other words, cracks are no longer seen as sharp discontinuities, but approximated over a finite length scale c . Furthermore, the loss of stiffness that comes along with failure is expressed by means of a degradation function depending on the phase-field. For small c , the regularised functional -converges against the sharp formulation [7], i.e. the fundamental theory of Griffith is recovered. Based upon the diffuse crack representation, phase-field models of brittle fracture have been proposed which include several advancements and a variety of improvements, see [8][9][10] and [11], or [12,13], inter alia. For a detailed review on different formulations for brittle fracture at small strains, the reader is referred to [14].

Phase-field models for inelastic materials
Rate-independent material behaviour. More recently, the phase-field approach is extended to materials which undergo an irreversible deformation before failure. In the literature, different modelling approaches to ductile fracture are introduced. For a detailed comparison of several formulations in the small strain context, see [15].
A first attempt to describe brittle fracture of elasto-plastic solids is made by Duda et al. [16]. Although the inelastic permanent deformations are considered, the coupling between plastic deformation and damage is not taken into account. However, the implementation of such a coupling is crucial if typical crack phenomena which are experimentally observed in ductile materials shall be reproduced, cf. for instance [17].
For this purpose, a non-energetic ductile fracture driving force based on the accumulated plastic strain is introduced by Miehe et al. [18,19]. Ambati et al. [17,20] considered an enhanced degradation function which, in addition to the phase-field variable, depends on plastic deformation. This coupling resulted in a distinct plastic contribution to the fracture driving force. Yin and Kaliske [21] recently adapted a concept which is well-established in the phase-field modelling of fatigue, cf. [22][23][24], to the description of ductile failure. Instead of defining a fracture driving force related to plastic deformation, the fracture toughness is assumed to decrease with increase of inelastic deformations. In all these contributions, plastic deformations play an essential role during the evolution of the phase-field variable. However, the phase-field variable does not explicitly influence the evolution of inelastic strain.
In several other models of ductile fracture, a total energy functional is considered in which both elastically stored energy and plastic work, or plastic energy, are assumed to be affected by fracture. In other words, degradation functions decrease the respective contributions to the energy functional. Consequently, the plastic energy which actually corresponds to the dissipation due to inelastic deformation and possibly includes hardening terms, naturally becomes part of the fracture driving force which is deduced from the functional. Alessi et al. [25,26] introduced this approach in the context of gradient damage models. Related phase-field rate-independent fracture models are proposed in [27,28]. Based on similar considerations, and furthermore, in order to capture the influence of failure on the evolution of inelastic deformation, a yield condition which depends on the phase-field is postulated by Borden et al. [29]. To the best of the authors' knowledge, only models which can be attributed to this group of formulations establish a two-way coupling between plastic deformation and failure.
Rate-dependent material behaviour. In recent years, few efforts have been devoted to describe the ratedependent material fracture behaviour using the phase-field fracture models.
A phase-field fracture model for viscoelastic solids is proposed by Schänzel [19]. Considering a viscoelastic material model and a phase-field driving force based on a generalised principal stress criterion, fracture response in rubber polymers is investigated. Furthermore, the viscosity assumed for the evolution of the phase-field which originally is solely numerically motivated, cf. [9,10], is understood as a material parameter and identified from experimental data.
Similar to the models of fracture in rate-independent materials which are based on energetic considerations, more recently, some formulations adapted to a rate-dependent response of the bulk were proposed. Shen et al. [30] presented a phase-field model for viscoelastic solids in the small strain context. Therein, the evolution equation for the viscous internal variables is assumed a priori and an energy functional which incorporates a stored viscous free energy is postulated. This term corresponds to a portion of accumulated viscous dissipation and is degraded in the same manner as equilibrium and non-equilibrium part of the elastic strain energy density. From the free energy functional, the phase-field equation is deduced by means of a microforce balance.
Loew et al. [31,32] considered a finite linear viscoelastic material model [33] which does not allow to distinguish between the non-equilibrium part of actually stored strain energy and accumulated viscous dissipation. Accordingly, a viscous contribution to the free energy is defined as the sum of the two quantities and assumed, as well as the equilibrium part of elastic strain energy, to be degraded in case of fracture. A similar finite linear viscoelastic material model was investigated by Liu et al. [34]. Yin and Kaliske [35] combined the phase-field approach to fracture with a model of finite viscoelasticity [36]. Unlike [30] or [31,32], in [35], the effectively stored strain energy is assumed to promote crack propagation, i.e. a dissipative contribution is not accounted to the crack driving force. Recently, a similar formulation is adopted by Brighenti et al. [37] based on statistical mechanics-based equations for the response of the bulk material.
While the aforementioned models address a rate-dependent response of the bulk material, Miehe et al. [18] and later Yin et al. [38] studied a strain rate-dependent resistance against fracture combined with a rate-independent model of deformation.

Intent of the present contribution
The aim of the present contribution is to present a unified, thermodynamically consistent phase-field model for fracture of viscoelastic solids. With the central purpose of describing a two-way coupling between damage and viscous effects, we introduce an energy functional which consists of the free energy that includes a contribution which is related to viscous dissipation, and the fracture contribution. From the proposed functional, the governing equations of the phase-field and viscous internal variables are derived in a thermodynamic consistent manner. Depending on the specific choice of the degradation functions and model parameters, respectively, the modelling approaches of [30][31][32] or [35,37] can be retained as limiting cases of the present model. With this unified formulation at hand, we study the coupling between viscous effects and fracture. This paper proceeds as follows. In Sect. 2, the proposed phase-field model of fracture in viscoelastic material and its thermodynamic consistency proof are illustrated. Furthermore, algorithmic aspects are addressed. In Sect. 3, the analysis of coupling between viscoelastic material response and fracture by means of numerical examples is presented. The paper concludes with a short summary and an outlook regarding the future work.
Within this paper, the summation convention applies for tensor quantities. In index notation, partial derivatives are denoted by (•) , j = ∂•/ ∂x j .

Derivation of the unified phase-field model
In this section, the unified phase-field description of fracture in viscoelastic materials is presented. A general energetic formulation of fracture adapted to viscoelastic materials is presented first. Subsequently, the specific constitutive assumptions for the bulk are specified and thermodynamic consistency is proven. Finally, governing equations and algorithmic aspects are provided. Without loss of generality, for the present contribution, we limit ourselves to the small strain setting and neglect inertia effects.

The energetic approach to fracture
The hypothesis of Griffith [4] defines the starting point for the phase-field approach to fracture. Dissipation due to crack growth is referred to as surface energy which is assumed to increase proportionally to newly created crack surface. In a simplified setting, crack growth is assumed to take place if the strain energy stored in the material equals the amount of surface energy corresponding to the increase in fracture surface. Following [3], this fundamental theorem can be reformulated by means of an energy functional wherein ψ and G c denote the density of free energy stored in the material and the fracture toughness, respectively. The two constituents to the functional Π correspond to stored free energy Π sd and dissipation due to formation of crack surface or surface energy Π fr . The domain under consideration is denoted by Ω ⊂ R N and Γ ⊂ R N −1 is the crack surface. For sake of brevity, terms arising from external loads are omitted. In order to make this energetic approach to fracture accessible to numerical implementation, a regularisation of the functional Π is introduced [5]. For this purpose, a phase-field variable which continuously varies from the intact (d = 0) to the fully broken (d = 1) material state is defined. Using this variable, the crack surface is approximated by means of the functional with the crack surface density 1 cf. [9], which incorporates the regularisation parameter c . This parameter defines the characteristic width of the crack approximated by the phase-field variable d. A part from its geometrical role, it also influences the critical stress and strain and can be seen as material parameter, cf. [40]. Figure 1 illustrates the concept of diffuse crack approximation in comparison to the sharp representation. With the diffuse approximation of crack surface at hand, dissipation due to crack evolution can be expressed as wherein Φ is defined as density of fracture dissipation. In order to account for the decrease of free energy due to evolution of cracks, the degradation function is defined which has to fulfil the conditions Based on the degraded free energy density the regularised functional of free energy is given by and the regularised total energy, or pseudo-energy, functional reads 2 In case of linear elasticity, it can be discussed that Π c γ -converges to Π for c → 0, see e.g. [41]. For γ c according to (4) and a specific choice of the degradation function g(d), an analytical proof of γ -convergence can be provided, cf. [7].
Generalisation for inelastic material response. In case of an non-elastic response of the bulk material, in the total energy functional Π c , an expression for the free energy density Ψ can be considered which adequately accounts for coupled inelastic deformations. For the description of combined fracture and viscous mechanisms, the proposed free energy density functional consists of two essential constituents: Naturally, the first one is the effectively stored strain energy which, in case of damage, is decreased by means of the degradation function g st (d). In addition, a free energy contribution related to accumulated viscous dissipation is assumed. The definition of ψ st and ψ vi can take any form suitable for describing viscoelastic deformations. The specific definition considered in this work is given in Sect. 2.2. In order to keep the formulation as general as possible, similar to [29] and [30], the parameter β vi ∈ [0, 1] is introduced as a weight of the contribution related to inelastic deformation. For the viscous term, the degradation function g vi (d) is prescribed. In contrast to e.g. [30], coincidence of g st (d) and g vi (d) is not a priori supposed. Therefore, the present formulation enables to account for an influence of damage on energy storage and dissipative viscous mechanisms, and vice versa, which may be not completely identical. The free energy density contribution Ψ vi in (11), related to viscous dissipation, can be interpreted as follows. Firstly, mechanical energy can be released due to the relaxation mechanisms present in the viscoelastic material, where an amount of strain energy which has been set free can convert into heat, i.e. lead to an increase in temperature. As experimental evidence has shown that resistance against failure of many rate-dependent materials decreases with temperature, see e.g. [42,43], viscous dissipation can be assumed to promote fracture. By means of proposed Ψ vi , the formulation enables to capture this effect, although the thermomechanical coupling is not explicitly modelled. Secondly, micro-mechanical mechanisms which lead to viscous dissipation can directly diminish the amount of strain energy that has to be stored in the material so that crack propagation can take place. Consider, for example, rubbery polymers under tensile loading. Viscous dissipation can be attributed to relaxation, i.e. reorientation and disentanglement of chain segments. These mechanisms go along with a redistribution of loading from entangled to crosslinked polymer chains, as discussed in [19]. As a result, in case of fracture, crosslinked chains are broken, essentially. The amount of strain energy necessary for crack propagation thus is lowered compared to the case where no viscous dissipation has taken place, cf. [42] for experimental results.
Choice of the degradation functions. In the literature, the quadratic degradation function in which a small residual k is included in order to enhance numerical stability, is adopted frequently, cf. for instance [5,9,27,30]. Nevertheless, several alternatives were proposed. In order to obtain a more brittle response, quartic and cubic expressions were introduced by [44] and [45]. With the same objective, [46] suggested an exponential function and [35,38] studied a sinusoidal expression. The functions of [45] and [46] as well as the expressions proposed by [47,48] and [49] are parametric, i.e. they include additional parameters which can be fitted to the physical behaviour of the material under consideration. The present formulation for the free energy allows to apply this concept for both the actually energy storage and dissipative mechanisms by means of g st (d) and g vi (d), separately.
Within the scope of this contribution, we focus on the development of the phase-field formulation rather than modelling the response of a particular material. To that end, herein, the quadratic function defined in (12) is considered for g st (d). In Sect. 3, regarding g vi (d), we study two limiting cases. Firstly, g vi (d) = g st (d) is Secondly, assuming that the nature of relaxation mechanisms remains unaffected even if damage compromises the stiffness of the material, g vi = 1 = const. is set. 3

Constitutive modelling of the bulk response
For the viscoelastic behaviour of the bulk material, a generalised Maxwell model, as depicted in Fig. 2, is adopted. Kinematic assumptions. For the response of viscoelastic materials, often significant differences between volumetric and isochoric deformation are observed, cf. [50]. Accordingly, the strain tensor is additively decomposed into the volumetric and isochoric parts Furthermore, for each non-equilibrium branch m ∈ [1, n] of the generalised Maxwell element, the strain is split into a viscous portion m α i j = vol m α i j + m a i j and an elastic fraction. Descriptively, vol m α i j and m a i j can be interpreted as the strain of the dampers depicted in Fig. 2, while vol ε i j − vol m α i j and e i j − m a i j correspond to the strain of the spring in the respective non-equilibrium branch.
Elastic strain energy density. The strain energy density Ψ st which is stored in the material is additively decomposed as where the equilibrium part is ψ st,eq ε i j = vol ψ st,eq ε pp + iso ψ st,eq e i j (17) and the non-equilibrium or over-stress contribution reads 3 For a more detailed discussion of the assumption g vi = 1 = const. from a physical point of view, see Sect. 3.1.
with either contribution additively decomposed into a volumetric and isochoric fraction. 4 For both the equilibrium and non-equilibrium part, linear elasticity is assumed. Accordingly, the respective terms are given by vol with K eq , m K ov > 0 and μ eq , m μ ov > 0 denoting the corresponding compression and shear moduli, respectively. 5 Note that in the above definitions, the summation convention does not apply for the index m.
Free energy contribution related to viscous dissipation. Similar to Ψ st , the dissipative contribution is decomposed into the volumetric and isochoric part Assuming constant viscosities vol The viscous term ψ vi can be interpreted in analogous manner to the virtually undamaged strain energy density ψ st : Descriptively, it can be understood as the viscous dissipation which would have been accumulated over time in the absence of damage.

Thermodynamic consistency
In order to prove that the proposed model is conform to the second law of thermodynamics, the density of dissipation powerḋ is considered which corresponds to the difference between external stress power and rate of free energy. The Clausius-Duhem inequalityḋ ≥ 0 ( 2 4 ) demands a locally non-negative dissipation power. Inserting the above definition of Ψ , Eq. (11), this constraint expands to whereinḋ fr and volḋ vi , isoḋ vi correspond the density of dissipation power due to crack evolution and viscous mechanisms, respectively. Standard arguments [52,53] lead to the definition of stress and the residual inequalitiesḋ As the elastic constants as well as β vi are defined to be positive, and both ψ st and ψ vi are quadratic, together with (7) 3 , holds, i.e. the thermodynamic conjugate to the phase-field variable is non-negative. Accordingly, constraint (27) reduces toḋ i.e. irreversibility of damage. In Sect. 2.5, the fulfilment of this constraint will be addressed.
Evolution of viscous strain. With (28) and (29), the evolution equations for the viscous internal variables are obtained from Note that, with 0 m s ov i.e. evolution of viscous strain is essentially related to the over-stress in the corresponding non-equilibrium branch. In the absence of damage, i.e. g st (d = 0) = g vi (d = 0) = 1, the equations resemble the classical linear-elastic assumption ηα = σ ov . Inserting the above definitions, the residual inequalities (28) and (29) read which is fulfilled for arbitrary (m)α i j due to (6) Remark on the definition of Ψ vi . With the above constitutive equations at hand, a short remark on coherency between the viscous contribution to free energy density Ψ vi and the time integral of stress power D vi in the damper elements of the non-equilibrium branches of the generalised Maxwell model is given, cf. Fig. 2. Classically, the damper elements are assumed to account for inelastic mechanisms and D vi is referred to as density of accumulated viscous dissipation. Considering, exemplary, the volumetric response of a single non-equilibrium branch m on time t, the former reads vol while the latter is given by where vol m σ ov i j = 0 m σ ov pp /3 δ i j denotes the volumetric over stress contribution. Note that for Ψ vi (t), degradation upon time t is assumed, while the permanently degraded over stress enters accumulated viscous dissipation vol m D vi . Applying the mean value theorem, vol m D vi can be rewritten as vol witht ∈ [0, t]. From (7) and (31), it follows which leads, together with β vi ∈ [0, 1], to vol Descriptively, Ψ vi recovers a certain fraction of accumulated viscous dissipation. The fact that for Ψ vi (t), in contrast to vol m D vi , degradation is prescribed not before current time t, does not lead to an overestimation of the respective fraction of vol m D vi .

Governing equations
The evolution of the phase-field variable is deduced from the total energy functional Π c by means of the variational derivative where η f is introduced as a kinetic fracture parameter for purely numerical purposes, i.e. for enhancing the stability of the solution scheme, cf. [27,54]. The constant η f is chosen such small that its influence on the simulation results vanishes which is verified by means of a comparative study of different values. 6 Inserting the definition of Ψ (11) into Π c , (42) 1 takes the form from which it becomes clear that, depending on the specific choice of β vi and g vi , fracture is driven by stored strain energy and the free energy contribution related to a portion of accumulated viscous dissipation. Note that for g vi (d) = g st (d), the phase-field equation resembles the formulation of [30], while the assumptions on the fracture driving force of [31] and [35] are preserved for g vi (d) = g st (d) and β vi = 1 or β vi = 0, respectively. The governing differential equations of the proposed model are summarised in Table 1.

Algorithmic aspects
Weak forms of the governing equations. For the derivation of the weak forms of balance of linear momentum and phase-field equation, the test functions δ u j ∈ W u and δ c ∈ W c are defined with Therein, H 1 is the Sobolev space of functions possessing square integrable derivatives and Ω u j denotes the part of the boundary where the displacementū j is prescribed, cf. [55]. Then, (a) and (b) from Table 1 are multiplied by δ u j and δ c, respectively. Integration by parts and making use of the divergence theorem yields wheret j denotes the traction vector given on ∂ Ω t j with ∂ Ω t j ∪ ∂ Ω u j = ∂ Ω and ∂ Ω t j ∩ ∂ Ω u j = ∅, and n k is the outward-pointing unit normal vector on ∂Ω. Time discrete forms of (46) and the evolution equations (c) and (d) from Table 1 are obtained by approximating the respective rates using an Euler backward scheme. For spatial discretization, Galerkin's method is applied. Then, the discretized equations are implemented into a standard finite element framework.
Irreversibility of fracture. In order to guarantee irreversibility of fracture, we follow the approach of [10,56] and apply Dirichlet boundary conditions to the phase-field on all nodes where the phase-field variable has reached a critical value d crit : Alternatively, in line with [9], one may introduce a history variable which comprises the maximum of virtually undamaged fracture driving force which has occurred and consider a phase-field equation which is altered accordingly. In contrast to the former approach for crack irreversibility, the latter enables to fulfil constraint (31) exactly, see [9]. However, it has the important disadvantage that it leads to an overestimation of crack surface, cf. [41]. Therefore, we adopt (48) and set d crit = 0.97 unless stated differently.
Staggered solution strategy. Due to its robustness, the staggered solution strategy widely spread in the context of phase-field modelling, cf. [9,14], is adapted. For every increment, initially, the phase-field variable is frozen and the update of the displacement field is determined from the balance of linear momentum. Subsequently, the evolution of phase-field is evaluated assuming constant displacements. In order to avoid inaccuracy of the numerically determined solution due the decoupled approach, multiple staggered iterations are performed at each increment and a convergence criterion for the staggered scheme is established based on the maximum change of the nodal degrees of freedom. The time step is adaptively controlled by means of a heuristic method.

Numerical examples
In this section, numerical examples are presented in order to analyse the fundamental principles of the present model and to discuss the consequences of different modelling assumptions. Firstly, the local response of the material is studied with a focus on the coupling between viscous mechanisms and damage. Subsequently, crack propagation within a domain where stress localisation is triggered by a preexisting notch is considered. While the presented examples can serve for qualitative verification and a better understanding of the coupling between fracture and rate-dependent response of the bulk material, the comparison with experimental results goes beyond the scope of this publication. For all the numerical simulations, a two-dimensional setting is considered. Plane-strain conditions are assumed. Unless stated differently, the kinetic fracture parameter η f is set to zero, i.e. no viscous regularisation of crack growth is considered.

Local response of the material under monotonic loading
In order to study the local material response, a homogeneous domain of unit size under tension is investigated. Displacement boundary conditions are applied as depicted in Fig. 3 which leads to spatially constant  In order to facilitate discussion, a generalised Maxwell model with only one non-equilibrium branch is considered. Identical relaxation times vol 1 T = iso 1 T = 1 s are assumed and K eq = 1 K ov = 175kN/mm 2 , μ eq = 1 μ ov = 80.77kN/mm 2 , G c = 2.7N/mm, c = 0.5mm, k = 10 −5 are defined. Various strain rateṡ ε ∈ [10 −8 , 10] s −1 are investigated. In what follows, the scalars σ , ε and α denote the normal components of the respective tensor quantities in the direction of external load.

Influence of viscous dissipation on damage.
With the goal to investigate the influence of the viscous contribution to the fracture driving force, we set g vi = g st and compare the two limiting cases β vi ∈ {0 ; 1}. For a representative selection of strain ratesε, the corresponding stress-strain curves are depicted in Fig. 4.
Regardless of the choice of β vi , the overall rate-dependency of the material response coincides qualitatively. As in classical viscoelasticity, for very high strain rates, the stress-strain curves converge against the upper elastic limit corresponding to the instantaneous moduli respectively. For lowε, the material response likewise approaches the lower elastic limit where the contribution of the non-equilibrium branches to the resistance of the material vanishes, i.e.
K |ε →0 = K eq and μ|ε →0 = μ eq . It becomes clear from Fig. 5, almost no viscous dissipation takes places until failure for very large as well as for very low strain rates. Accordingly, for both β vi = 0 and β vi = 1, there is no notable viscous contribution Ψ vi to the free energy. Independent of β vi , in the limiting cases, the respective stress-strain curves thus coincide. For intermediate strain rates, the critical stress continuously increases withε, whereas the critical value of strain decreases. This phenomenon arises from the fact that the effective stiffness of the viscoelastic material under consideration is the higher, the larger the rate of strain. Accordingly, the amount of energy stored in the material at a certain level of deformation is more significant the higherε. In contrast, the maximum amount of free energy density which can be stored remains constant, cf. Fig. 6.
When there is a notable amount of viscous dissipation, i.e. at intermediate strain rates, the viscous contribution to the crack driving force leads to lower critical values of both stress and strain compared to the case where only elastically stored energy is assumed to promote fracture. For the material parameters under consideration, the influence of β vi is not prominent. Regarding the critical stress, the relative difference between the respective values for β vi = 0 and β vi = 1 does not exceed 20 % foranyε. As the parameter β vi does not explicitly enter the evolution equations of the viscous internal variables, the dependency of α on ε is unaffected by the choice of β vi . Exemplary, an α-ε curve is depicted in Fig. 8a for an intermediate strain rate and β vi = 0.
Comparison of assumptions for the degradation function g vi . For the purpose of studying the influence of different assumptions for the degradation of the viscous contribution to the free energy, the two limiting cases g vi = g st with β vi = 0 and g vi = 1 = const. are compared. Figure 7 depicts the respective stress-strain curves. Regarding the overall rate-dependency of the material response, the above statements remain valid. For intermediateε, compared to g vi = g st with β vi = 0, the choice of g vi = 1 leads to slightly higher values of critical stress. This effect can be understood by analysing the evolution of the viscous internal variable α. To this end, α-ε and d-ε curves are depicted in Fig. 8 for either g vi and a representative, intermediateε. For g vi = g st , the phase-field eliminates from the evolution of viscous strain. Accordingly, the α-ε curves resemble the shape which would have been obtained in classical linear viscoelasticity. The assumption g vi = g st thus can be seen as analogy to the hypothesis of strain equivalence which can be made within the framework of classical continuum damage models, cf. [57]. In contrast, for g vi = 1, the increase in viscous strain α slows down with growth of the crack phase-field variable. In what follows, this effect is discussed by means of the evolution equation for the viscous strain. Exemplary, the volumetric response is considered. For g vi = 1, (c) in Table 1 takes the form It is evident that damage leads to an increase of the effective relaxation time which can be defined as vol m T /g st (d). From a physical point of view, g vi = 1 corresponds to the assumption that the nature of the viscous mechanisms remains unaffected by fracture. For g vi = 1, the activation energy of the underlying viscous relaxation processes is not supposed to decrease with damage whereas the stiffness of the material or the elastically stored energy degrade. Therefore, in this case, the characteristic time scale associated with the viscous or relaxation mechanisms is assumed to increase with damage.
In particular, when compared to the overall rate-dependency, the influence of the choice of g vi on the response of the material has revealed to be rather small. Nevertheless, from a physical point of view, an a priori assumption of identical degradation functions g vi = g st may not always be straightforward. A thorough analysis which takes the specific behaviour of different materials into account can give further insight. This, however, goes beyond the scope of the present study.

Local response of the material under temporally constant strain or stress
For viscoelastic materials, in addition to monotonically increasing loading, it is of considerable interest to study a temporally constant load as creep phenomena can occur or stress may reduce due to relaxation. Therefore, in the following, the evolution of the local response of the material when a certain level of stress or strain, respectively, is held constant is investigated. In order to account for the rate-dependency of the response, different strain or stress rates are considered during loading stage.
Creep phenomena. In order to study creep phenomena, different stress ratesσ are applied until a certain level of stress σ ref is reached which is then held constant, see Fig. 9a. Various values of maximum stress σ ref ∈ [0.2, 0.3] kN/mm 2 and stress ratesσ ∈ 10 −5 , 10 kN/mm 2 /s have been investigated. As a representative example, the response for σ ref = 0.25kN/mm 2 is depicted in Fig. 10 for g vi = g st , and β vi ∈ {0, 1}. In order to enhance the stability of the solution scheme when it comes to creep failure, η f = 10 −9 kNs/mm 2 is set. From Fig. 10, it becomes clear that for σ ref = 0.25kN/mm 2 , in the absence of a viscous contribution to the fracture driving force, i.e. for β vi = 0, the response does qualitatively coincide for all stress ratesσ . Regardless of the rate of stress during the loading stage (t < 0), with creep time t, strain converges towards an identical, finite equilibrium level. Generally, the strain at the end of the loading stage, ε(t = 0), increases witḣ σ . However, very low ratesσ 10 −3 kN/mm 2 /s correspond to almost quasi-static application of load, i.e. the lower purely elastic limit case, while forσ 10 −3 kN/mm 2 /s, during the loading stage, the upper purely elastic limit is reached for which the increase in strain due to creep, ε(t → ∞) − ε(t = 0), becomes maximal. Increase in strain due to creep goes along with a growth of strain energy density Ψ st and the phase-field variable In that case, the time which elapses until it comes to failure increases withσ as the critical amount of strain energy is reached the earlier the smallerσ , cf. Fig. 10c.
In contrast to β vi = 0, when a viscous contribution to the fracture driving force is assumed, e.g. β vi = 1, both σ ref andσ influence if it comes to creep failure, as Ψ vi decisively depends onσ . In particular, for an identical maximum stress σ ref prescribed, it is possible to predict failure for higherσ , only: The material response is identical for any β vi in case of the lower elastic limit case, i.e. forσ → 0, and β vi > 0 only leads to a slender increase in the equilibrium strain with respect to β vi for small, finiteσ , since there is almost no viscous dissipation for small rates. In contrast, for higherσ , there are more considerable values of Ψ vi . Accordingly, in the latter case, Ψ vi can essentially trigger creep failure, e.g. forσ ≥ 0.05kN/mm 2 /s in Fig. 10. As Ψ vi is the greater the higherσ , in contrast to β vi = 0, for β vi = 1, the time which elapses until failure is the smaller the higherσ . If a critical value of free energy Ψ is reached, it suddenly comes to failure. From Fig. 10c, it can be seen that the critical amount of Ψ is identical for anyσ and coincides with the value obtained for monotonically increasing load, cf. Fig. 6. It is obvious that there are smaller σ ref which are not depicted here, where failure can not be observed for any finiteσ or assumption on the coupling since they do not go along with a sufficiently large amount of Ψ .
Unlike the assumption on β vi , the choice of the degradation functions does only quantitatively influence the material behaviour which can be predicted. If, for example, in case of β vi = 0, viscosity is not assumed to be degraded, i.e. g vi = 1, this can lead to an increase of the time which elapses until it comes to creep failure with respect to g vi = g st , since for g vi = 1, the effective relaxation time grows with the phase-field variable, cf. Eq. (49).

Stress relaxation.
With the goal to demonstrate the capability of the present model to adequately account for relaxation phenomena, different strain ratesε are applied until a certain level of strain ε ref is reached which is then held constant, see Fig. 9b. Various values of maximum strain ε ref ∈ [0.001, 0.004] and rateṡ ε ∈ 10 −8 , 10 s −1 have been investigated. As a representative example, the response for ε ref = 0.002 is depicted in Fig. 11 for g vi = g st , and β vi ∈ {0, 1}. In order to exactly fulfil constraint (31) within the scope of this analysis of the local response, the history variable 7 is considered instead of irreversibility constraint (48). 7 From the definition of vol ψ vi and iso ψ vi in Eq. (22), it directly follows that ψ vi is monotonically increasing. However, an increase in damage could be predicted for decreasing stress if only ψ st would be incorporated into H. Independent of the strain rate or the assumption on the fracture driving force, when strain is held constant, stress relaxes towards a certain equilibrium level while the fracture phase-field remains constant. Since the amount of damage which has been reached at t = 0, i.e. when relaxation begins, depends on bothε and β vi , the equilibrium stress σ | t→∞ changes accordingly. As there is almost no viscous dissipation while ε increases for very high or very smallε and free energy Ψ does not increase when stress relaxes, σ | t→∞ becomes independent of β vi forε → 0 andε → ∞. For g vi = 1, qualitatively identical results are obtained. Compared to the reference g vi = g st for β vi = 0, if g vi = 1 is adopted, a slight increase in effective relaxation time with phase-field is stated, cf. Eq. (49). 8

Heterogeneous case: SENT benchmark
As a representative example for the response of non-homogeneous structure, we follow [8,28,31,38], inter alia, and consider the single edge notched tension test (SENT), i.e. a notched specimen under tensile loading. For the geometry and the displacement boundary conditions which are prescribed, see Fig. 12a. In contrast to the previous section, we adopt the material parameters E eq = 0.327 N/mm 2  with K = E/(3 (1 − 2 ν)) and μ = E/(2 (1 + ν)) from [58] and set G c = 10 −3 N/mm, c = 0.04mm, k = 10 −4 . For optimal computational efficiency, Isogeometric analysis combined with adaptive refinement is employed, see [59], guaranteeing a characteristic element size h < c /4 in the vicinity of the crack. For β vi and g vi , the same limiting cases are investigated as in the previous section.
Monotonically increasing load. Similar to the local response, a monotonically increasing external displacement is investigated first. A large variety of displacement ratesu ∈ [10 −8 , 10] mm/s is considered. The same considerations observed in the previous example apply here as well. Independent of the choice of β vi or g vi , qualitatively coinciding simulation results are obtained. In each case, straight crack growth from the tip of the notch towards the upper edge of the domain is predicted when a certain amount of load is reached, see Fig. 12b. The corresponding force-displacement curves are depicted in Figs. 13 and 14, respectively.
For high displacement ratesu, all the curves converge against the upper elastic limit. For very lowu, the response of the structure likewise approaches the lower limiting case where the non-equilibrium branches do not contribute to the resistance. However, crack propagation leads to finite strain rates even foru → 0. Therefore, in contrast to the stress-strain curves depicted in Figs. 4 and 7, a slight difference is stated between the F-curves for different β vi and g vi , respectively, even if the rate of the external load approaches zero. For intermediate displacement rates, the critical force level continuously grows withu, whereas the displacement at failure decreases. As it has been outlined in the previous section, this effect is explained as follows. The effective stiffness of the viscoelastic material under consideration increases with the rate of strain. In contrast, the amount of fracture driving force necessary for crack growth remains constant.
The significant increase in displacement at failure with decreasing rate of external load can be interpreted as a kind of brittle-to-ductile transition: For the lower elastic limit case, the maximum displacement is more than twice as high as for the upper.
Similar to its influence on the local response of the material, for g vi = g st , the viscous contribution to the crack driving force, i.e. the assumption of β vi = 1, leads to a lower level of both critical force and displacement compared to β vi = 0, see Fig. 13. In particular for intermediateu, its influence is clearly noticeable. In contrast, the definition of g vi = 1, compared with g vi = g st and β vi = 0, leads to a slightly higher critical force level, see Fig. 14. However, for the specific material parameters under consideration, the differences between the response of the notched structure predicted for various modelling assumptions reveal to be not too pronounced. Depending on the maximum force F ref prescribed, the simulation results can be divided into three groups. Firstly, for low F ref , there is an increase in displacement with creep time towards an identical equilibrium value for allḞ and assumptions on the coupling. This increase does not go along with crack propagation. In that case, almost no difference is stated between the results for g vi = g st with β vi = 0 and g vi = 1. Secondly, for intermediate F ref , this phenomenon is observed in absence of a dissipative fracture driving force contribution, only. In contrast, β vi > 0 can trigger creep failure for sufficiently largeḞ, see Fig. 15a, b. Then, in terms of crack path, the result is similar to the case of strictly monotonically increasing external load. Thirdly, for high F ref , failure is observed for any assumption on the coupling between viscous mechanisms and fracture, see Fig. 15c, d. For smallḞ, in this case, failure can also occur while force is still increasing (t < 0), whereas for higherḞ, crack propagation does not start before a sufficiently large amount of creep time t > 0 has passed. Viscous fracture driving force then can lead to an earlier occurrence of failure. In contrast, the assumption of a non-degraded viscosity can effect a notably larger amount of time which has to elapse until it comes to crack propagation compared to g vi = g st for β vi = 0.
The present numerical study demonstrates that in case of creep loads, varied assumptions on the coupling between viscous mechanisms and fracture can lead to significant differences between the predicted responses. From the comparison to experimental data, it remains to be verified which choice for the fracture driving or the degradation functions, respectively, may be most appropriate.

Conclusion and outlook
A unified phase-field model of fracture in viscoelastic materials is presented which enables a two-way coupling between viscous mechanisms and crack growth. For this purpose, a pseudo-energy functional is introduced which consists of the free energy and the dissipation due to crack evolution, and a contribution to the free energy that is related to viscous dissipation is defined. From the energy functional, the differential equations governing the response of the material are deduced in a thermodynamically consistent manner. The approaches of recently published models [30,31,35] are obtained as limiting cases.
By means of numerical examples, the consequences of different modelling assumptions on the local material response as well as on failure of a notched domain are discussed. Especially in case of strictly monotonically increasing external load, predictions are qualitatively similar and the overall rate dependency of failure is captured for all model variants, whereas quantitative differences arise, for instance regarding the critical load level. In contrast, in case of creep, the simulation result can decisively depend on the assumption made for the fracture driving force: A viscous contribution can essentially trigger failure which is not predicted for the same load level if only strain energy is assumed to drive fracture. Furthermore, in the latter case, occurrence of failure does only depend on the maximum value of force, whereas in the former, a possible influence of the rate during the application of load can also be taken into account.
Further work will focus on the comparison of model predictions and experimental data. In addition, an extension to the framework of finite viscoelasticity, cf. [36], is subject of current research. In that case, larger derivations away from the thermodynamic equilibrium path can be captured, possibly leading to a more pronounced influence of coupling between viscous dissipation and fracture.