Phase-field modelling and analysis of rate-dependent fracture phenomena at finite deformation

Fracture of materials with rate-dependent mechanical behaviour, e.g. polymers, is a highly complex process. For an adequate modelling, the coupling between rate-dependent stiffness, dissipative mechanisms present in the bulk material and crack driving force has to be accounted for in an appropriate manner. In addition, the resistance against crack propagation can depend on rate of deformation. In this contribution, an energetic phase-field model of rate-dependent fracture at finite deformation is presented. For the deformation of the bulk material, a formulation of finite viscoelasticity is adopted with strain energy densities of Ogden type assumed. The unified formulation allows to study different expressions for the fracture driving force. Furthermore, a possibly rate-dependent toughness is incorporated. The model is calibrated using experimental results from the literature for an elastomer and predictions are qualitatively and quantitatively validated against experimental data. Predictive capabilities of the model are studied for monotonic loads as well as creep fracture. Symmetrical and asymmetrical crack patterns are discussed and the influence of a dissipative fracture driving force contribution is analysed. It is shown that, different from ductile fracture of metals, such a driving force is not required for an adequate simulation of experimentally observable crack paths and is not favourable for the description of failure in viscoelastic rubbery polymers. Furthermore, the influence of a rate-dependent toughness is discussed by means of a numerical study. From a phenomenological point of view, it is demonstrated that rate-dependency of resistance against crack propagation can be an essential ingredient for the model when specific effects such as rate-dependent brittle-to-ductile transitions shall be described.


Introduction
The mechanical behaviour of many engineering materials depends on rate of deformation.For example, the response of polymers can be much more stiff or brittle when the loading rate is increased, see [1,2].The same applies for natural materials such as cheese [3] or confections [4].In order to reduce experimental effort for design and testing of engineering products as well as for the optimisation of production processes of foods, the computational modelling and simulation of crack phenomena in rate-dependent materials is of increasing interest.
For the modelling of crack phenomena, the phase-field approach to fracture has become a well-established concept.Different from classical finite element approaches (FE), it enables to simulate crack growth without the need for remeshing.Furthermore, complex crack patterns which are not a priori known can be simulated in a straightforward manner, which especially makes the concept attractive compared to alternative approaches such as Cohesive zone elements [5] or the Extended-finite-element-method (X-FEM) [6].The phase-field fracture approach goes back to the variational formulation of brittle fracture of Francfort and Marigo [7], who recast the Griffith criterion [8] for crack propagation into a variational setting.Bourdin et al. [9,10] introduced a diffuse crack representation by means of the phase-field variable, which continuously varies from the intact to the fully broken material state.In other words, cracks are no longer seen as sharp discontinuities, but approximated over a finite length scale ℓ c .Making use of this smeared crack representation, a regularisation of the pseudo-energy functional is carried out.Based upon the fundamental work of Bourdin, numerous phase-field models of brittle fracture have been proposed which include several advancements within both the infinitesimal strain regime [11,12,13,14,15] as well as finite deformation [16,17,18].Further extensions have been proposed to also include fatigue effects, see [19,20,21,22], inter alia.Very recently, the variational approach to fracture also is combined with machine learning and data-driven approaches [23,24,25,26,27].
Furthermore, fracture phase-field modelling has been advanced towards elasto-plastic materials, see [28] for an overview on several approaches within the infinitesimal strain setting.For the performance of ductile fracture models, the description of interaction between inelastic dissipative mechanisms and crack growth has revealed crucial.In particular, in the absence of an adequate coupling, crack patterns that are experimentally observed in metals, for instance, can not be reproduced, see e.g.[29,Fig. 14].Different manners of introducing such a coupling are proposed, including non-energetic ductile fracture driving forces based on accumulated plastic strain [30,31], and an enhanced degradation function which, in addition to the phase-field variable, depends on plastic deformation and results in a distinct plastic contribution to the fracture driving force [29,32].Furthermore, instead of a fracture driving force related to inelastic mechanisms, degradation of fracture toughness depending on equivalent plastic strain is introduced [33].Several other phase-field models of ductile fracture are based upon a pseudo-energy functional in which both elastically stored energy and a plastic quantity, which is referred to as plastic work or plastic energy, are assumed to degrade upon fracture.Depending on the specific formulation, the plastic contribution to free energy actually corresponds to hardening terms [34,35] or accumulated plastic deformation [36].
More recently, the approach is combined with rate-dependent models for the deformation of the bulk material.A first phasefield fracture model for viscoelastic solids is proposed by Schänzel [31], where a non-energetic fracture driving force based on a generalised principal stress criterion is adopted.Alternative driving forces based on energetic or thermodynamic arguments are introduced by Shen et al. [37] as well as Liu et al. [38] within the kinematically linear regime and by Loew et al. [39,40] within the linear viscoelasticity framework [41] at finite deformation.In these models, a viscous dissipative contribution is incorporated into the degraded free energy and thus enters fracture driving force.Different from the aforementioned models, only equilibrium and over-stress parts of the strain energy density are assumed to promote crack propagation by Yin and Kaliske [42], who combined the phase-field approach to fracture with a model of finite viscoelasticity [43].Recently, similar formulations are adopted by Brighenti et al. [44] based on statistical mechanics-based equations for the response of the bulk material, and in [45] where the rate-and temperature-dependent behaviour of polymer nanocomposites is investigated.In some of these models based on either nonenergetic or energetically motivated driving forces, [31,39], the viscosity assumed for the evolution of phase-field which originally is solely numerically motivated, cf.[12,13], is understood as a material parameter and identified from experimental data.In the recent work of Dammaß et al. [46,47], a unified energetic phase-field model for fracture of viscoelastic solids has been presented in the kinematically linear regime.Depending on the specific choice of the degradation functions and model parameters, respectively, the modelling approaches of [37], [39,40] or [42,44] are retained as limiting cases of the present model and by means of representative numerical studies, the coupling between viscous effects and fracture is analysed.
Compared to the rate-dependent behaviour of the bulk material, less efforts have been devoted to the study of strain rate-dependent resistance against fracture.Miehe et al. [30] suggested a phenomenological ansatz for the rate-dependent toughness in order to investigate the brittle-to-ductile fracture mode transition observed in the Kalthoff-Winkler experiment, i.e. for shear-loaded metals.Yin et al. [48] assumed the toughness of a linear elastic material to depend on rate of deformation.In their formulation, dissipation due to crack formation is incorporated into the free energy so that additional stress contributions are obtained from the rate-dependent fracture toughness.In these two models [30,48], rate-independent models for the deformation of the bulk are considered.To the best of the authors' knowledge, so far, there are no phase-field models that consider both a rate-dependent toughness and a model of rate-dependent deformation.
In the present contribution, a thermodynamically consistent phase-field model for fracture of materials with rate-dependent behaviour is presented.For this purpose, the previously introduced pseudo-energy functional [47], which consists of the free energy that includes a contribution related to viscous dissipative mechanisms, and the fracture contribution is advanced towards the finite viscoelasticity setting of Reese and Govindjee [43].Depending on the specific choice of the model parameters, the modelling approaches of [37], [39,40] or [42,44] can be retained as limiting cases.Based on experimental data for an Ethylene Propylene Diene Monomer (EPDM) rubber from the literature [39], the model parameters for the response of the bulk and the fracture behaviour are identified and model predictions are qualitatively and quantitatively verified on experimental results.In doing so, two assumptions for the fracture driving force, i.e. whether there shall be a contribution related to viscous dissipation or not, are investigated.With the aim of studying the possible influence of such a driving force component on the crack path, an asymmetrical setup is studied in addition to the symmetrical ones considered in recent publications, e.g.[39,42].Furthermore, based on experimental evidence on strain rate-dependent fracture toughness, cf.[3,49,50], and motivated from a phenomenological point of view, a rate-dependent fracture toughness is introduced.A numerical study on the coupling between rate-dependent resistance against crack propagation and viscoelastic bulk response is then performed.An overview on the structure of the proposed unified model is given in Fig. 1.
The paper proceeds as follows.In Sect.2, the proposed phase-field model of fracture in rate-dependent materials at finite deformation is presented and its thermodynamic consistency is proven.Subsequently, in Sect.3, algorithmic aspects are addressed.In Sect.4, the model parameterization is described and various numerical examples serve for validation and analysis of the model.A short summary and an outlook regarding the future work is given in Sect. 5.In the Appendix, information on the tangent for the local Newton iteration and the global material tangent is given.
Within this paper, italic symbols are used for scalar quantities (,  ) and bold italic symbols for vectors ().For Secondorder tensors, bold non-italic letters (T, τ) are used, whereas fourth-order tensors are written in Blackboard bold (C).

Phase-field formulation
In this Section, the phase-field model of fracture in materials with rate-dependent behaviour is presented.At first, the general energetic formulation of fracture in viscoelastic materials derived in [47] is extended to the finite deformation setting.Subsequently, the specific constitutive assumptions are outlined.

Finite viscoelasticity
Fracture phase-field Viscous driving force contribution Finally, governing equations are provided and thermodynamic consistency is proven.

Pseudo-energy functional
The variational approach to fracture.Following the pioneering work of Griffith [8], the dissipation due to crack growth  fr can be understood as an energetic quantity, which increases proportional to the crack surface.Accordingly, a pseudo-energy functional can be defined [7], wherein  0 ⊂ R  is the reference or undeformed configuration of the -dimensional domain under consideration and  0 ⊂  0 denotes the corresponding crack surface.The stored free energy is given by  sd and its density with respect to the reference configuration is denoted by .
The proportionality coefficient G c > 0 typically is referred to as fracture toughness.While G c is assumed to be a constant in the classical theory, in the recent literature, it is assumed to change during fatigue life, see e.g.[21], or due to plastic deformation [33,51].Furthermore, it can explicitly depend on the position in space in heterogeneous materials [52,53].In the following, the variational phase-field framework is set up for the case that G c is constant, first.Subsequently, the model is extended to account for a fracture toughness that depends on rate of deformation in Sect.2.4.For a given external load, the deformation of the domain and the crack surface  0 then can be determined from the equilibrium condition  → stat . ( In order to make this energetic approach to fracture accessible to numerical implementation, a regularisation of the functional  is introduced [9].For this purpose, cracks are described in a diffuse manner by means of a phase-field variable For sake of brevity, terms arising from external loads are omitted in ( 1) and what follows.
Diffuse representation of a crack within a domain that undergoes finite deformation which continuously varies from the intact ( = 0) to the fully broken ( = 1) material state.Using this variable, a crack surface density can be defined, cf.[12], in which the regularisation parameter ℓ c defines the characteristic width of the diffuse crack and the Nabla operator with respect to the reference coordinate is defined to with   denoting -th basis vector of the Cartesian reference coordinate frame.Fig. 2 illustrates the concept of diffuse crack representation.With this approximation of crack surface at hand, dissipation due to crack evolution can be expressed as wherein  fr is defined as density of fracture pseudo-energy with respect to a volume element in the reference configuration.
For the crack surface density  ℓc , several choices are possible, see e.g.[54].The expression adopted here typically is referred to as AT-2 model-with reference to the fundamental work of Ambrosio and Tortorelli [55].

Preprint submitted to Springer arXiv
In the regularised setting, the decrease of free energy due to fracture is expressed by means of the degradation function which has to fulfil the conditions Based upon the degraded reference free energy density the regularised functional of free energy is given by and the regularised counterpart of the pseudo-energy functional  reads Generalisation for inelastic material response.Following the previous work [47] and similar to phase-field fracture models for elasto-plastic materials, the free energy density is assumed to be additively decomposed into two essential ingredients.Naturally, the first one is the effectively stored strain energy st .In addition, in order to adequately account for the coupling between inelastic deformation and fracture mechanisms, a free energy contribution  vi related to accumulated viscous dissipation is assumed.In the latter, in order to keep the formulation as general as possible, the parameter  vi ∈ [0, 1] is introduced as a weight, cf.[36] and [37].The specific definitions of  st and  vi considered in this work are given in Sect.2.2.For a rigorous motivation and interpretation of  vi from a physical point of view, and a numerical investigation in the kinematically linear regime, the reader is referred to [47].Contributions similar to  vi are also considered in other recent phase-field models of fracture in viscoelastic materials [37,39].Furthermore, analogue terms are widely spread in modelling of failure in elasto-plastic materials [36,56,35], where a free energy contribution related to inelastic deformation, which Note that different from [36,56] and in line with e.g.[12,35], dissipation due to evolution of crack surface is not assumed to contribute to the free energy  , yet included as a distinct contribution  or  fr ℓc , respectively, to the regularised pseudo-energy functional  ℓc .
st and  vi can be understood as virtually undamaged densities of free energy with respect to a volume element in the reference configuration, i.e. the respective free energy which would be stored in such a reference volume element in the absence of damage.
is degraded in case of crack growth can be essential for the description of ductile fracture, cf.[29,28].
For the two contributions to the free energy, any degradation functions  st and  vi satisfying the conditions (8) can be considered, which, in general, do not have to coincide.In the literature, different approaches have been taken, e.g.quartic and cubic expressions [57,58], a sinusoidal ansatz [33,42], and parametric functions that include additional parameters, which can be fitted to the behaviour of a specific material [58,59,60,15].Without loss of generality, within the scope of this publication,  st () ≡  vi () ≡ () is assumed.Furthermore, the frequently adopted [9,12,34,37] quadratic function in which a small residual  is included in order to enhance numerical stability, is considered.

Viscoelastic bulk response 2.2.1 Kinematics
The displacement of a material point with coordinate  ∈  0 in the reference configuration is denoted by wherein is the motion function.Due to the diffuse approximation of crack topology, ( , ) can be assumed to be bĳective and continuous in space and time.The deformation gradient F and its determinant  are then given by For the rate-dependent deformation behaviour of the material, the approach of Reese and Govindjee [43] is pursued and a generalised Maxwell model is adopted as shown in Fig. 3.
In the non-equilibrium, or over-stress branch, deformation is assumed to consist of an elastic and an inelastic viscous portion, and the deformation gradient is multiplicatively decomposed into accordingly.Furthermore, following Flory [61], a decomposition of the deformation gradient into volumetric and isochoric parts is applied.For the equilibrium branch, the split is given by It has to be noted that there are alternative concepts for the phase-field modelling of ductile failure, also.For example, a degradation function  which, in addition to the fracture phase-field, depends on a measure of plastic deformation [29,32], and a fracture toughness that diminishes with accumulated inelastic strain [33], have been proposed.
Herein, without loss of generality, only one non-equilibrium branch is considered, which is sufficient for the material investigated in Sect. 4. The extension to multiple non-equilibrium branches can be done in a straightforward manner, though.

Preprint submitted to Springer arXiv
wherein I designates the second-order unit tensor and F is the isochoric portion of the deformation gradient.For the the nonequilibrium branch, F el and F vi are decomposed separately.Considering, for example, the elastic portion of deformation, its isochoric part is given by in which   ∈ {1, 2, 3} and  el  ∈ {1, 2, 3} are the number of pair-wise different principal stretches   and elastic principal stretches  el  , respectively.The second-order projection tensors, or eigenvalue-base tensors, are obtained from with the Kronecker delta    given by T and equivalent relations for P  , p el  and Pel  [62,63].

Specification of free energy densities
Strain energy.The strain energy density is additively decomposed into an equilibrium and over-stress part.Accordingly, for the virtually undamaged quantity  st , is defined, wherein C and F vi form the set of independent thermodynamic state variables considered here, in addition to the phase-field variable .Each contribution splits further into a volumetric portion vol  st,eq and vol  st,ov , and an isochoric part iso  st,eq and iso  st,ov , respectively.A compressible Ogden model [64] is assumed for both the equilibrium and nonequilibrium branches and the respective strain energy density contributions are defined to  st,eq = vol  st,eq + iso  st,eq wherein λ  =  −1/3   and λel  =  el −1/3  el  are the isochoric total and elastic principal stretches following from ( 18) and (19).Their algebraic multiplicity is given by   ∈ {1, 2, 3} and  el  ∈ {1, 2, 3}, respectively.Furthermore, the compression moduli are denoted by  eq > 0 and  ov > 0, and  eq O ,  eq  ,  eq  > 0, as well as  ov O ,  ov  ,  ov  > 0 are parameters of the Ogden models.From these constants, the initial shear moduli and the according Poisson's ratios can be defined.For example, for the equilibrium branch, they read and similar relations hold for the non-equilibrium branch.
Viscous contribution.The degraded free energy contribution  vi related to inelastic mechanisms is designed such that a certain portion of accumulated viscous dissipation can enter the phase-field fracture driving force.Before defining If there are three pair-wise different principal stretches, i.e.   = 3 or  el  = 3, the projection tensors can also be represented by means of the eigenvectors in a straightforward manner, e.g.p  =   ⊗   with   denoting the -th eigenvector of b.
Preprint submitted to Springer arXiv the respective virtually undamaged quantity  vi in the finite viscoelasticity framework, the simple setting of a uniaxial deformation in the kinematically linear regime is considered for motivational purpose.Then, in the absence of damage, viscous dissipation in a material described by means of the generalised Maxwell model takes the form in which  vi is the rate of inelastic deformation and  designates the viscosity of the material.In order to generalise  vi,1D , the tensor is introduced as a measure of the rate of inelastic deformation in the finite viscoelasticity setting, wherein is the Lie derivative of b el .Furthermore, a fully symmetric, positive definite, isotropic fourth-order tensor is defined, in which I D , is the fully symmetric fourth-order deviator projection tensor.Therein, iso , vol  > 0 are viscosities with respect to the isochoric and volumetric portion of deformation, respectively.The virtually undamaged free energy density contribution related to viscous mechanisms is then defined to which is positive and monotonically increasing in time.
Remark on the measure of rate of inelastic deformation.The definition of d vi (28) can be written in an alternative form, which may be more intuitive.For this purpose, the inelastic is introduced.It refers to the intermediate configuration defined by F vi .The counterpart of lvi transformed to the current configuration reads Assuming that there is no inelastic spin, i.e. lvi = dvi with dvi = sym lvi denoting the rate of inelastic deformation with respect to the viscous intermediate configuration, (29) 1 can be rewritten as and holds, from which, together with the transformation rule (34), the definition of d vi as an Eulerian measure of rate of inelastic deformation becomes clear.For more details, the reader is referred to [65], where similar kinematic relations are derived in the context of plasticity.

Evolution of phase-field
The equation governing the evolution of the fracture phasefield variable is deduced from the pseudo-energy functional  ℓ c by means of the variational derivative wherein  f is introduced as a kinetic fracture parameter in order to avoid discontinuity of the field variables in time and for numerical purposes, i.e. for enhancing the stability of the solution scheme, cf.[66,34] and  denotes the outwardpointing unit normal vector on  0 .For the simulations presented in Sect.4,  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.
Inserting the definitions made in the previous Sections into  ℓ c , the evolution equation (37) 1 takes the form from which it becomes clear that, depending on the specific choice of  vi , fracture is driven by stored strain energy and the free energy contribution related to a portion of accumulated viscous dissipation.It has to be noted that in the present form (38), the evolution equation enables the phase-field variable to decrease, i.e. crack healing is not prohibited.Therefore, a modification is adopted which overcomes this issue, see Sect.3.1.

Rate-dependent fracture toughness
For various materials, in addition to or instead of the deformation behaviour of the bulk material, the resistance against fracture has been reported to depend on rate of deformation.
In several other models, e.g.[67,39],  f is assigned a finite value and thus considered as an additional material parameter.On the one hand, such a direct coupling of rate effects into the evolution of phase-field by means of  can enable more modelling flexibility especially regarding the post-critical stage of a response.On the other hand, when it comes to damage, incorporation of a finite  f is equivalent to assuming a pseudo-viscous dissipation in addition to proper viscous effects and fracture dissipation.However, for fracture dissipation, according to the fundamental modelling hypothesis (1), the fracture toughness G c is assumed to be the essential parameter.Therefore, a toughness depending on rate of deformation is presumed to be more consistent from an energetic point of view if a direct coupling of rate-effects into phasefield evolution is necessary.Furthermore, a finite  f would also incorporate some redundant information which should rather be taken into account by the viscoelastic model for deformation.

Preprint submitted to Springer arXiv
For instance, in elastomers, at low rates of deformation, chain entanglements can be resolved, which is not the case at high rates of deformation.Therefore, the number of chemical bonds that are broken when a crack propagates can be assumed to rise with rate of deformation and the fracture toughness increases accordingly, cf.[68,2] for a more detailed discussion and experimental results.Furthermore, for several natural materials and foods, where the underlaying microscopic mechanisms can be more complex, a rate-dependency of G c has been reported [3,49,50].Therefore, as an extension of the variational phase-field equation (38), fracture toughness is considered to depend on deformation rate, which enables a maximum of flexibility in modelling rate-dependent fracture processes.For this purpose, is introduced as a scalar measure of effective rate of deformation.Furthermore, without loss of generality, in line with [69], the sigmoid-shaped function is adopted, see Fig. 4. The extended phase-field evolution equation can then be written as In Sect.4.3, for different parameters G 1 c , G 2 c , ,  ref , coupling between rate-dependent deformation and toughness is analysed.
Fracture pseudo-energy and rate-dependent toughness.If fracture toughness is a function of rate of deformation and thus implicitly depends on time, density of fracture pseudo-energy  fr has to be rewritten as (42) in order to account for possible a posteriori changes of G c following the evolution of fracture phase-field at a given time t.Therefore, the phase-field equation ( 41) can be seen as a nonvariational extension of (38), similar to the suggestions made in [31,70], for instance.Furthermore, it is noted that, different from e.g.[48], in the proposed model, rate-dependency of fracture toughness does not affect the density of free energy  , since dissipation due to crack growth is not supposed to enter  .Accordingly, no additional stress terms arise from rate-dependent toughness, see the evaluation of the entropy inequality below in Sect.2.5.
It has to be noted that, unlike here, in some publications dealing with fracture of inelastic materials, G c is not only regarded as a measure of dissipation directly coming along with breaking of bonds but also comprises dissipative mechanisms of the bulk material.
2.5 Stress tensor, viscous evolution and thermodynamic consistency Under isothermal conditions, the second law of thermodynamics can be stated by means of the density of dissipation power  as cf. [71], with T denoting the second Piola-Kirchhoff stress tensor.For  =  (C, F vi , ), this inequality can be rewritten to wherein the contributions to dissipation power density due to fracture,  fr , and viscous effects,  vi , can be identified.The standard argument that  ≥ 0 shall hold for arbitrary processes leads to the definition of stress with the virtually undamaged equilibrium and over-stress tensors denoted by 0 T eq and 0 T ov , respectively, and the residual inequalities Stress tensor.Inserting the definitions of  st,eq and  st,ov , Eqs. ( 25) and (26), into (45), the contributions to the second

Preprint submitted to Springer arXiv
Piola-Kirchhoff stress tensor take the form Residual inequalities.As both  st and  vi are positive,  vi ∈ [0, 1], and due to (8) 3 , the condition  fr ≥ 0 reduces to  > 0, i.e. irreversibility of fracture.The fulfilment of this demand will be addressed in Sect.3.1.
Due to (7),  vi ≥ 0 reduces to Making use of the relations outlined in Sect.2.2, after some lengthy manipulations, the first term can be rewritten as is the virtually undamaged Kirchhoff over-stress.For a more detailed derivation see also [43].Then (49) takes the form from which, in line with [43], the equation governing viscous is defined.By reason of  vi ∈ [0, 1], the quadratic form obtained from inserting (53) into ( 52) is compatible with the second law of thermodynamics.

Irreversibility of fracture
In order to guarantee irreversibility of fracture, the history variable approach of Miehe et al. [12] is pursued.For this purpose, the phase-field equation ( 41) is rewritten to Table 1: Governing equations for the present model following the total Lagrangian approach: Balance of linear momentum (a), phase-field equation (b), viscous evolution (c), ratedependent fracture toughness (d).Without loss of generality, volume forces are neglected in (a).
wherein the history variable comprises the maximum of virtually undamaged fracture driving force which has occurred.With this form of the phase-field evolution at hand, the governing equations of the model are summarised in Tab. 1 considering the total Lagrangian approach.Alternatively, in line with [13,72], Dirichlet boundary conditions can be applied to the phase-field on all nodes where the phase-field variable has reached a critical value  crit : For the setups analysed in Sect.4, the two approaches have been compared, exemplary, and no relevant differences could be noticed.

Viscous evolution
For the integration of viscous evolution equation ( 53), an operator split scheme of predictor-corrector type is adopted For the simulation of relaxation-dominated load cases together with  vi > 0, special attention has to be paid to the fact that  vi is incorporated into H according to (55).Therefore, in these cases, either an altered definition of the history variable or the Dirichlet boundary condition approach would be more reasonable, cf.[47].
as proposed in [43].Within the scope of this well-established approach, the evolution of elastic deformation is split into the contributions from change in total deformation, which is considered in the predictor step, and viscous evolution, which is accounted for in the inelastic corrector.For the predictor step, viscous deformation F vi or C vi is frozen, giving a trial state of elastic deformation at time step   to Subsequently, within the corrector step, ( 58) is evaluated for the total deformation assumed to be constant, i.e. l = 0, which, with evolution equation ( 53) and kinematic relations ( 28) and ( 29) can then be written as Due to isotropy, the principal directions of b el , b el tr and 0 τ ov coincide, which makes the evaluation of (60) in terms of elastic principal stretches  el  attractive.For the viscosity tensor V defined according to (30), this leads to  are determined from an iterative solution of the system of non-linear algebraic equations   = 0 with  ∈ [1, ] ⊂ N.However, in case of two-dimensional plane stress setups as considered in Sect.4, in addition to these three equations, it has to be ensured that the out of plane stresses vanishes, i.e. 0 τ ov 3 = 0 must hold.In these cases, in addition to  el  , the out of plane stretch  3 has to be determined from the system of equations Regardless of whether a plane stress state is considered or not, the respective system of equations ( 62) or ( 63) is solved by means of a local Newton iteration scheme at each quadrature point.In the following, the procedure is briefly described for the case that a plane stress state has to be guaranteed.With the vector of unknowns then written as Initialisation at each increment : the linearisation of (63) around    is given as For the specification of the derivatives   /  , the reader is referred to Appendix A. Based on the linearisation, the Newton procedure is carried out as summarised in Algorithm box 1.
For this, at each increment   , the iteration scheme is initialised by means of with Within the iterative solution procedure, special attention has to be paid to tr  el 3 as it needs to be updated after each local iteration  according to due to the change of  3 = ln  3 .

Weak forms of the governing equations
For the derivation of the weak forms of equilibrium and phase-field equation, the test function spaces and are defined.Therein, H 1 ( 0 ) is the Sobolev space of square integrable functions possessing square integrable derivatives in  0 , and  0   denotes the parts of the boundary where the -component of the displacement vector  is prescribed.Then, (a) and (b) from Tab. 1 are multiplied by and  ∈ W  , respectively.Integration by parts and making use of the divergence theorem yields wherein p denotes the Piola traction vector with its components p  prescribed on  0 \  0   , and d = 0 .
(74) Time discrete forms 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.The coupled problem is solved by means of a staggered approach.Furthermore, adaptive control of the time step size is employed based on a heuristic scheme.Information on the material tangent that is required for the iterative solution of ( 73) is given in Appendix B.

Representative simulations
In this Section, several numerical examples are presented in order to analyse the characteristics of the present model and to demonstrate its flexibility in describing different responses.Furthermore, the comparison of numerical predictions to experimental results of Loew et al. [39] serves for validation of its predictive capabilities.

Parameter identification from experimental data
Within this publication, the viscoelastic behaviour of EPDM rubber is considered that has been experimentally analysed in [39].
Bulk response.At first, the parameters describing the deformation of the bulk material are determined.For this purpose, the averaged stress-stretch curves from [39, Fig. 6] are considered as depicted in Fig. 5.For three different rates of deformation, these curves have been identified from uniaxial tension tests with dumbbell specimens.A homogeneous uniaxial stress state is assumed and damage is not taken into account, here.Furthermore, as no information on deformation in transversal direction is available,  eq =  ov = 0.48 is set in order to account for the high resistance against volumetric deformation that is typically observed for rubber.For isochoric and volumetric deformation, an identical relaxation time is assumed.The Ogden parameters  eq  ,  eq  ,  ov  ,  ov  as well as  are then identified by means of minimising the deviation between experimental data and model prediction.In doing so, following [64, p. 305], it is demanded that the constants satisfy the requirements for any  and similar constraints for the non-equilibrium branch.
To this end, in Matlab R2020b, the GlobalSearch strategy together with the fmincon algorithm for constrained optimisation problems is employed.For an adequate approximation of the material behaviour, two Ogden exponents have revealed necessary for both the equilibrium and over-stress branch, respectively, i.e.  eq O =  ov O = 2.The parameters obtained are summarised in Tab. 2. From Fig. 5 it becomes clear that the finite viscoelasticity formulation together with the Ogden approach allows for a very good approximation of the experimental results over the entire range of stretch  ∈ [1, 2.5] that has been experimentally investigated.Furthermore, the present model enables to capture the rate-dependent response in a more reliable manner then the linear viscoelasticity model based on the Yeoh-type strain energy density [39].
Identification of the fracture parameters.With the calibrated bulk deformation model at hand, the fracture phase-field is parameterized from SENT experiments, i.e. specimens with A proof of uniqueness of the parameters identified, i.e. a global minimum of discrepancy between model and experiment, can not be provided.Nevertheless, GlobalSearch involves minimisation for a huge number of different start values in order to obviate local minima.
An increase of the number of Ogden branches to  eq O =  ov O = 3 did not lead to a perceptibly better approximation.
It has to be noted that the rate-dependency perceptible in Fig. 5 is not too pronounced.Accordingly, it could be worth investigating a broader range of stretch rates, since the rate-dependent behaviour of EPDM rubber can play a crucial role when it comes to failure, e.g. in case of creep fracture.Furthermore, additional experiments such as relaxation or creep tests could allow for differentiating between equilibrium and non-equilibrium contributions to stress in a significantly more accurate manner.However, within this contribution, we proceed with the experimental results available in the literature.
Preprint submitted to Springer arXiv   [39] vs. present model for three different stretch rates  a single pre-existing notch under tension.These experiments have been conducted at two rates of prescribed displacement [39].The according specimen geometry is depicted in Fig. 6.
For the numerically motivated kinetic fracture parameter and the residual stiffness, the values  f = 10 −4 Ns/mm 2 and  = 10 −10 , respectively, are chosen.In a convergence study, these values have revealed sufficiently small so that the influence of  f and  on the simulation results vanishes.The regularisation parameter is set to ℓ c = 0.275 mm, which is identical to [39].In order to enable a step-by-step analysis of the model, a constant fracture toughness is assumed, here, and G c (d) according to (40) is investigated in Sect.4.3.Furthermore, with the aim of performing a thorough analysis of viscous fracture driving force contribution in Sect.4.2, the two limiting cases  vi = 0 (approach A) and  vi = 1 (approach B) are considered.Under these two assumptions, the respective values of G c are identified from experimental data.For this purpose, regarding the critical deformation in SENT for the two rates experimentally investigated, deviation between simulation and mean values from the experiments is minimised by means of a gradient-free approach.Since the specimens are of low thickness, plane stress conditions are assumed and two-dimensional simulations are performed, here.Due to symmetry, only one half of the SENT specimen is considered.The mesh consists of quadratic triangular elements and is refined along the crack path.ℎconvergence is verified.The optimal simulation results are compared to the range of experimental data in Fig. 7 and the parameters of the fracture model are summarised in Tab. 3.For both  vi = 0 with optimal G c = 10.7 N/mm, and  vi = 1 with optimal G c = 12.0 N/mm, good agreement between simulation and experiment can be stated.With  vi = 1, a marginally better approximation is obtained for this setup.However, in both cases, the critical force is slightly overestimated.Furthermore, especially for the higher rate ū = 3.328 mm/s, the simulated - curves do not completely reproduce the smooth decrease experimentally observed in the post-critical stage preceding complete failure.Instead, the critical point is followed by a sudden drop of reaction force that, interestingly, does not come along with complete failure yet.It corresponds to crack initiation at the tip of the pre-existing notch, see Fig. 8, and is succeeded by a smoother decrease of force for which crack propagation through the specimen involves a slight increase of external load before, finally, it comes to complete failure.To the best of the author's knowledge, such a phenomenon does not arise in hyperelastic models, whereas it also has been reported for linear viscoelasticity [39,47].The effect is the more pronounced the lower ū.Obviously, it is provoked by the rate-dependent behaviour of the bulk material that involves an increase of effective stiffness as well as the effective load bearing capacity of the material when, locally in the vicinity of the crack, rate of deformation suddenly raises up due to the initiation of fracture.For a rigorous analysis within the small strain context, the reader is referred to the previous work [47].

Model validation and analysis of viscous driving force
For further model validation and analysis, double notched specimens under tension (DENT) with varying length of the pre-existing notch  are considered as depicted in Fig. 9.
A straightforward way for tuning the model such that it would better reproduce this specific experimental observation could be defining a finite  f 0, which leads to a smooth decrease of post-critical  - curve instead of a sudden jump, see e.g.[67,Fig. 9].However, as outlined in Sect.2.3, this approach has some important drawbacks which is why it is not pursued here.For a more expressive investigation, it may be eligible to thoroughly elaborate on crack initiation mechanisms.For example, cavitation or void formation are often observed in rubbery polymers, see e.g.[73], and modified fracture phase-field models that explicitly take these mechanisms into account have recently been proposed in [74,75], wherein hyperelasticity is assumed for the bulk.For comparison of model prediction with experimental data from [39], symmetrical specimens are considered, i.e.  = 75/2 mm.[39] for  ∈ {9, 5} mm and a constant rate ū = 75 mm/min in Fig. 10.For both approaches A and B, model predictions fit the experimental results well, which is also true for  ∈ {7, 3} mm (not depicted).The good agreement demonstrates the predictive capability of the present model and the suitability of the parameter identification from experiments with homogeneous and single-notched specimens.
With the aim of more thoroughly analysing the ratedependency of responses and elaborating on the driving force contributions, additional simulations are performed for  = 7 mm and various rates ū ∈ [12.5, 400] mm/min.The numerical predictions for the two approaches  vi = 0 and  vi = 1 are compared in Fig. 11.Regardless of the approach for the driving force, for high displacement rates, the responses converge against an upper elastic limit for which there is almost no viscous dissipation until failure.For very low ū, the responses of the structure likewise approach a lower elastic limiting case where over-stresses do approximately vanish during entire simulation.In between, for intermediate displacement rates, the critical displacement level diminishes with ū for both approaches A and B. In contrast, regarding the rate-dependency of critical force level, the model predictions do significantly differ depending on whether a viscous fracture driving force contribution is assumed or not.On the one hand, critical force monotonically increases with rate when there is no such contribution, i.e.  vi = 0 (A).On the other hand, for  vi = 1 (B), critical force becomes minimal for intermediate ū, for which the greatest critical values of  vi are observed, see Fig. 12.
Although no experimentally-determined force-displacement curves are available, it can be stated that the former is in agreement with experimental observations [39], whereas the latter contradicts experimental experience.At least when modelling fracture of elastomeric materials under monotonic loading, in some cases, fracture driving force contribution associated to accumulated viscous dissipation can thus lead to erroneous model predictions.In other words, modelling approach A has revealed more plausible, which, in a sense, is different from plasticity, where a fracture driving force related to inelastic mechanisms has revealed advantageous [29,36].Interestingly, such an observation has not been made in the previous study within the small strain framework [47], where a less pronounced influence of viscous effects on crack propagation has been observed.This can probably be attributed to the fact that the present formulation enables to describe larger deviations away from thermodynamic equilibrium, resulting in considerably greater viscous contributions to fracture driving force.
As it has been comprehensively described in [47], it essentially is the change of effective stiffness and the amount of dissipation until failure that lead to the change of critical force and displacement level with rate of external load.While the Figure 12: Symmetrical DENT-free energy contribution related to viscous dissipation (approach B) for various rates and a fixed size of pre-existing notch  = 7 mm amount of fracture driving force necessary for crack growth remains constant, the fracture driving force available for a constant level of deformation can change with rate.On the one hand, effective stiffness of the viscoelastic material monotonically increases with increasing rate of deformation.For a certain external displacement ū prescribed, the density of strain energy raises with ū, accordingly.On the other hand, in case of monotonic loads, the amount of viscous dissipation and thus, in case of  vi > 0, the level of  vi at failure becomes maximal for intermediate rates.
Although viscous fracture driving force contribution has revealed not suitable for describing failure of elastomers under monotonic loads, it might be suitable for other classes of materials, e.g.thermoplastics, and especially for the modelling of fatigue fracture, e.g. with 0 <  vi 1.In composites and thermoplastic materials, for instance, viscous dissipation and self-heating mechanisms can have an important influence on fatigue life, cf.[76].
Crack patterns in asymmetrical specimens.In addition to the symmetrical specimens, simulation results are presented in the following for an asymmetrical DENT geometry as depicted in Fig. 9 with  = 27.5 mm and  = 9 mm.Since for ductile fracture of metals, where instead of viscoelasticity another For example, in the phase-field fatigue fracture model [77], which is applied to a rubbery polymer, a fatigue fracture driving force is introduced that also incorporates viscous dissipation.However, similar to [39], a model of linear viscoelasticity at finite deformation is used which does not allow for separation of accumulated viscous dissipation and non-equilibrium part of stored strain energy.As a consequence, entire viscous dissipation is also included in the quasi-static fracture driving force contribution.class of dissipative materials is involved, the choice of fracture driving force revealed crucial the appropriate numerical description of asymmetrical crack patterns, cf.[29], simulations are performed for both approaches A and B. The corresponding force-displacement curves are depicted in Fig. 13.The overall rate-dependency of the specimen response is identical to what has been described above for the symmetrical geometry.In particular, for  vi = 1, the numerically predicted critical force becomes minimal for an intermediate rate of external displacement, which does hardly coincide with what would be observed in experiments.In Fig. 14, the final crack patterns are compared for ū = 200 mm/min.In order to ease comparison, the phase-field is shown with respect to the reference domain  0 .For both fracture driving forces A and B, the crack pattern predicted for the viscoelastic material is essential different from what is typically observed when metals fail in a ductile manner.Instead of a single crack that connects the two pre-existing notches, two cracks independently propagate through the specimen.At a certain length, one of the two stops to propagate, resulting in an asymmetrical final crack pattern, see Fig. 15.Regardless of  vi and ū, qualitatively identical crack paths are predicted.However, depending on  vi , slight differences concerning the final length of the shorter crack can be stated especially for intermediate rates.Interestingly, when critical force is reached, the two cracks suddenly propagate over a finite width, which comes along with a significant abrupt drop of force.For intermediate and higher rates, similar to SENT geometry, a slight increase of external displacement ū is necessary to make one of the cracks propagate further, resulting in a less heavy slope of the force-displacement curve before it finally comes to catastrophic failure.For these higher rates, in the simulations there is a stage that can be seen as a kind of stick-slip-like crack propagation, where the crack tip suddenly advances over a finite distance and then arrests over and over again.These effects also lead to a non-smooth - curve in the post-critical range.Interestingly, for very small ū, such a behaviour is not simulated.In the literature on dynamic crack growth, comparable phenomena have been reported, cf.e.g.[79].However, it has to be noted that regarding this particular aspect, the predictive capabilities of the present model are somewhat limited, as inertia effects are not taken into account.
For the EPDM rubber for which the model has been parameterized here, no experimental results are available for crack propagation in asymmetrical specimens.Nevertheless, the crack patterns simulated with the present model are in excellent agreement with what has been observed in experiments for other viscoelastic materials, see e.g.[78].It is obvious that, when specimen geometries are similar, these crack patterns in viscoelastic materials can differ from the ones that form in elasto-plastic ones, since the inelastic mechanisms are essentially different.For example, there typically is no zone of inelastic localisation in viscoelastic materials whereas localisation of plastic deformation can play an important role For all the simulations performed, it is always the right crack tip which stops propagating at a certain length.It is deemed likely that this is due to the non-symmetric mesh that has been used for all the computations.when it comes to ductile fracture of metals.
Creep fracture.In addition to fracture under monotonically increasing loads, a qualitative analysis of creep fracture is performed by means of one representative example.For this purpose, the symmetrical DENT geometry with  = 7 mm is revisited.Instead of ū, a traction force  is applied that linearly increases with time until a certain value  max is reached and is hold constant, subsequently.For two different values of  max , boundary conditions and model predictions are depicted in Fig. 16 for both approaches A and B. It can be stated that, generally, creep fracture can be captured regardless of the value of the assumption made on fracture driving force.In case of  vi > 0, failure can occur for lower  max and after a shorter amount of creep time than for  vi = 0. Furthermore, if a fracture driving force contribution from viscous dissipation is assumed, it can also depend on the rate  if creep fracture is predicted, since viscous dissipation vanishes for very small , see [47] for a discussion in the small strain context.

Investigation of rate-dependent fracture toughness
In the foregoing Section and the previous work [47], it is demonstrated that within the scope of an energetic phase-field fracture approach, a rate-dependent material model for the bulk induces a certain relationship between critical load and rate of deformation when G c is constant.Therefore, in addition to experimental indication [68,2,3,49,50], there also is a clear motivation for assuming a rate-dependent toughness from a phenomenological point of view.Assuming The case of G c increasing with rate of deformation is investigated first.As a representative example, the specimen response is depicted in Fig. 17  .Accordingly, in the post-critical range of the - curves, a slightly less sharp slope can be observed with respect to G c = const.However, this effect is not pronounced compared to the effects arising from the rate-dependent toughness when pre-critical rate of deformation  becomes close to the threshold value  ref .
In that case, deformation at failure begins to raise with rate similar to stiffness and critical force.Experimentally, similar effects can be observed for some natural materials, see e.g.[80] for an overview, as well as viscoelastic silicone elastomer based model systems [81].
For the discussion of G c decreasing with , G 2 c = G 1 c /4,  ref = 200 s −1 ,  = 10/ ref , are considered, exemplary.From the force-displacement curve depicted in Fig. 18 it appears that for ū 300 mm/min, the responses do again coincide with the case G c = G 1 c = const.Naturally, the initiation of the phase-field crack at the notch tips is immediately followed by complete failure, since in this moment, the sudden increase in rate of deformation comes along with a drop of toughness.Nevertheless, for the DENT geometry, similar behaviour is obtained as simulation result for G c = const., which is in agreement with experiments.For high rates ū ≥ 400 mm/min, where   ref also holds in pre-critical range, the decrease of deformation of failure that stems from the rate-dependent stiffness of the viscoelastic material is further intensified by the rate-dependent fracture toughness.In addition, critical force does no longer raise up with ū yet also decreases.For sugar-based confections [82], a similar characteristic behaviour has been observed very recently.For high displacement rates, these materials fail in a brittle manner, i.e. at small deformation as well as low external force, whereas at low rates, they can It has to be noted that, when crack propagation takes place, quantitative predictive capability regarding the rate of deformation is somewhat limited for present formulation, since inertia effects are not taken into account.

Conclusion and outlook
For the simulation of fracture of materials with ratedependent behaviour, a flexible phase-field model is presented.To this end, the theory of finite viscoelasticity [43] is adopted for the deformation of the bulk material.The phase-field model is formulated such that, depending on the choice for the parameters, a portion of viscous dissipation can enter the fracture driving force.Moreover, in addition to the viscoelastic A publication on experimental and numerical investigation of this brittleto-ductile fracture mode transition is in preparation.In order to analyse the coupling between different rate effects, a gradual analysis of the model is performed.The model of finite viscoelasticity is parameterized for an EPDM rubber based upon stress-deformation curves from the literature.Ogden-type strain energy densities are considered for both the equilibrium and over-stress parts of the response and very good agreement of the model with experimental data is obtained.Assuming a constant fracture toughness for the EPDM rubber, two limiting cases are studied regarding the fracture driving force and the respective values of toughness are identified from experimentally-determined SENT force-displacement curves.In doing so, either entire viscous dissipation or only effectively stored strain energy is assumed to enter the fracture driving force, respectively.In the absence of a driving force contribution related to viscous dissipative mechanisms, very good agreement between model predictions and experiments can be stated for different setups.In this case, plausible results are obtained over a broad range of rates of external load and deformation, respectively.On the contrary, if viscous dissipation is assumed to enter fracture driving force, erroneous model predictions can arise, here.In this case, agreement with experimental data is obtained for some specific rates, only.Accordingly, different from e.g.phase-field modelling of ductile fracture in metals, a distinct fracture driving force contribution related to inelastic dissipative mechanisms as proposed in [39,40] or [37] has revealed not favourable for viscoelastic materials, in particular not for rubbery polymers.Furthermore, comparing the crack paths predicted in asymmetrical DENT specimens, it is demonstrated that such a driving force con-tribution is not necessary in order to predict non-symmetric crack patterns in an appropriate manner.

Preprint submitted to Springer arXiv
By means of a numerical study, it is demonstrated that a rate-dependent fracture toughness can significantly increase the capability of the phase-field model in capturing varied experimentally-observable responses.In particular, it seems suitable to describe rate-dependent brittle-to-ductile fracture mode transitions.In contrast, in case of a constant toughness, the rate-dependent model of bulk deformation induces a certain rate-dependency of critical stress and deformation, which does not coincide with experimental evidence for some specific materials.At least from a phenomenological point of view, rate-dependent fracture toughness thus seems to be an essential tool for modelling of rate-dependent fracture phenomena.While this contribution clearly demonstrates the potential of a rate-dependent fracture toughness within the proposed model, a quantitative description of rate-dependent brittleto-ductile fracture mode transitions is beyond its scope.A thorough experimental analysis of these effects in materials with rate-dependent deformation behaviour, e.g.caramelbased confections [82], as well as a quantitative description based upon the framework presented in this contribution are the subject of current work.have to be evaluated.

B Material tangent
The consistent Lagrangian material tangent =: () 0 C eq + 0 C ov is determined in order to enable the iterative solution of the weak form of balance of linear momentum (73).In line with e.g.[83], the derivation of the tangent is performed assuming   =  el  = 3 and the case of identical principal stretches or elastic principal stretches is then a posteriori addressed by means of L'Hôpital's rule.The equilibrium part of the virtually undamaged tangent is given by 0 C eq =2  0 T eq C = ∑︁ ∈ {1,2,3}  ∈ {1,2,3} wherein   denote the orthonormal eigenvectors of C, see e.g.[62,63] or [84] for a derivation of the derivatives of principle stretches and projection tensors.Into this expression (84), for the specific model under consideration, It is assumed that an appropriate orthonormalization method is used in case of multiple principal stretches.can be inserted.In case of multiple principal stretches, i.e. ∃ ≠  :   =   , the second term in (84) can be evaluated making use of L'Hôpital's rule [83] lim and the two contributions in (97) specified by ( 78) and (79).
In case of multiple elastic principal stretches, i.e. ∃ ≠  :  el  =  el  , for the treatment of the second term in (92), the same procedure applies as outlined above for the case of multiple principal stretches.In particular, L'Hôpital's rule reads (98) with the respective derivatives given by (93).

Figure 1 :
Figure 1: Modular structure and flexibility of the proposed model for rate-dependent fracture phenomena

wherein
el = det F el .For the specific definition of the material model, the positive definite left and right Cauchy-Green deformation tensors, b = F • F and C = F • F, as well as their elastic counterparts b el = F el wherein 0 τ ov,dev  denote the principal components of the over stress deviator dev 0 τ ov .Within the scope of the FE framework, differential equation (61) is integrated in an approximate manner by means of an exponential mapping ansatz and rewritten in terms of logarithmic elastic principal stretches  el  = ln  el  as 0 =  el  + Δ 1 2 iso  0 τ ov,dev  + 1 9 vol  tr 0 τ ov −  el tr  =:   .(62) Generally,  el

𝑛𝑗 + 1 endAlgorithm 1 :
=   =0    , =0     =    =0    while     ∞ > tol do Solution of the linearised system of equations: Local Newton iteration scheme in case of plane stress state and the local tangent matrix

Figure 5 :
Figure 5: Stress response of the bulk material under homogeneous uniaxial tension-Experimental data [39] vs. present model for three different stretch rates

Figure 6 :
Figure 6: SENT-Setup considered for the identification of G c mm min −1 ) : range of experimental data approach A approach B average critical point in experiments increments depicted in Fig.8

Figure 13 :Figure 14 :
Figure 13: Asymmetrical DENT-comparison of specimen responses for approaches A and B for various rates, pre-notch position  = 27.5 mm and notch length  = 7 mm

Figure 15 :
Figure 15: Asymmetrical DENT-phase-field crack initiation and propagation through the specimen for ū = 200 mm/min and  vi = 0 (approachA).The corresponding force-displacement curve is depicted in Fig.13.Qualitatively similar results are obtained for approach B and other ū.
G c to be a function of effective rate of deformation  = d F enables more flexibility in describing the rate-dependent failure of varied materials.In what follows, this is demonstrated by means of numerical studies considering both an increase and a decrease of G c with .For this purpose, the DENT setup with  = 7 mm and  vi = 0 is revisited.For G c , the sigmoid-shaped function (40) is assumed with G 1 c = 10.7 N/mm and  vi = 0 as parameterized for EPDM whereas the responses for different G 2 c > G 1 c as well as G 2 c < G 1 c are investigated.Apart from that, the parameters are identical to the ones listed previously.

Figure 16 :
Figure 16: DENT-boundary conditions and model predictions for the investigation of creep fracture

Figure 17 :
Figure 17: DENT-model prediction in case of fracture toughness G c assumed to increase with effective rate of deformation  undergo large deformation.

Figure 18 :
Figure 18: DENT-model prediction in case of fracture toughness G c assumed to decrease with effective rate of deformation

1 𝑛 2 N
], for the derivation of 0 C ov , a virtually undamaged over-stress tensor0 Tov = −1 F vi • 0 T ov • −1 F vi (87) is introduced with reference to the intermediate configuration defined by the viscous deformation gradient of the previous time step −1 F vi , i.e. based on the decomposition of the deformation gradient at increment  into  F =  F el tr • −1 F vi .stress part of the virtually undamaged material tangent then can be written as 0 C ov     = 2 −1  vi   −summation convention applies for double indices.From this, the over-stress tangent in terms of the intermediate configuration described by −1 F vi can be defined to  ⊗ N ⊗ N  ⊗ N + N ⊗ N  , (92) Preprint submitted to Springer arXiv wherein N  denote the orthonormal eigenvectors of Cel tr and 0 Tov  are the eigenvalues of 0 Tov .The first term in (92) the local Newton iteration has converged towards zero, is made.This assumption leads to by (77) and further specified in Appendix A. Accordingly, the derivative  0 τ ov  / el tr  in (93) can be expressed as  0 τ ov • F el and Cel = Isotropy of the material is assumed and the constitutive equations are specified in terms of principal stretches   and  el  , which are obtained from the spectral decompositions b

Table 2 :
Parameters of the finite viscoelasticity model for the deformation of the bulk material

Table 3 :
Parameters of the phase-field model calibrated for EPDM rubber with G c = const.assumed