Fatigue fracture characterization by cyclic material forces in viscoelastic solids at small strain

The study at hand introduces a new approach to characterize fatigue crack growth in small strain linear viscoelastic solids by configurational mechanics. In this study, Prony series with n-Maxwell elements are used to describe the viscoelastic behavior. As a starting point in this work, the local balance of energy momentum is derived using the free energy density. Moreover, at cyclic loading, the cyclic free energy substitutes the free energy. Using the cyclic free energy, the balance of cyclic energy momentum is obtained. The newly derived balance law at cyclic loading is appropriate for each cycle. In the finite element framework, nodal material forces and cyclic nodal material forces are obtained using the weak and discretized forms of the balance of energy momentum and cyclic energy momentum, respectively. The crack driving force and the cyclic crack driving force are determined by the nodal material forces and the cyclic nodal material forces, respectively. Finally, numerical examples are shown to illustrate path-independence of the domain integrals using material forces and cyclic material forces. The existence of the balance of energy momentum and cyclic energy momentum are also illustrated by numerical examples.


Introduction
A lot of structures or materials exhibit stress relaxation and creep. To describe this behavior, viscoelastic material models are employed. In the literature, many publications present approaches to characterize the behavior of linear viscoelastic solids. The foundation of linear viscoelastic models started with the works of Boltzmann (1878). The stress at a particle within a viscoelastic solid depends on both history and temperature at the particle. Further information regarding thermomechanical analysis of linear viscoelastic solids can be found in the work of Taylor et al. (1970). In the study of Holzapfel and Reiter (1995), a fully coupled thermomechanical finite element formulation of viscoelastic solids is derived. Further studies related to viscoelastic solids can be found in the work of Kaliske and Rothert (1997), Park and Kim (1999), Lakes (1998), Wineman (2009) and in many others. However, in the study at hand, we limit us to small strain linear viscoelastic solids using Prony series with n-Maxwell elements. Analysis of fracture of structures is of high interest for engineers because crack propagation leads to the failure of the structure at loads that are smaller than their real bearing capacity. The analysis of stress fields in linear elastic continua with a sharp crack is car-ried out in the work of Westergaard (1939), Sneddon (1946) and others. Although a lot of studies investigate cracks in elastic continua, however, the first who presented a proper application of the stress intensity factor is Irwin (1957). The stress intensity factor derived by Irwin (1957) is limited to problems with small yielding zones around the crack tip, or, in other words, to cases where the theory of linear elastic fracture mechanics is governing. The analysis of crack propagation in viscoelastic solids started earlier in the works of Schapery (1964), Knauss (1963) and others. The crack speed is defined as a criterion for unstable crack propagation in Knauss (1966). The stability of crack propagation is investigated by Knauss (1969). Furthermore, an energetic fracture criterion to characterize initiation and propagation of cracks in viscoelastic media is derived in the work of Williams (1965). Moreover, an appropriate crack growth law in Maxwell solids subjected to small-scale and large-scale yielding, based on energy and crack opening displacement fracture criteria, is derived in the study of McCartney (1979). A Barenblatt fracture criterion is used in the work of Willis (1967) to characterize crack propagation in viscoelastic solids.
Configurational forces or material forces are known as non-Newtonian forces that act on an inhomogeneity within a homogeneous body. These inhomogeneities can be either inclusions, voids or even cracks. The theory of material forces arises from Eshelby's thought experiment Eshelby (1951). In the aforementioned study, Eshelby considers that the energy variation within a homogeneous elastic body due to the movement of the inclusion is equal to the product between the material force vector, that is acting on the inclusion, and the displacement vector. In the framework of fracture mechanics, the inclusion can be considered as a crack. The material force acting on the crack tip is considered as the crack driving force. In linear elasticity, the component of the material force in crack direction is nothing else but the J -integral derived by Rice (1968). The use of material forces in the field of the finite element method starts with the work of Braun (1997) and Mueller and Maugin (2002). The formulation of material forces in viscoelastic materials can be found in the publications of Kaliske et al. (2005), Näser et al. (2009) and others. In the work Nguyen et al. (2005), material forces are derived for small strain linear viscoelastic solids, that are characterized by the Prony term with one Maxwell element.
At cyclic loading, the pioneering law that describes crack growth is Paris' law introduced by Paris and Erdogan (1963). This law correlates the cyclic stress intensity factor ΔK to the crack growth rate per cycle da d N . However, the cyclic stress intensity factor is limited to the cases, where the theory of linear elastic fracture mechanics (LEFM) is applicable. To overcome the limitation of LEFM, Dowling and Begley (1976) derived the cyclic J -integral ΔJ , which is appropriate for gross plasticity. Moreover, Ochensberger and Kolednik (2014) derived the elasto-plastic cyclic J -integral ΔJ ep using configurational forces, that is appropriate for fatigue crack growth at non-proportional loading. In a further work of Ochensberger and Kolednik (2015), an appropriate parameter known as the active plastic zone elasto-plastic cyclic J -integral ΔJ ep actPZ is derived. This parameter is able to depict the overload effect in elasto-plastic materials. However, this parameter is path-dependent and requires to know in advance the active plastic zone around the crack tip. In our previous study Khodor et al. (2021), a path-independent domain integral based on cyclic material forces is introduced. The domain-integral is able to characterize the overload effect in small strain elastic-plastic materials without any need to identify the active plastic zone.
The aim of the study at hand is to derive a pathindependent domain integral to describe fatigue crack growth in viscoelastic solids at cyclic loading. The method has been previously introduced for elasticplastic materials in Khodor et al. (2021). The approach is based on cyclic material forces. These forces are similar to conventional material forces used at monotonic loading, however, they are derived using the cyclic free energy density.
The paper at hand is structured as follows. In Sect. 2, the constitutive equations used to characterize the viscoelastic solids are presented. Moreover, the balance of energy momentum and the balance cyclic energy momentum for n-Maxwell elements are derived in Sects. 3 and 4, respectively. Furthermore, the numerical computations of nodal material forces and cyclic nodal material forces are shown in Sect. 5. Path-independence of the domain integrals using both material forces and cyclic material forces as well as the proof of the existence of the balance of energy momentum and cyclic energy momentum are presented in Sect. 6. Furthermore, the cyclic material force approach is validated in Sect. 6 by comparing it to experimental data. Finally, concluding remarks are drawn in Sect. 7.

Material description
Small strain linear viscoelasticity in the subsequent derivations is characterized by Prony series. The Cauchy stress tensor is obtained as where τ , e and Δ are the past time, deviatoric strain and volumetric strain, respectively. The parameters G(t) and K (t) are the Prony series shear-relaxation and bulk-relaxation moduli, respectively. The shearrelaxation and the bulk-relaxation moduli are formulated as where τ G i and τ K i are the shear relaxation time and bulk relaxation time, respectively. G 0 and K 0 are the shear and bulk moduli at t = 0. More details regarding the constitutive equations can be seen in the paper of Kaliske and Rothert (1997). n G and n K are the number of shear Maxwell elements and bulk Maxwell elements, respectively. The total number of Maxwell elements is equal to the maximum of n G and n K . α G i and α K i are the relative moduli. Using the relative moduli, α G eq and α K eq can be obtained as The linear viscoelastic material model can be described using a generalized Maxwell solid as shown in Fig. 1. The model consists of a set of n-Maxwell elements in parallel to a spring that is the equilibrium spring. Each Maxwell element consists of a spring and a dashpot, each spring within this element is a non-equilibrium spring. Due to infinitesimal strain, the total total strain tensor ε is an additive decomposition of the elastic strain and the viscous strain within each Maxwell element and can be written as where ε e i and ε v i are the elastic strain tensor in the ith Maxwell element. The total Cauchy stress tensor is a sum of the stress tensors in the equilibrium spring and the non-equilibrium springs and can therefore be obtained as where σ eq and σ neq i are the Cauchy stress tensors in the equilibrium spring and the non-equilibrium spring of the ith Maxwell element. The free energy is also the summation of the free energies in the equilibrium and non-equilibrium springs and can be calculated as where ψ eq and ψ neq i are the contributions to the free energy of the equilibrium spring and non-equilibrium spring in the ith Maxwell element, respectively. Let C eq and C neq i be the fourth order elasticity tensors in the equilibrium spring and the non-equilibrium spring of the ith Maxwell element. In isotropic material, the fourth order elasticity tensor is obtained using Young's modulus E and Poisson's ratio ν or using the bulk modulus K and shear modulus G. Therefore, C eq is obtained using G eq = α G eq G 0 and K eq = α K eq K 0 . On the other hand, G neq i = α G i G 0 and K neq i = α K i K 0 are taken to compute C neq i . By C eq and C neq i , the free energy and the Cauchy stress tensors reformulate into (2.11)

Balance of energy momentum
In Nguyen et al. (2005), the local balance of energy momentum is derived for small strain linear viscoelastic solids with one Maxwell element, in other words n = 1. In this section, the local balance of energy momentum is derived for n Maxwell elements. The inertia effects will be neglected and, therefore, the balance of linear momentum is expressed as where ∇ X · (•) is the divergence operator, and b are body forces. The starting point to obtain the balance is by calculating the gradient of the free energy density, which can be written as where the operator ∇ X (•) is the gradient operator and the term ∂ψ ∂ X exp is the explicit dependence of the free energy on the position X. This term is used when for instance the Young's modulus is a function of the position X, E := E(X). Inserting Eqs. (2.10) and (2.11) into the above equation yields (3.4) Using the additive decomposition of the strain tensor, Eq. (3.4) can be recast into (3.5) By the balance of linear momentum and due to the fact that ε = 1 2 (∇ X u + (∇ X u) ), the term σ : ∇ X ε can be written as where ∇ X u is the gradient of the displacement field and (•) denotes the transpose operator. The insertion of Eq. (3.6) into Eq. (3.5) yields (3.7) The gradient of the free energy can take the form ∇ X ψ = ∇ X ·(ψ I), and, therefore, leads to the balance of energy momentum that can be written as are the configurational volume forces.

Balance of cyclic energy momentum
In case of fatigue crack propagation, a continuum is subjected to cyclic loading, therefore, a new balance of energy momentum is derived in our previous work Khodor et al. (2021) and is called balance of cyclic energy momentum. In the aforementioned work, the balance of cyclic energy momentum is appropriate for elastic-plastic materials at cyclic loading. In the study at hand, this balance law will be derived for small strain linear viscoelastic solids. This balance is a combination of the balance of energy momentum at maximum load and minimum load. The starting point to obtain the balance of cyclic energy momentum is to evaluate the cyclic free energy ψ cycle that is defined as where ψ eq cycle and ψ neq i cycle are the cyclic free energies of the equilibrium spring and the non-equilibrium spring in the ith Maxwell element, respectively. The subscripts (•) max and (•) min designate the value of a specific parameter at maximum and minimum load, respectively. For instance, ε max is defined as the value of the total strain tensor at maximum load and this notation applies for all other quantities. To better understand the calculation of the cyclic energy, the reader is referred to Appendix B. Developing the terms of Eq. (4.1), the cyclic free energy takes the form (4. 2) The gradient of the cyclic free energy can be expressed as is the explicit dependence of the cyclic free energy on the position X. The partial derivatives ∂ψ cycle ∂ε cycle and ∂ψ cycle ∂ε e i cycle can be evaluated as (4.7) The subtraction of Eq. (4.7) from Eq. (4.6) leads to the modified form of linear momentum at cyclic loading, that can be expressed as where σ cycle = σ max − σ min is the cyclic Cauchy stress tensor and b cycle = b max − b min are the cyclic body forces. The cyclic Cauchy stress tensor can be calculated as σ cycle = σ eq cycle + n i=1 σ neq i cycle . The total cyclic strain tensor can be additively decomposed into where ε e i cycle is the cyclic elastic strain and ε v i cycle is the cyclic viscous strain. Using the additive decomposition of the cyclic strain tensor and inserting Eqs. (4.4) and (4.5) in Eq.(4.3) yields (4.9) Defining the cyclic gradient of the displacement field ∇ X u cycle = ∇ X u max − ∇ X u min , the total cyclic strain tensor can be written as ε cycle = 1 2 ∇ X u cycle + ∇ X u cycle . Analogously to σ : ∇ X ε, the term σ cycle : ∇ X ε cycle is computed as (4.10) Knowing that ∇ X ψ cycle = ∇ X · ψ cycle I and inserting Eq. (4.10) into Eq.(4.9) results in (4.11) Introducing the cyclic Eshelby stress tensor cycle = ψ cycle I − ∇ X u cycle σ cycle and the cyclic configurational volume forces , the balance of cyclic energy momentum is defined as It should be noted from the above equation that cycle = max − min and B cycle = B max − B min .

Nodal material forces and cyclic nodal material forces
Consider a cracked body B with traction free crack surfaces, volume V and boundary ∂B, as shown in Fig. 2. A finite element discretization over the whole body domain results in a body denoted by where B e is a finite element with volume V e and N e is the total number of elements in B h . The nodal material forces are equivalent to the internal nodal forces, however, the only difference is that the nodal material forces are obtained using the balance of energy momentum. More details related to nodal material forces can be found in Mueller and Maugin (2002). Carrying out a finite element discretization over the balance of energy momentum, shown in Eq. (3.8), yields the nodal material force where N I are the shape functions, k is the total number of elements connected to the node whereas ksur is the number of surface elements connected to the node. Let F node , F diss node and F sur node be the Eshelby, dissipative and surface parts of the nodal material force, respectively. These parts can be written as where V e is the volume of the element and A e is the area of the element's outer surface. At small strain, the dissipation in B according to Gurtin and Podio-Guidugli (1996) can be written as whereȧ is the crack tip velocity. The second term in Eq. (5.5) is the dissipation due to crack propagation, which leads to the crack driving force According to Runesson et al. (2009) andTillberg et al. (2010), the crack driving force, for the traction free crack surfaces, is obtained from the weak form of the balance of energy momentum in a cracked body and is defined as where δẊ =ȧ. After carrying out a finite element discretization over the whole domain, a discrete form of Eq. (5.7) recasts into (5.8) where N B h is the total number of nodes in B h . Consider the domain P around the crack tip and its discretized form P h , the crack driving force, in case of traction free crack surfaces, can be obtained by carrying out a domain integral over P and its discretized form P h as h is the number of nodes in P h . In the subsequent sections, the crack driving force will be called global material force. The global material force is divided into Eshelby F and F diss dissipative contributions. These contributions are obtained as Analogously to nodal material force, the cyclic nodal material force is obtained from the discrete form of Eq. (4.12) (5.13) Let F node,cycle , F diss node,cycle and F sur node,cycle be the Eshelby, dissipative and surface parts of the cyclic nodal material force, respectively. The contributions of the cyclic nodal material forces can be written as (5.14) In our previous work Khodor et al. (2021), the derived cyclic crack driving force is defined as Using the weak form of balance of cyclic energy momentum in a cracked body, the cyclic crack driving force is computed as The cyclic crack driving force is called cyclic global material force in the subsequent sections. The cyclic global material force can be computed in the domain P h as The path-independence of material forces and cyclic material forces is proved in Appendix A. Let F cycle and F diss cycle be the Eshelby and dissipative components of the cyclic material force, respectively. These contributions can be written as (5.22)

Numerical examples
In this section, several numerical investigations are discussed. At first, an investigation on the pathindependency of the domain integral using the material force approach is carried out for a specimen at monotonic loading. Furthermore, the path-independency of the domain integral obtained using cyclic material forces is investigated for a specimen at cyclic loading. Furthermore, the existence of the balance of energy momentum and cyclic energy momentum is shown. Finally, the cyclic material force approach is validated by experimental results. It should be noted that all the studies are carried out in pure Mode I of fracture. This means that the component of the cyclic global material force perpendicular to the crack vanishes. The simulations are performed using the finite element code ANSYS.
6.1 Path-independency of material forces at monotonic loading In this section, a numerical study is presented to illustrate the path-independent behavior of material forces.
The study is carried out on a two-dimensional precracked disc at plane strain conditions. The specimen used is similar to that in Nguyen et al. (2005) and Özenç  (2014). The disc is subjected to prescribed displacements at its boundaries that are consistent with a Mode I K -field. The applied displacements at the boundary in x-and y-direction are calculated as where κ ∞ = 3 − 4ν ∞ , μ ∞ is the shear modulus at the relaxed state, E ∞ is Young's modulus and ν ∞ is Poisson's ratio at the relaxed state. The material properties at the relaxed state are obtained as

4)
E ∞ = 9K eq G eq 3K eq + G eq , (6.5) ν ∞ = 3K eq − 2G eq 2 3K eq + G eq , (6.6) where K 0 = 16666 Nmm −2 and G 0 = 23076 Nmm −2 . In this numerical example, two Maxwell elements are used to characterize the viscoelastic behavior. The bulk and shear relaxation times as well as the shear and bulk ratio of each Maxwell element are shown in Table 1.
The finite element mesh of the disc is shown in Fig. 3). The disc radius is R disc = 4000h, where h is the element size of the mapped mesh surrounding the crack tip as depicted in Fig. 3). A pre-crack is located at y = 0 and x ≤ 0, with the crack tip located at the centre of the disc at x = 0 and y = 0. In this example, the disc is subjected to an energy release rate G = 10 Nmm −1 , (6.7) Since viscoelastic solids are time-dependent, therefore, the load is applied in the following sequence: at t = 0, K I = 0 and for t ≥ 10 s, K I = G E ∞ . t represents the time. Let F mat 1 , F 1 and F diss be the components in crack direction of the global material force, the Eshelby part of the global material force and the dissipative part of the global material force, respectively. The path-independent feature is illustrated by evaluating the value of the material forces using different domain integrals surrounding the crack tip. The domain integrals are shown in Fig. 4.
In this example and the subsequent examples, 1 is considered as the first contour that includes only the crack tip node. Moreover, the global material force evaluated on the second contour 2 is a summation of the crack tip nodal material force and the nodal material forces of the first set of nodes surrounding the crack tip. Furthermore, the global material forces evaluated at the third contour 3 are equal to the summation of the crack tip nodal material forces and the nodal material It can be seen from Fig. 5 that the material forces are path-independent. It can be noticed that the value of the global material forces at t = 300 s reaches a value of 10 Nmm −1 . The variations with respect to time of the global material force and the dissipative part of the global material force, evaluated at Contour 10, are plotted in Fig. 6.
To better understand this behavior of material forces, the σ yy -component of the total Cauchy stress tensor as well as the σ yy -components of the σ The reason behind the material forces reaching a value almost equal to 10 Nmm −1 is due to the fact that the total Cauchy stress tensor after relaxation is equal to the stress in the equilibrium spring σ eq . Since after relaxation, stress is only in the equilibrium spring, therefore, the model at t = 300 s acts like a linear elastic material model with a Young's modulus E ∞ and a Poisson's ratio ν ∞ . To verify the previous explanation, the same boundary conditions are prescribed to the same disc but using a linear elastic material model instead of the viscoelastic material. The linear elastic material model exhibits E ∞ and ν ∞ as material properties. The global material force value obtained using The purpose of the cyclic material force approach is to yield a path-independent integral, therefore, the pathindependency of cyclic material forces is illustrated by subjecting the same crack disc with the same material parameters represented in Sect. 6.1 to cyclic displacements. The cyclic displacements are applied with a load ratio R load = u min u max = 0.4, where u max and u min are the maximum and minimum applied displacements, respectively. u max is equivalent to a prescribed energy release rate G = 10 Nmm −1 and is therefore calculated using Eqs. (6.1) and (6.2) in x-and y-direction, respectively. In this study, the disc is subjected to 50 cycles. The cyclic global material force is evaluated at  Fig. 10.
Considering the second cycle as shown in Fig. 10, u max is applied at t = 30 s and u min at t = 40 s. Let N cycle be the number of a specific cycle. In terms of N cycle , u max and u min are applied at t = (2N cycle − 1) × 10 whereas u min is applied at t = (2N cycle ) × 10. The evaluation of the cyclic global material force at different contours at N cycle = 2 and N cycle = 50 is shown in Fig. 11. In Fig. 11, the component in crack direction of F mat cycle is denoted by F mat 1 cycle . It can be seen from Fig. 11, that the cyclic material force approach results in a path-independent domain integral. The cyclic nodal material forces and their dis- It should be noted that Eq. (6.8) cannot be used to calculate the cyclic global material force if J ep max and J ep min are replaced by F mat 1 max and F mat 1 min . The reason of Eq. (6.8) being not appropriate to calculate F mat 1 cycle is that F mat 1 min yields negative values. The cyclic material force approach overcomes the issue of negative material forces at unloading stages. Moreover, the cyclic material force results in a path-independent domain integral and, therefore, it can be used as a fatigue crack growth criterion. In linear elastic materials, the cyclic global material force is nothing else as the cyclic Jintegral derived in Tanaka (1983).  Table 1. The lower boundary of the block is fixed in both x-and y-direction, whereas at the upper boundary the x-direction is fixed and a displacement u is prescribed in y-direction. The prescribed displacement is timedependent, where u(t = 0) = 0 and u(t ≥ 10) = 0.1 mm. The nodal material force distribution within the block at t = 10 s and t = 300 s are depicted in Fig.  17. It can be seen from Fig. 17 that the nodal material forces at internal nodes are negligible. Furthermore, the summation of all nodal material forces is zero and, therefore, the balance of energy momentum at monotonic loading is fulfilled.
The aim now is to show the existence of the balance of cyclic of energy momentum, represented by Eq. (4.12). The validation is carried out by subjecting the block to a cyclic displacement with ratio R = u min u max = −1 and the same cyclic pattern as shown in Fig. 10. The cyclic nodal material force distribution at N cycle = 2 and N cycle = 50 are depicted Fig. 18. It can be seen from the sketch, that the cyclic nodal Paris' pioneering law characterizes fatigue crack propagation. This law correlates the rate of crack length per cycle to the cyclic stress intensity factor ΔK as where a is the crack length, C and m are material parameters that are obtained experimentally. ΔK is however limited to the case of LEFM, therefore, Dowling and Begley (1976) introduced the cyclic J -integral for gross plasticity. In viscoelastic fracture mechanics, an analytical cyclic J -integral is derived in Kuai et al. (2009) andLee et al. (2015) in order to model fatigue crack growth in asphalt concretes. The aforementioned cyclic J -integral is denoted in this study by ΔJ visc . ΔJ visc is the difference between the generalized Jintegral derived by Schapery (1984) at maximum and minimum load. ΔJ visc is derived using the cyclic stress intensity factor and the creep compliance. Moreover, it requires in advance to have the cyclic stress intensity factor as an explicit function of the crack length a.
In a complex geometry, ΔK cannot be obtained as an explicit function of a and, therefore, ΔJ visc is not valid for the case of complex geometries. The purpose of this example is to illustrate that the cyclic global material force can be used as a fatigue crack growth criterion, leading to the new form of Paris' law where C and m are material parameters. To validate that Eq. (6.10) exists, a numerical study is conducted on the same specimen as presented in the work of Kuai et al. (2009). In the aforementioned study, a fatigue crack growth test is carried out on an asphalt concrete specimen. Asphalt concrete is a very complex material that exhibits relaxation and creep. A simplified description of the asphalt concrete is obtained using a small strain viscoelastic material model. The specimen has a thickness t = 50 mm with an initial crack length a 0 = 30 mm and width W = 110 mm. The specimen is subjected to a periodic load P that is defined as where P 0 is the load amplitude and f is the load frequency. Numerically, only half of the specimen is used due to symmetric boundary conditions. The specimen dimensions and finite element mesh are shown in Fig. 19. A mapped mesh with an element size of 1 mm × 1 mm along the crack path is used. The specimen is simulated at plane stress conditions with 50 mm thickness. One volumetric Maxwell element and one isochoric Maxwell element with α G 1 = 0.4, τ G 1 = 0.02 s, α G 1 = 0.5 and τ K 1 = 0.05 s have been assigned to the specimen. The initial shear modulus G 0 and initial bulk modulus K 0 are obtained from the initial Young's modulus E = 2000 Nmm − 2 and initial Poisson's ratio ν = 0.3. The cyclic global material forces are evaluated at the end of each cycle. From Eq. (6.11), it can be noticed that the load reaches its maximum value P max and its minimum value P min at t max = 2N cycle −1 2 f and t min = N cycle f , respectively. The cycle ends at t min , therefore, the quantities needed to evaluate the cyclic global material force are obtained from t max and t min of each cycle.
Viscoelastic materials are time-dependent which means that the loading frequency or rate have an impact on the behavior of the material. Different material behavior stands for different states of stresses and strains and, therefore, different values of the cyclic global material forces. Therefore, the starting point in this study is to evaluate the effect of the load frequency and the load amplitude on the cyclic global material force. The behavior of the cyclic global material force is investigated by applying the load for 50 cycles at different loading frequencies and at a fixed loading amplitude P 0 = 200 N.
The variation of F mat 1 cycle at different frequencies is depicted in Fig. 20. It can be seen from this figure that the value of F mat 1 cycle decreases with increasing the frequencies. This means that for higher frequency the fatigue life of the material increases. In the experimental work of Lee et al. (2015), it has been shown that the fatigue life of asphalt concrete increases with increasing the loading frequency. Therefore, it can be concluded that the cyclic global material force takes into consideration the effects of loading frequency.  Moreover, the effect of the loading amplitude on the cyclic global material force is investigated by varying the loading amplitude and fixing the loading frequency to f = 10 Hz. The results of this investigation are plotted in Fig. 21. It can be seen from the obtained results that the cyclic global material force values increase with increasing the loading amplitude. Higher values of F mat 1 cycle means that the fatigue life of the specimen is shorter. It is also observed in the work of Lee et al. (2015), that higher loading amplitudes resulted in shorter fatigue life. Therefore, the cyclic global material force characterizes the effect of the loading amplitudes.
It  Crack length values at different cycles obtained from experimental data and best-fit function fifth cycle. An additional study is carried out in order to review the variation of the cyclic global material force with crack propagation. The specimen is subjected to the same periodic load P 0 = 200 N at f = 10 Hz. The crack propagation is achieved by applying a node release algorithm with a crack increment Δa = 1 mm.
The crack propagates at the end of the fifth cycle of each crack tip. The reason behind releasing the node at the end of the fifth cycle is that the cyclic global material force reaches a constant values after 5 cycles as shown in Fig. 21. This behavior is also observed in the previous example, see Fig. 14. Moreover, the values of the cyclic global material force are plotted against the experimental crack growth rate da d N cycle . da d N cycle is calculated using the crack length at different cycles obtained from the experimental work of Kuai et al. (2009). The experimental crack length increase with respect to cycles as well as the best-fit function are plotted in Fig. 22.
The best-fit function of the crack length and the crack crack growth rate per cycle are defined as a(N cycle ) = 10 −16 N 3 cycle − 6 × 10 −11 N 2 cycle + 2 × 10 −5 N cycle + 29.5, (6.12) da d N cycle = 3 × 10 −16 N 2 cycle − 12 × 10 −11 N cycle + 2 × 10 −5 . (6.13) The variation of F mat 1 cycle is depicted in Fig. 23 Using Eq. (6.14), the Paris' parameters of Eq. (6.10) are obtained, where C = 3.21 × 10 −5 and m = 0.55. Therefore, it can be concluded that the cyclic global material force can be used as fatigue crack growth parameter in Paris' law instead of the cyclic J -integral and cyclic stress intensity factor.

Conclusions
The focus of the study at hand is mainly on small strain viscoelastic solids, that are characterized by Prony series with n-Maxwell elements. The main aim of this paper is to present a path-independent domain integral using cyclic material forces in viscoelastic solids to evaluate the cyclic energy release rate. As a starting point, the local balance of energy momentum for n-Maxwell is derived using the free energy density. Furthermore, at cyclic loading, the free energy density is replaced by the cyclic free energy of a cycle. The cyclic free energy is computed from the cyclic stresses and the cyclic strains. The balance of cyclic energy momentum is obtained using the gradient of the cyclic free energy.
In the framework of numerical methods, a finite element discretization is carried out leading to nodal material forces and cyclic nodal material forces in viscoelastic solids. The crack driving force at monotonic loading is obtained by evaluating a domain integral around the crack tip. Numerically, the crack driving force is computed by summing up all the nodal material forces within a region surrounding the crack tip. Analogously to the crack driving force at monotonic load, the cyclic crack driving force is evaluated numerically by performing a summation of cyclic nodal material forces within a region surrounding the crack tip.
Several numerical studies are carried out in this research. At first, it has been shown that the domain integral obtained using material forces at monotonic loading is path-independent, in other words the results obtained using two different domains are equal. Moreover, it has been shown that after relaxation, the value of the global material force in crack direction is equal to the prescribed energy release rate on a linear elastic material having the relaxed Young's modulus E ∞ and Poisson's ratio ν ∞ as material parameters. Furthermore, path-independency of the domain integral using cyclic material forces is validated. The cyclic global material force overcomes the problem of negative values obtained using the global material force at minimum load. Moreover, the existence of the balance of energy momentum and cyclic energy momentum is shown. Furthermore, the cyclic material force approach is validated using experimental results. It has been shown, that the cyclic global material force characterizes the effects of both the loading frequency and the loading amplitude. Finally, the relation between the crack growth rate per cycle and the cyclic global material force illustrates that the cyclic material force approach can characterize fatigue crack growth.
Acknowledgements The authors gratefully acknowledge the financial support of ANSYS Inc., Canonsburg, PA 15317, USA.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Code availability The data are generated using the finite element code ANSYS and the cyclic global material forces are computed using an in-house post-processing subroutine.

Material availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendix A: Path-independent domain integral
In this section, the path-independence of material forces and cyclic material forces is validated. Consider the cracked body shown in Fig. 2, the path-independent J -integral derived by Kishimoto et al. (1980a) for small strain elastic-plastic materials is defined aŝ where ε p is the plastic strain tensor. For more details related to the proof of path-independency of the contour integral, the readers are referred to the work of Kishimoto et al. (1980a). According to Shih et al. (1986), the surface integral in Eq. (A.1) can be transformed into domain integral asĴ = − P ∇ X · dV + P σ : ε p − (∇ X u) b dV.

(A.2)
In visco-elastic materials, the crack driving force can be evaluated 3) an additional term is accounted for the explicit dependence of the free energy density on the position within the body. Therefore, it can be concluded that the crack driving force computed by the material force method is path-independent. In the work of Kishimoto et al. (1980b), the field quantities used to evaluate thê J -integral are considered to change only due to crack propagation, in other words the crack length is used as time-like variable. Due to the aforementioned assumption, theĴ -integral is proven to be path-independent. At cyclic loading, the cyclic crack driving force is evaluated at each cycle. Furthermore, the cyclic stresses, strains and energies needed to compute the cyclic crack driving force are considered to be constant within one cycle, if the crack does not propagate. This statement means that the cyclic quantities at one cycle depend only on the crack length and, therefore, the crack length can be used as a time-like variable for the cyclic quantities. Using the crack length as a time-like variable for the cyclic quantities, a path-independent domain integral for evaluation of the cyclic crack driving force is derived as Using Eq. (A.6), it can be concluded that the integrals evaluated at any two arbitrary domains are equal. This means, the domain integral is path-independent. At cyclic loading, the same mathematical proof is repeated, however, all the quantities are replaced by their cyclic counter parts.