{'fraction_non_alphanumeric': 0.03773791692543747, 'fraction_numerical': 0.014868666203923575, 'mean_word_length': 4.646452145214521, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fatigue phase-field damage modeling of rubber using viscous dissipation: Crack nucleation and propagation

By regularizing sharp cracks within a pure continuum setting, phase-damage models offer the ability to capture crack nucleation as well as crack propagation. Crack branching and coalescence can furthermore be described without any additional efforts, as geometrical descriptions of the cracks are not required. In this contribution, we extend our previous phase-field model for rate-dependent fracture of rubbers in a finite strain setting (Loew et al. 2019) to describe damage under cyclic loading. The model is derived from the balance of mechanical energy and introduces a fatigue damage source as a function of the accumulated viscous dissipation under cyclic loading. We use uniaxial cyclic tension to present the influence of the fatigue material parameters and to confirm the model’s energy balance. The parameters are subsequently identified using monotonic and cyclic experiments of a plane stress nature. Finally, the model is validated by separate experiments, which demonstrate that the model accurately predicts (fatigue) crack nucleation as well as propagation.


Introduction
Virtual prototypes can be used to reduce the financial burden of physical prototypes in design processes. They also allow more designs to be evaluated but require accurate mechanical models to describe the three stages of (fatigue) failure: crack initiation, crack propagation and final fracture.
The lifetime of rubber parts is affected by many factors ( [19], [20]), which can be subdivided into four classes: mechanical load history, environmental conditions (such as temperature), rubber formulation and constitutive behavior [20]. This contribution focuses on the mechanical load history and constitutive behavior.
Hence, we test one rubber compound and ignore temperature differences and chemical aging.
Several fatigue damage models for rubber have been published: [31] presented a continuum damage model and defined the fatigue life as a function of the strain amplitude. Multi-axial loading was considered by [4], who used the cracking energy density criterion of [18]. [11] and [25] extended the theory to thermo-viscoelastic solids to describe self-heating under cycling loading. All these continuum damage models show a mesh sensitivity [27]. Therefore, [26] presented a gradient-enhanced fatigue damage model for rubber which is mesh independent. A disadvantage of traditional gradient enhanced damage models is spurious damage growth caused by large deformation in damaged elements.
This results in a large diffused process zone at complete failure [28] 1 . Similar to gradient enhanced damage models, phase-field damage models for fracture ( [6], [21], [23]) treat the sharp discontinuities of cracks in a continuous manner by introducing a finite damage zone, governed by a length scale parameter. Consequently, they are able to handle crack propagation, branching and coalescence. The extension to finite strains and rubber was first published in [22], while [16] presented a rate-dependent phase-field damage model for rubbers.
Gradient-enhanced and phase-field damage models are similar [9], as they both rely on a length scale parameter to produce mesh insensitive results. However, only phase-field damage models do not show the mentioned spurious damage growth, since the damage driving force vanishes for a complete loss of the stiffness [9].
To the best of the authors' knowledge, [3] was the first to consider fatigue damage in a phase-field model by introducing a fatigue history variable in the Ginzburg-Landau equation. [7] additionally introduced damage caused by aging, while [5] defined the internal fatigue history variable with a differential constitutive law to be found. [2] used a different approach and reduced the fracture toughness with cyclic loading. All these contributions show promising results for small strains in one-dimensional settings. Our model is similar to the fatigue model of [3] as we also introduce a fatigue history variable. The differences are, however, that we start the derivation from the balance of mechanical energy, define a fatigue history variable based on the viscous dissipation and allow for finite strains.
The length scale can be interpreted as a material parameter that controls the process zone in which damage occurs. It also affects the material strength in practice, i.e. the stress required to nucleate a crack [30]. [16] has shown that incorrect length scale parameters are identified if the identification is only based on the global mechanical response. To overcome this issue, local strain measurements should be included in the identification process. Recently, [33] proposed a phase-field damage model for linear elastic materials in which the length scale does not affect the global mechanical response, but merely the process zone.
Although the main objective of this work is the presentation of a fatigue phasefield damage model for rubber, this contribution combines the following novelties: • We use the degradation function and crack density function of [33] and show that this extension also leads for finite strains and a visco-hyperelastic material to a global force response, which is not affected by the length scale parameter.
• Applying a cyclic load to the rate-dependent phase-field damage model of [16], we show that this model describes fatigue damage but yields poor accuracy.
• Therefore, we propose an extension of our previous constitutive model [16] so that fatigue damage under cyclic loading can accurately be described.
This framework also works with rate-independent materials.
• We experimentally identify all material parameters: the bulk parameters, the fracture parameters for monotonic loading and the fatigue fracture parameters for cyclic loading.
• Crack nucleation as well as crack propagation can be accurately predicted, which is presented using independent validation experiments. This is, to the best of the authors' knowledge, the first time a phase-field damage model for rubber matches uniaxial tensile test data.
• The fatigue damage model is able to predict Woehler lines and the Paris fatigue crack growth law.
The paper is organized as follows: In Section 2, we derive the generalized fatigue phase-field damage model. Section 3 presents the numerical and experimental results for monotonic loading, while the results for cyclic loading are shown in Section 4. We conclude the contribution in Section 5.
In this work, we denote scalars by lowercase and capital letters (a and A), vectors by bold, lowercase letters (a) and second-order tensors by bold capitals (A).
1 Recent gradient enhanced damage approaches are able to describe sharp damage profiles [28].

Fatigue phase-field damage model
In this section, we derive our fatigue phase-field damage model for a body Ω 0 in the reference configuration with its external boundary denoted by ∂Ω 0 . We start by defining a scalar damage variable d ∈ [0, 1] so that d = 1 on an internal discontinuity Γ 0 and d = 0 on Ω 0 \ Γ 0 (Fig. 1 a)). The motion and deformation of the body are described by displacement u, deformation gradient F = I+∇ 0 u and Green's strain tensor E = 1/2(F T ·F−I).
I denotes the unit tensor and spatial derivatives associated with the reference configuration are denoted by ∂ ·/∂X = ∇ 0 (·). The balance of mechanical energy requires:Ė whereĖ,Ḋ andṖ ext denote the internally stored and dissipated energy and the externally supplied energy per time unit, respectively.
The internally stored energy in the bulk reads: where we have introduced the degradation function g d = g d (d) and the strain energy density ψ bulk . The degradation function controls the mechanical re-sponse with respect to the virgin (i.e. undamaged) state. We use a generalized degradation function, which was recently proposed in [32] and [33]: . (3) The advantage of this generalized degradation function is that the global mechanical response can be made independent of the length scale. For this, the material parameter a must be related to the selected length scale l 0 [33].
Rate-dependent effects are, as in [16], incorporated by splitting the strain energy density ψ bulk into an elastic and viscous contribution: Φ α denotes an internal strain-like tensor, measuring the dissipation in the bulk.
Assuming a material model based on m Maxwell elements, Φ α can be considered as the 3D extension of the 1D strain γ α in a dashpot (see Fig. 2). The first Piola-Kirchhoff stress can consequently be expressed as follows: P ∞ denotes the time-infinity stress and Q α denote the non-equilibrium stresses.
The rate of dissipation is attributed to three dissipative phenomena and can therefore be expressed as 2 : whereḊ visc ,Ḋ crack andḊ crack,visc denote the rate of dissipation due to the viscosity of the bulk, the rate of dissipation due to crack growth and ratedependent crack growth dissipation, respectively.
The energy dissipated by crack growth reads: where G c denotes the energy dissipated by the formation of a unit crack area.
To avoid integrating over the fractured surface Γ 0 and a sharp discontinuity, we approximate Γ 0 ≈ Γ l = Ω0 γ l dV , with the crack density function γ l = γ l (d, l 0 ) [6]. The width of the damage zone is controlled by the length scale l 0 (see Fig.   1 b)). Thus, multiplying Eq. (7) with γ l , we can change the integration to a volume integral over the domain Ω 0 : We furthermore use the following crack density function [32]: where w denotes a geometric crack function characterizing the homogeneous evolution of the phase-field d and c 0 denotes a scaling parameter to ensure Γ 0 = Γ l for l 0 → 0. The geometric crack function w = w(d) ∈ [0, 1] must satisfy w(d = 0) = 0 and w(d = 1) = 1. As in [33], we set w = 2d − d 2 , which yields An extension to a generalized fatigue phase-field damage model is achieved by R can be interpreted as the energetic contribution of an extrinsic volumetric micro-force triggering the fatigue damage growth. Processes on a smaller scale feed this extrinsic micro-force. The relations forĖ,Ḋ andṖ ext are discussed in more detail in Appendix A, whileṘ is presented in the following subsection.

Fatigue damage source
Fatigue damage is assumed to be equivalent to the formation of micro-cracks, which occur for small cyclic loads. These micro-cracks coalesce under cyclic loading and appear as damage d at the macroscale (see Fig. 3). We describe this micro-crack growth with a history variable H, so that the micro-crack formation and coalescence with each cycle is equal to the increase of H. The more micro-cracks have been formed (or the higher the value of H), the more likely macro-crack growth.
R =Ṙ(H) is the energetic fatigue contribution to the macroscopic mechanical energy balance (Eq. 10), linking the micro-cracks growth (H) to the macrocrack growth (d). We assume that micro-cracks only influence the macroscale when they have coalesced to a macro-crack. Therefore,Ṙ only contributes to the energy balance in case of damage growth (ḋ ≥ 0) and we define: ሶ , ሶ The amount of micro-cracks H = f (h) depends on a load history-dependent variable h, which has to be specified. Further, we introduce two fatigue damage material parameters ζ d and ζ e , and a fatigue degradation function g f = g f (d) so that: The fatigue degradation function g f = g f (d) ensures that micro-cracks only form in regions without a macro-crack and therefore must fulfill 4 : In the gradient enhanced damage model of [26], the load history-dependent variable is a function of the strain energy density, i.e. h = f (ψ elas ). The fatigue phase-field damage models of [3] and [5] use h = h(ψ bulk = P :Ḟdt, ), while [2] sets h = f ( ,˙ ). We see micro-crack formation as a dissipative mechanism and therefore propose a load history variable depending on the energy dissipated under cyclic loading in the bulk: i.e. the hysteresis. The advantage of an energetic fatigue failure criterion for rubbers over a maximum stress or strain criterion was shown by several studies , [17]). For simplicity, we set g f = −∂g d /∂d so that: The specific choice of g f has the advantage that h and ψ bulk are equally reduced with growing damage (see Eq. 19). This leads to the most straightforward implementation. Another choice could be proposed to further improve the results but it is not the subject of the current work.

Balance of mechanical energy
Inserting Eq. (6), (15), (A.2) and (A.11) into Eq. (10), we obtain: With Eq. (A.7) and the Neumann boundary conditions: we extract the governing equation for the displacement field u 5 : and for the phase-field damage field d: We want to point out that Eq. (19) contains, in addition to load history h, a rate-dependent driving force: This driving force grows with each load cycle, thanks to its viscoelastic contribution (see Fig. 13). Therefore, fatigue damage can also occur without h (see . Although our original model [16] is able to describe fatigue damage due to the viscosity, we incorporate the history variable h to improve the model's agreement with experimental measurements (compare e.g. the red and blue lines in Fig. 17) and to generalize the framework so that it is applicable to rate-independent models.
Eq. (18) and Eq. (19) are transformed into their respective weak form using the standard Galerkin procedure. We discretize the problem in space with linear, plane stress, isoparametric, quadrilateral elements. Details on the numerical implementation can be found in [16].

Results: Monotonic loading
In this section, we focus on monotonic loading. After presenting our experimental set-ups, we present the identification of the visco-hyperelastic material parameters. Subsequently, we study the influence of the length scale l 0 on the global mechanical response, as well as on the formed damage profile. The following fracture parameter identification considers test configurations in both crack nucleation and propagation. We close the section with the presentation of validation cases.

Experiments
All experiments are performed on an EPDM rubber with a constant temperature of 20 • C. The displacements are measured using a laser extensometer and we use digital image correlation (DIC) to measure local strain fields.

Identification of visco-hyperelastic material parameters
The bulk material parameters are identified using uniaxial tensile tests according to ISO 37 with dumbbell specimens and three clamp velocities (see Tab. 1 for the parameters according to Eq. (A.3) and (A.5) and Fig. 9 for the fit). 6 3.3. Influence of the length scale [33] showed that the combination of the degradation function and crack density function according to Eq. (3) and (9) results in a length scale independent global mechanical response. Since the results in [32] are limited to small strains and linear-elastic materials, we investigate the behavior for large strains and visco-hyperelasticity.

Influence length scale: Uniaxial tension
We perform uniaxial tensile tests (specimen length 20mm and width 4mm, Length scale l 0 = 2mm results in a wider damage zone (see Fig. 4), but the global force to stretch ratio response is nearly identical (see Fig. 5). We conclude that the model of [33] produces also for finite strains and non-linear material models length scale insensitive global mechanical results. 6 Cyclically loading rubber, the load on reloading is less than when loading for the first time [24]. This so-called Mullins effect must be especially considered when modeling fracture and fatigue experiments. The used material parameters (Tab. 1) are a compromise between the response for the first loading and the softened response after multiple load cycles. By using the phenomenological model of [24] to incorporate the Mullins effect, the accuracy of our predictions could be increased. 7 The length scale l 0 is selected large enough (relative to the specimen size) to avoid snapback behavior. Smaller length scales, for the same specimen dimensions, require path-following methods to solve the model. Possible implementation strategies in the context of phase-field damage models can be found in [29].

Influence length scale: Single edge notched tensile test
In this subsection, we present the results for the SENT tests with an initial crack length of 5mm. The dimensions are depicted in Fig. 6 a) In Fig. 6, we compare the finale damage pattern for the two length scales. The damage process zone is wider for an increasing length scale. In contrast, the maximum tensile strength, as depicted in Fig. 7, remains effectively the same.

Identification of fracture material parameters
Single edge notch tensile tests (SENT) (see Fig. 8 a) for the dimensions) are performed to identify, together with the uniaxial tensile test, the phase-field fracture parameters. We measure local strains near the crack tip using DIC, but only for a clamp velocity of 25mm/min due to the specifications of the camera.
The local strain fields are computed with GOM Correlate software. Fig. 8 b) presents the yy-component of Green's strain tensor at a clamp displacement where w = 0.25 is a scalar that controls the influence of the force-displacement data relative to the measured strain fields. Subscript mes denotes experimen-  This approach diverges from the one presented in [16], since the measurements of double edge notch tensile tests are omitted, and uniaxial tensile tests are included. Tab. 2 presents the identified values, which are used for the following computations. The numerical results for uniaxial tension are presented in the next subsection, while the results for crack propagation are shown in Appendix B. For a more detailed investigation of the rate-dependent fracture of rubbery polymers, we refer to our previous publication [16].

Uniaxial tensile test
In this subsection, we show the fitted results of the uniaxial tensile test. The numerically predicted stress-stretch ratio response for the uniaxial tensile test at loading rate˙ = 0.05 1/s is compared to the experimental ones in Fig. 9. The fit between the numerical and experimental results is satisfactory, including the maximum tearing stress and the maximum stretch ratio.

Single edge notched tensile test
We validate the identified fracture material parameters (Tab. 2) with simulations of a SENT test according to ISO 34-1 type C. The geometry and the final damage field is presented in Fig. 10  Applying uniaxial cyclic tension to a slender bar (length 4mm, width 1mm and thickness 0.25mm), we observe hysteresis in the calculated force versus stretch ratio plot (Fig. 12). If we furthermore visualize the development of the energy of the various model components (with parameters a = 50, ζ d = 5 and ζ e = 1), it can be observed that the visco-elastically stored and dissipated energy: increases with every load cycle (Fig. 13). We can deduce from Fig. 13 that the model is thermodynamically consistent since the balance of mechanical energy (Eq. 10) is fulfilled at any given time.

Parametric study: Fatigue damage
Now, we vary the maximum clamp displacement for uniaxial cyclic tension and measure the cycles until failure. The results are presented in Fig. 14 and 15 and so-called Woehler lines are produced. Unloading takes place until 20% of the maximum clamp displacement. Fig. 14 and 15 show that the multiplier (fatigue material parameter ζ d , Eq. (12)) yields a shift of the response, whilst the exponent (fatigue material parameter ζ e , Eq. (12)) also yields a slope difference.
We conclude from the results that the introduction of the history term with two additional fatigue damage material parameters, leads to a general framework which can be used to fit a desired experimental response (see chapter 4.2).
The parameters ζ d and ζ e only affect the result for monotonic loading if both parameters are relatively large (see the response for ζ d = 5 and ζ e = 1 at N = 1 in Fig. 14).

Single edge notched tensile fatigue test
In this subsection, we focus on single edge notched tensile fatigue tests with specimens according to ISO 34-1 type C (see Fig. 16). We apply a cyclic, displacement controlled load until the sample is totally fractured. During each load cycle, maximum stretch ratio λ max is reached and unloading takes place until λ = 1.075. The numerically predicted crack is shown in Fig. 16 for   Figure 17: Fatigue SENT test according to ISO 34-1, type C: Numerically predicted and experimentally measured lifetime. 8 In this contribution we do not consider rate-dependent effects on fatigue damage and perform all cyclic load cases with strain rate˙ = 0.05 1/s. High loading frequencies can lead to generation of heat, resulting in a change of the material properties and thermal degradation of rubbers [10]. Therefore, experiments are normally conducted so that the temperature rise is only a few degrees and negligible. Nevertheless, [35] found that even for a constant temperature, the loading frequency affects the lifetime of non-crystallizing rubbers like EPDM.
[35] observed a decreasing lifetime with a decrease of the frequency and explains this effect with continuous crack growth under static loading. As presented in [16], our model describes this phenomenon and if we increase the frequency of the cyclic loading, we indeed observe a longer lifetime. However, the effect is of a much higher magnitude in our calculations as in [35], which is mainly attributed to a reduced hysteresis of the bulk viscosity. Experiments show a nearly rate-independent hysteresis and [15] proposed to describe this with elasto-plasticity. An extension of the model in terms of plasticity would increase the accuracy of our predictions, but is out of the scope of this publication.

Pure shear fatigue crack growth
Finally, we present predictions for fatigue crack growth in a pure shear setting.
The pure shear specimen is depicted in Fig. 18 and follows ISO 27727. This set-up is characterized by a large width relative to its height. Consequently, the tearing energy G can be calculated independent of the crack length: h 0 = 3mm 36mm Figure 18: Fatigue pure shear test: Specimen dimensions and boundary conditions. The initial crack length is 9mm.
The simulations are displacement controlled with a minimum strain of 20% of the maximum applied strain. 9 By tracking finite elements with d > 0.95, we measure the crack length. In Fig. 19, we plot the current crack length c as a function of the load cycle N for a tearing energy G = 1.283N/mm. As expected from the literature (see for example ISO 27727 or [10]) the initial crack growth rate dc/dN is linear.
Performing the simulation for several tearing energies, we can plot the crack growth rate dc/dN as a function of the tearing energy G (Fig. 20). The crack growth of rubber can be described with a power-law relationship between the crack growth rate and the tearing energy: where a 1 and a 2 denote material parameters.

Computational cost
One disadvantage of phase-field damage models is the high computational cost.
The model was implemented in MATLAB. Because the code is not parallelized, computing on an Intel Xeon processor at 2.4GHz for a pure-shear test case datapoint (17000 nodes) takes about 24h.
The main reason for the high computational costs is the necessity to calculate each load cycle explicitly. We do not apply a cycle jump technique, but we will examine the acceleration of the computations in the future.

Concluding remarks
A rate-dependent fatigue phase-field model for finite strains is proposed based on a load-history dependent fatigue damage source. The general framework allows a flexible choice of the source term and the history variable according to the material being examined. For our rubber, we define the history variable as a function of the dissipation due to the bulk viscosity and introduce two additional fatigue damage material parameters. A reduced polynomial hyperelastic model [34] is used to describe the bulk response, while the rate-dependency of the bulk is incorporated with the well-known material model of [13]. Applying the degradation function and crack density function according to [33], we show global force-displacement results, which are independent of the length scale parameter. All phase-field fracture parameters can be identified from uniaxial The internally stored energy is defined in Eq. 2 with the degradation function g d = g d (d) and the strain energy density ψ bulk . The degradation function controls the stiffness of the bulk and requires: Most phase-field damage models (e.g. [6], [21] and [22]) use a quadratic degra- In the work of [16] multiple degradation functions are assessed in their performance to match experimental crack propagation data of rubber and the quadratic one performs best. We use the generalized degradation function (Eq. 3) from [32] and [33]. Note that for a = 2 the quadratic degradation function is recovered.
The rate of the internally stored energy then reads: Defining the elastic energy according to [34] : where C i are the material parameters, the time-infinity stress reads: The non-equilibrium stresses are calculated according to [13]: where the relaxation times τ α and scalar free energy factors β α are the viscoelastic material parameters. Subscript n denotes converged solutions of the previous time step t n and ∆t is the temporal step size. Further details on the viscoelastic material model and its numerical implementation can be found in [13], [14] and [16].

Appendix A.2. Rate of dissipation
The rate of dissipation is attributed to three dissipative phenomena (Eq. 6).
According to [13], the rate of dissipation due to the bulk viscosity reads: where we have introduced: The energy dissipated by crack growth is defined in Eq. (8). Note that w = d 2 yields c 0 = 2, so that the crack density function as defined in [6], [16] and [21] is recovered. By differentiating with respect to time, the dissipation rate due to crack formation can be written as: The dissipation rate due to the rate-dependency of the crack growth is defined as:Ḋ crack,visc = Ω0 κḋ 2 dV, (A.9) where scalar κ denotes a viscosity parameter [16].
The second law of thermodynamics requires: D =Ḋ visc +Ḋ crack +Ḋ crack,visc ≥ 0. (A.10) Appendix B. Monotonic loading: Results parameter calibration In this appendix, we present results of the fracture parameter identification tests.
The maximum tearing force for clamp velocities of 25mm/min to 200mm/min is successfully predicted (Fig. B.21), as well as the local strain field near the crack tip along paths x 1 and y 1 (see Fig. B.22 and Fig. 8 b) for the definitions of the paths). The correlation between the measured and predicted local strains indicates that the length scale parameter is sufficiently accurately identified.
Note that the calibration is only done with the measurements at 25mm/min and 200mm/min, so that the results at 50 − 100mm/min are true prediction.   Fig. 8 b).

Conflict of Interest and Authorship Conformation Form
Please check the following as appropriate • This manuscript has not been submitted to, nor is under review at, another journal or other publishing venue.
• The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript