Accounting for viscoelastic effects in a multiscale fatigue model for the degradation of the dynamic stiffness of short-fiber reinforced thermoplastics

Under fatigue loading, the stiffness decrease in short-fiber reinforced polymers reflects the gradual degradation of the material. Thus, both measuring and modeling this stiffness is critical to investigate and understand the entire fatigue process. Besides evolving damage, viscoelastic effects within the polymer influence the measured dynamic stiffness. In this paper, we study the influence of a linear viscoelastic material model for the matrix on the obtained dynamic stiffness and extend an elastic multiscale fatigue-damage model to viscoelasticity. Our contribution is two-fold. First, we revisit the complex-valued elastic models known in the literature to predict the asymptotic periodic orbit of a viscoelastic material. For small phase shifts in an isotropic linear viscoelastic material, we show through numerical experiments that a real-valued computation of an “elastic” material is sufficient to approximate the dynamic stiffness of a microstructure with a generalized Maxwell material and equal Poisson’s ratios in every element as matrix, reinforced by elastic inclusions. This makes standard solvers applicable to fiber-reinforced thermoplastics. Secondly, we propose a viscoelastic fatigue-damage model for the thermoplastic matrix based on decoupling of the time scales where viscoelastic and fatigue-damage effects manifest. We demonstrate the capability of the multiscale model to predict the dynamic stiffness evolution under fatigue loading of short-fiber reinforced polybutylene terephthalate (PBT) by a validation with experimental results.


State of the art
Different aspects of the material behavior of short-fiber reinforced thermoplastics under cyclic loading for a high number of cycles have been studied by many authors. Their matrix materials, i.e., crystalline polymers are known to exhibit viscoelastic material behavior due to time-dependent chain oscillations or grain boundary relaxation at the molecular scale [1]. These molecular motions can be thermally activated. Thus, polymers show a strong temperature dependence in their viscoelastic properties. In fact, in linear viscoelastic materials, an increase in temperature corresponds to a shift of the material properties, i.e., has a similar effect as a decrease in the applied loading rate or frequency. The empirical model proposed by Williams, Landel and Ferry [2] is the standard model to represent temperature-rate superposition, see, e.g., Maurel-Pantel et al. [3] and Dorléans et al. [4] for recent works on thermoplastics. Additionally, the viscoelastic properties may be influenced by the fillers [5,6] or blending of different thermoplastics [7]. At low stress levels, the polymers are well-described by linear viscoelasticity, like generalized Maxwell models, which were applied to, e.g., polyamide [8] or polystyrene [9]. At higher stress levels, nonlinear viscoelastic effects become apparent in some materials, e.g., polypropylene [10].
The viscoelastic nature of the thermoplastic matrix is reflected in the macroscopic properties of polymer-based composite materials and were reported in various fiber reinforced polymers, e.g., polybutylene terephthalate (PBT) [11] or natural fiber-reinforced polyethylene [12]. Abdo et al. [13] studied frequency-dependent loss and storage moduli of fiber reinforced PBT. Stadler et al. [14] investigated the life time of short fiber reinforced polymers under combined cyclic and creep load.
Under cyclic loading, thermoplastic materials exhibit fatigue damage. In composite materials under fatigue loading, the evolving damage in the matrix translates to an anisotropic stiffness loss of the material. This stiffness decrease of the material is a measurable quantity for characterizing its progressive degradation. A number of authors [15][16][17][18] studied the dynamic stiffness of the material to obtain insights into the fatigue processes in both neat and reinforced polymer materials.
For viscoelasticity, phenomenological models include power-law models developed for, e.g., glass-fiber reinforced polypropylene [19], poly-propylene-wood composites [20] and sheet molding compounds (SMC) [21], or the HRZ [22]model developed for, e.g., short-fiber reinforced materials and also applied to long-fiber composites [23], describe the creep behavior accurately. Dacol et al. [24] recently discussed an approach to transfer the relaxation properties analyzed in creep experiments to the dynamic behavior, valid for pure and reinforced polymers. Phenomenological progressive damage or fracture models, directly formulated for the composite material [25][26][27][28][29][30][31][32], may be used to model and predict the material behavior or the material degradation. Models based on stress invariants, e.g., Puck's criteria [33], were used to predict failure [34][35][36] of continuous fiber reinforced materials. Another approach to model the composite's material behavior starts from the material behavior of the constituents and predicts the resulting effective properties of the composite by homogenization techniques. In addition to analytic techniques such as the mean-field approaches, computational homogenization techniques may be used. In the context of viscoelastic material models, several authors used mean field approaches to predict the relaxation [37,38] and cyclic behavior of composites [39]. Computational methods were developed by, e.g., Tran et al. [40] and Curtois et al. [41,42]. In the context of damage, several authors [43][44][45][46] used homogenization techniques to predict the effective stiffness degradation of the material based on the evolution of its constituents.
Under constant loading, viscoelastic materials exhibit a creep or relaxation behavior reaching a steady state at infinity. In the context of fiber reinforced structures, the relaxation behavior of the constituents translates to the composite's behavior. The Laplace-Carlson transform (LCT) is often used to transfer techniques developed for linear elasticity to lin-ear viscoelastic materials [47][48][49][50]. Schapery's collocation method [51,52] can be used to find approximations of the composite's relaxation function.
Similarly, under constant cyclic loading, viscoelastic materials exhibit a relaxation to a cyclic steady state, which can be described as a linear elastic material with a complexvalued stiffness tensor [53,54]. Gallician et al. [55,56] extensively studied the asymptotic relations of linear and fractional viscoelastic materials. Figliuzzi et al. [57] applied the scheme to carbon black filled rubber materials, using a Maxwell model for the matrix material.
For neat thermoplastic materials, Zerbe et al. [58] as well as Praud et al. [59] proposed viscoelastic-viscoplasticdamage models. Coworkers of Praud et al. [59] later extended the model to include thermo-mechanical effects [60]. For nanoparticle-filled epoxy, Arash et al. [61] proposed a viscoelastic damage model. He et al. [62] discussed a phenomenological viscoelastic-viscoplastic-damage material model for fiber reinforced polymers. Recently, Chen et al. [63] developed a rate-dependent material model for short glass-fiber reinforced materials taking into account viscoplasticity, interphase effects and matrix damage at the microscale, with a Mori-Tanaka scheme.
In the context of fatigue, Guedes [64] investigated the life time of short-fiber reinforced polymers under creep and fatigue loading. Hammoud et al. [65] suggested a viscoelastic damage model for neat polymers. Despringe et al. [66] applied the Mori-Tanaka scheme to viscoelastic damage models in short-fiber reinforced polymers.

Contributions
In contrast to metals, polymer-based materials exhibit stiffness degradation under fatigue loading. This observation is the starting point of many damage materials in the fatigue literature [67][68][69]. They aim to model the material's degradation under cyclic loading based on measurements of the dynamic stiffness. At the same time, it is well-accepted in material science that polymers are viscoelastic in nature [70]. Thus, the observed dynamic stiffness of polymer-based materials depends on the loading frequency [13,71].
In the work at hand, we extend existing computational multiscale approaches for elastic fatigue-damage to viscoelasticity. We choose a linear viscoelastic material model and study its influence on the dynamic stiffness under cyclic loading conditions. Subsequently, we incorporate the most relevant effect, i.e., the frequency dependence at high cycles, into a fatigue-damage material model.
In Sect. 2, we start by a systematical investigation of the influence of viscoelasticity on the dynamic modulus and its evolution in experiments. Under cyclic loading, we observe that the material reaches an asymptotic periodic orbit within a fairly low number of cycles. Subsequently, we introduce a linear viscoelastic material model to describe the observed effects. For a linear viscoelastic material, the shape of the hysteresis curve at high cycles depends on the applied frequency, but does not evolve further. In contrast, fatigue-damage evolves at all stages of high-cycle fatigue. We thus propose to treat the two effects separately and utilize a two-step analysis of the current dynamic stiffness that is both accurate and efficient. In a first step, we discuss how to compute the dynamic stiffness for the periodic orbit, i.e., the asymptotic dynamic stiffness of the undamaged material (Sect. 3). In a second step, we integrate the obtained fixed-point stiffness into a cycle-based fatigue-damage model (Sect. 4).
In Sect. 3, we discuss an elastic model with complexvalued stiffness [53,54] to compute the periodic orbit for cyclic loading of viscoelastic materials. We show that for components consisting of linear elastic materials and just one viscoelastic constituent with a small phase shift, a real-valued scheme is sufficient to approximate the materials dynamic stiffness. Thus, for short-fiber reinforced thermoplastics with elastic fibers, a real-valued scheme is applicable.
In Sect. 4, based on the derived real-valued elastic scheme, we propose an approach to integrate the frequencydependence of the dynamic stiffness into existing fatiguedamage models. The method modifies the undamaged stiffness according to the current applied frequency. Apart from a simple precomputation of the effective matrix properties dependent on the current frequency, there is no need for any increase in the computational effort when going from an elastic-damage to a viscoelastic-damage model. We demonstrate the capability of the model to reproduce experimental results in a validation step.

In fatigue experiments
Polymer materials degrade under cyclic loading. To characterize the evolving fatigue damage in experiments, multiple quantities may be monitored. Apart from self-heating [72,73], the most common approach is an analysis of the strainstress hysteresis at the current cycle. As scalar measures of the hysteresis curves, stiffness measures are widely used in literature. They may be chosen in different ways, e.g., the slope of the curve close to the point of reversal. One of the most common choices is the dynamic modulus, defined by where σ max and σ min are the maximum and minimum stress level reached within the current hysteresis cycle, respectively. With ε max and ε min , we refer to the maximum and mini-mum strain level within the same cycle. To understand the evolution of this quantity in polymer-based materials, we first investigate the evolution of the stress-strain hysteresis, which serves as the basis for the evaluation of the dynamic stiffness in experiments. Therefore, we subject a short-fiber reinforced PBT specimen to cyclic loading at constant stress amplitude and a stress ratio of R = 0, where The hysteresis curves and their evolution in the experiment are shown in Fig. 1. For the evolution of the hysteresis under fatigue loading, we observe two main effects: a shift of the hysteresis to larger strain values (creep) and a rotation of the hysteresis in the stress plane (damage). The shift effect is especially pronounced in the first few cycles: in Fig. 1a, we observe a visible shift of the hysteresis between cycle number 17 and 18. For higher cycle numbers, e.g., cycle number 27 and 28, the shift is barely noticeable. Additionally, we observe that the shift from 17 to 27 is of similar magnitude as the shift from 27 to 67 cycles. Thus, the shift effect decreases with increasing number of cycles. The second observation is a rotation of the hysteresis. In Fig. 1b, the hysteresis curves at cycle 140 and 9677 are compared. The dashed lines in gray and black indicate the slope of the two hysteresis curves. To simplify a comparison on a visual basis, we shifted both curves to zero. We observe that the slope of the hysteresis loops decreased under fatigue loading, i.e., the hysteresis loop was rotated clock-wise.
Another effect which we observe in fatigue experiments on short-fiber reinforced PBT is its frequency dependence. To realize the fatigue experiments in the shortest possible time, experiments at different loading amplitudes are often performed at different frequencies [74,75] to avoid self-heating of the material at high loading amplitudes. The choice of a different frequency also influences the measured dynamic modulus. Goméz et al. [71] published experiments on short-fiber reinforced PBT, in which they characterized the dynamic stiffness at different frequencies. The dynamic modulus of PBT with 20 weight-% of reinforcements at 50.2 • C is shown in Fig. 2a. A clear increase of the material's dynamic modulus with increasing frequency is observed. Similarly, for short-fiber reinforced PBT at room temperature, we observe a dependence of the dynamic modulus on the loading frequency in fatigue experiments, see Fig. 2b. The specimens used for fatigue testing were cut from injection molded plates parallel (0 • ) and perpendicular (90 • ) to the flow direction of the mold. For each of the load cases listed in Table 1, three specimens were tested. The semi-transparent bands around the data plotted in Fig. 2b mark the standard deviation of the experimental results. We observe, both for the 0 • -oriented and the 90 • -oriented specimens, that the load cases at higher  frequencies exhibit a higher dynamic modulus. This is in accordance with the results from Gomez et al. [71]. The evaluated dynamic modulus is affected by all of the observed effects, i.e., shift of the strain hysteresis to higher strain values, rotation and frequency dependence. The effects originate from two material properties: viscoelasticity and damage. Viscoelastic behavior of the material is observed in the shift (creep) of the stress-strain hysteresis to higher strain values as well as the frequency dependence of dynamics modulus. Damage behavior is reflected in the clock-wise rotation of the stress-strain hysteresis with increasing load cycle. We investigate the consequences of a continuum mechanical model for linear viscoelasticity on the evaluated dynamic modulus in Sect. 2.2.

In viscoelastic materials
To investigate the effect of a linear viscoelastic material model on the measured dynamic modulus, we assume that the material follows the behavior of a generalized Maxwell model [76,77]. Its rheological model is schematically shown in Fig. 3. The set of internal variables in this model consists of N viscoelastic strains ε v, j . In the context of generalized standard materials [78,79], the generalized Maxwell model is defined via the free energy density where C j denotes the stiffness of the j-th Maxwell element, and the dissipation potential Here, D j refers to the viscosity tensor of the j-th Maxwell element. We assume that both the stiffnesses and the viscosities are isotropic and their Poisson's ratios coincide in every element. Then, for a given function and a Young's modulus E j , viscosity η j and Poisson's ratio ν j in every Maxwell element, we define the viscosity tensor and the stiffness tensor where I is the second-order identity tensor and I s denotes the symmetric part of the fourth-order identity tensor. Additionally, we assume that ν j = ν for all j ∈ {0, . . . , N }. (2.8) For concise notation, we define the time constant τ j and its inverse β j given by for every Maxwell element.
To understand the evolution of the dynamic modulus for a generalized Maxwell model, we first investigate the onedimensional case. We choose to use an example with N = 1 and the material constants E 0 = 1 GPa, E 1 = 10 GPa and η 1 = 1 GPa·s. We subject the material to loading at a constant stress amplitude of σ a = 30 MPa, a stress ratio of R = 0 and a frequency of f = 1 Hz. The results are shown in Fig. 4. In the evolution of the stress-strain curve in Fig. 4a, a shift of the hysteresis curve to higher strains is visible. As observed in the experiments discussed in Sect. 2.1, the largest shift manifests in the first cycle and the magnitude of the shift decreases steadily with the number of cycles. When evaluating the dynamic modulus at every cycle, see Fig. 4b, we observe a strong increase within the first few cycles. For cycles larger than ten, the material response reaches a steady state, i.e., the dynamic modulus approaches an asymptotic value. This observation is independent of the applied load amplitude and stress ratio. Indeed, in Fig. 5a, we show the material subjected to cyclic extension at different stress amplitudes. We observe that the stress-strain curve approaches a periodic orbit, i.e., the difference between the hysteresis curves of two subsequent cycles become smaller with every cycle. In the sixth cycle, we compare the dynamic modulus marked by the dashed line in the plots at every stress amplitude. We observe that the slope of the dashed line, i.e., the dynamic modulus, is independent of the applied stress amplitude. Moreover, we plot stress-strain curves at different stress ratios R in Fig. 5b. Similar to the hysteresis evolution under different load amplitudes, the hysteresis loops converge to periodic orbits, and the evaluated dynamic modulus in the sixth cycle is independent of the applied stress ratio.
In contrast to stress amplitude and stress ratio, the steady state reached for a viscoelastic material is not independent of the applied frequency. In Fig. 6, the material is subjected to cyclic loading at different frequencies. In Fig. 6a, the stress-strain curves of the material for the first four cycles are plotted. The hysteresis evolves around a center point in the stress-strain diagram which is independent of the frequency. However, the shape of the hysteresis and, in particular, its symmetry axis, depends strongly on the frequency. The dynamic moduli are plotted in Fig. 6b. They reflect the rotation of the symmetry axis of the hysteresis curve at different frequencies. The modulus values rapidly converge to an asymptotic value. The reached asymptotic value depends on the frequency.

Material model
The constitutive behavior of heterogeneous materials is strongly influenced by the properties of the underlying microstructure. For short-fiber reinforced materials, these microstructures are characterized by the fiber orientation [80,81], the fiber volume content [71] and the aspect ratio [82,83] of the fiber inclusions. Thus, accurate material models need to account for the characteristics of the underlying microstructure. Due to the inherent anisotropy, characterizing ad-hoc phenomenological material models requires an expensive experimental program, in particular for long-term loading. In contrast, computational multiscale methods offer a possibility to derive the macroscopic material model from the material behavior of the constituents by computational means. Indeed, characterizing the individual phases typically We thus chose to model the material behavior of the constituents, i.e., matrix and fiber, and to derive the macroscopic material behavior by micromechanics.
We consider a cubic cell Y ⊆ R 3 , together with spatially varying positive definite stiffness tensors C m (x) for m ∈ {0, . . . , N }, and the coefficients β j (x) with j ∈ {1, . . . , N } are given at every microscopic point x ∈ Y . We seek the strain fields ε and ε v, j , the stress field σ and the displacement fluctuation field u satisfying together with the initial and boundary conditions u periodic on ∂Y ,  3. In formulas (3.1)-(3.6), we suppress the dependence of the fields on x ∈ Y and the time t for simplicity of notation. For later use, we record that for an initially stress free material, i.e., σ (x, 0) = 0, the differential form of the constitutive law (3.2) and (3.3) can be equivalently represented in integral form [84] σ with the material function

Choice of parameters
Similar to Krairi et al. [43], to identify matrix material parameters valid for the composite's material behavior, we consider tensile tests at different strain-rates. For short-fiber reinforced PBT, Mortazavian and Fatemi [85] reported on such experiments. However, these authors did not characterize the fiber orientations within the specimen. We assume the fiber orientation inside the injection molded specimen to be layered with the fiber orientation tensor components shown in Fig. 8a, as suggested by the analysis of CT-scans of similar materials [86][87][88][89][90]. We use a commercial PBT with 35 weight-% of fiber inclusions and a respective fiber volume content of 17.8% [91]. The aspect ratio of the fibers is set to 29, see, e.g., Hessman et al. [89] for an analysis of typical aspect ratios. Based on the strain-rate dependent experiments on the composite published by Mortazavian and Fatemi [85], we identify the moduli E j and relaxation time constants τ j = 1/β j of the polymer material. The authors used the ISO-527 standard test method [92] to obtain the Young's modulus of the tested specimen from its stress-strain curve. In this standard, the secant between a strain value of 0.05% and of 0.25% is used to identify the Young's modulus of the specimen. The experiments were reproduced computationally and the Young's modulus of the material was identified from the computed stress strain curve in the same manner, as shown in Fig. 7. The computed stress-strain curves for the parameters given in Table 2 are shown using dashed lines. The evaluated Young's modulus in the direction of loading, both for the 90 • -and the 0 • -oriented specimen, is indicated by the solid line, i.e., the secant on the computationally generated stressstrain curve. In Fig. 7, the Young's moduli obtained from the computationally generated stress-strain curves and from experimental measurements are compared. We observe that for both the 0 • -oriented specimen and the 90 • -oriented specimen, the computational prediction of the material modulus matches the experimental data rather well. Thus, the dependence on both strain rate and orientation is well-reproduced by the linear viscoelastic model with material parameters given in Table 2. We fix the identified parameters for PBT for the remainder of the manuscript, unless stated otherwise (Fig. 8).

Dynamic stiffness of linear viscoelastic composites for long-term cyclic loading
We are interested in the long-term response of the viscoelastic composite subjected to a periodic excitation. Thus, we consider a periodic loading σ (t) of the system (3.1)-(3.6) and seek solutions which are periodic in time with a period T = 2π/ω. The system (3.1)-(3.6) possesses a unique solution, and we are concerned with obtaining computationally feasible approximations. For this purpose, we expand the total strain and the viscoelastic strain fields into Fourier series w.r.t. the time variable.
For the time derivative, we obtaiṅ The evolution equations (3.3) for the viscoelastic strains can be reassembled as (3.14) As differential equations in real space get transformed to algebraic equations in Fourier space, we may solve for the fieldsε v, j With this formula at hand, we may recast the constitutive law (3.2) in a form reminiscent of linear elasticity, but with complex coefficientŝ Suppose that the tensor fields C j are positive definite. Then, provided the Fourier coefficientsσ k of the macroscopic stress loading are contained in a specific subset F of the integers, i.e., we concludeû Put differently, the displacement fluctuation field is then supported on this reduced set F of Fourier frequencies, as well.
We obtained a complex-valued elastic model [53] for the viscoelastic material's asymptotic periodic orbit. Typically, fiber reinforced polymer materials comprise just one viscoelastic constituent, the matrix, and one type of elastic inclusions, the fibers. We will discuss this special case in the following and derive a real-valued approximation for the dynamic stiffness of such materials.
To simplify further discussion, we assume a stress-driven sinusoidal load in the form with a stress-tensor valued "amplitude" σ a .
With the help of the expressions for the cosine and the sine in terms of the exponential function cos y = e iy + e −iy 2 and sin y = e iy − e −iy 2i , y ∈ R, (3.21) this loading may be written in the form Thus, the Fourier coefficients compute aŝ (3.23) By the previous argument (3.19), we thus seek a strain field ε supported only on the ±1 Fourier frequencies. Although complex-valued expressions are typically manipulated more easily, physical considerations enforce us to work with realvalued solutions for the macroscopic stress loading (3.20). For this purpose, it is more convenient to seek real-valued are valid. We will also use a similar representation for the local stress fieldŝ In the fiber material, the linear elastic constitutive law (3.16) gives within the linear viscoelastic matrix. With the abbreviations we deduce the strain-explicit representation need to be determined via computational homogenization. Due to the tensorial nature of the strain and the stress tensor, the convenient splitting of a superimposed cosine and sine wave in terms of a shift angle and an amplitude, well-known in the scalar case, is not available in this general setting. However, it appears plausible-at least in the case q/r 1-that we may approximate the effective response 3.35 by an effective model with shiftangle-amplitude form where C is a scalar and the corresponding sine-cosine form is Then, the relations hold. To approach this problem in a systematic manner, we use approximations of the local compliance tensors (3.30) that lead to a homogenized response of the appropriate form (3.36). Subsequently, we will use computational experiments to determine which approximation is appropriate. We consider three approximations of this stress-strain relation for small phase shifts in the viscoelastic material, i.e., for q/r 1.
1. Matrix-dominated approximation MA: For small fraction q/r , we consider the approximation (3.39) We write the approximated stress-strain relation as We define the approximated dynamic compliance as The maximum effective strain is obtained at We define δ fiber = 0 (3. 48) and the approximated dynamic compliance as (3.50) 3. Limit shift approximation LA: In the limit case of a neat matrix material or a neat fiber material under cyclic loading, the above approximations violate the expected phase shift of the homogenized material. Indeed, using MA, we obtain for neat fiber material. For FA, we obtain

Summary of the surrogate elastic computation (SEC) approach
Assumptions: • single isotropic viscoelastic material, reinforced by elastic fibers • Poisson's ratio ν equal in all stiffnesses and viscosities of the viscoelastic material Input: elastic stiffness C f ; viscoelastic properties of the matrix E i , β i for all i ∈ {1, . . . , N }; Poisson's ratio ν; angular frequency ω of the loading Computation: 1. Evaluate the dynamic material constants of the viscoelastic material for neat matrix material. To circumvent this, we propose to use a combination of the above approximations (3.53) The three approaches MA, FA, and LA are compared to full linear viscoelastic computations (LVE) in a numerical study in Sect. 3.4.
With the dynamic stiffness of the viscoelastic material and the elastic stiffness of the inclusions at hand, we may approximate the dynamic stiffness of the composite material based on elastic computations. We will refer to this approach as surrogate elastic computation (SEC) approach in the remainder of the manuscript. For convenience of the reader, we summarize the necessary computational steps in the following box.

Computational investigations of the dynamic stiffness of SFRPs
In this section, we validate the derivations of the previous sections. Indeed, in Sect. 3.3, we proposed to approximate the dynamic stiffness with an efficient scheme. Instead of computing the material response of a viscoelastic material subjected to cyclic loading, we may compute a single "elastic " stiffness.
To verify the accuracy of the approximation, we proceed in two steps. For the approximation proposed to be valid, we assumed the phase shift of the matrix material to be small. In a first step, we thus test the effect of large phase shifts to quantify what "smallness" means in engineering practice. Here, we use artificial parameters for the matrix material that are designed to show large phase shifts to thoroughly test the approximation quality. Since the phase shift in the material can be designed independently of the number of Maxwell elements, we chose to use a single Maxwell element material in this step. In a second step, we return to realistic material parameters identified for PBT in Sect. 3.2 to validate that typical phase shifts in engineering materials are within the range of applicable phase shifts identified in the first step.
For the first step, we choose a simple microstructure with fibers of aspect ratio ten, shown in Fig. 9a. To compare the effect of the shift approximations MA, FA and LA. A viscoelastic material with large phase shift δ matrix at the applied loading frequency is chosen. We subject the material to stress-driven cyclic loading with a stress amplitude of σ a =15 MPa and a frequency of 3 Hz. To keep the numerical effort at a minimum, only one Maxwell element for the matrix material model is used. Its material parameters are E 0 =1 GPa, E 1 =10 GPa and ν =0.4. We study the material behavior for two different viscosities, η 1 =160 GPa·s leading to a phase shift of δ matrix =57.6 • at the applied load frequency and η 1 =25 GPa·s leading to a phase shift of δ matrix =27.4 • , respectively. We vary the Young's modulus of the fiber material between 10 1 −10 7 MPa, its Poisson ratio is set to ν =0. 22.
For the linear viscoelastic computation (LVE), we evaluate the strain of the homogenized material ε(t) and identify the maximum ε max x x and minimum value ε min of the linear viscoelastic material (LVE). In Fig. 9, the predicted dynamic moduli are compared to the SEC model, using the three different approximation approaches MA, FA and LA. For a large phase shift, shown in Fig. 9b, the FA overestimates the dynamic stiffness of the composite material significantly. Both MA and LA approximate the composite behavior comparatively well. In the region of typical stiffness contrast between matrix and fiber, i.e., for log 10 (E fiber /GPa) = 10 4 − 10 5 , the LA model matches the composite behavior more accurately than MA. For smaller phase shifts, the deviation between the models decreases. Comparing the predicted dynamic moduli in Fig. 9c for a smaller phase shift to the moduli in Fig. 9b, the deviation between MA and LA is barely visible. Both models approximate the linear viscoelastic model with reasonable accuracy.
To summarize, the deviations from the full computation using the LVE model are smaller for the LA than for the FA and the MA approach. Thus, we fix the LA approach to study the material behavior in the remainder of the manuscript. With these insights at hand, we return to the industrial fiberreinforced PBT. The material parameters of the matrix are chosen as described in Sect. 3.2. The material parameters were identified using tensile tests on the composites at different strain rates. The fiber inclusions are modeled as an isotropic linear elastic material with the Young's modulus E = 72 GPa and the Poisson's ratio ν = 0.22. The fiber structure is layered, with the fiber orientations in each layer chosen according to the structure plotted in Fig. 8a and its realization shown in Fig. 15.
The results for excitation of the fiber structure in 0 • -and 90 • -direction are shown in Fig. 10. The hysteresis loops plotted are the result of the linear viscoelastic computation. Here, for each load case the fifth hysteresis loop is shown. For better visualization, each of the hysteresis loops is shifted by some strain value ε 0 . For the 0 • -direction, the strain at the first time step of the plotted loop starts at 0%, 0.1% and 0.2% for the computations at 3 Hz, 5 Hz and 10 Hz, respectively. For the 90 • -direction, the strain at the first time step of the plotted loop starts at 0%, 0.2% and 0.4% for the computations at 3 Hz, 7 Hz and 12 Hz, respectively. The dashed lines represent the Young's moduli computed by the SEC model in the respective direction. In Fig. 10a, different stress amplitudes at a frequency of 10 Hz are shown. The hysteresis loops are shown in blue, green and red. As discussed in Sect. 2.2 for the one-dimensional case, we observe that the dynamic stiffness of the composite material is independent of the load amplitude. In Fig. 10b, the hysteresis curves for excitation in 0 • -direction of the fiber structure are shown for different frequencies. The colored dashed lines, which represent the dynamic stiffness of each hysteresis loop computed by the SEC model, are shifted to higher strain values for better comparison. With increasing frequency, the dynamic stiffness increases. In Fig. 10c, the hysteresis curves and their dynamic stiffness is shown for the 90 • -direction. The SEC model represents the dynamic stiffness of the linear viscoelastic computations quite well. To study the approximation of the viscoelastic model quantitatively, the axial dynamic stiffness obtained from both models is shown in Fig. 11. We observe that the SEC model represents the dynamic properties of the material model in this fifth cycle very well. The deviation between the stiffness computed by the time resolved vis-  To get an idea of the speed up obtained by using the SEC model instead of the viscoelastic computation, we computed one load cycles discretized by 40 time steps for the viscoelastic model on the layered fiber structure. The number of Maxwell elements is 4 and we used the FeelMath-solver [93] on 8 nodes. The resulting runtime is approximately 2.4h. In contrast, a single elastic time step on the same fiber structure takes approximately 2min on the same hardware. To obtain the material stiffness, we need to compute six load cases both in the SEC and the viscoelastic model. It is necessary to compute a few cycles to reach the periodic orbit. Assuming that six cycles are enough to reach the fixed point using the viscoelastic model, we get a speed-up factor of about 432. Additionally, the SEC algorithm highly increases the memory efficiency of the computation. For a viscoelastic material with 4 Maxwell elements, i.e., 4 internal strain fields, we need to store 24 scalar fields. On the other hand, the SEC computation does not need any additional internal variables. Taking a microstructure discretized by 256 3 voxels as an example, storing one scalar field of double precision values takes 0.125 GB in memory. Thus, we can reduce the memory use by 3 GB in comparison with the viscoelastic material model.

Influence of fiber volume content and fiber orientation
To get further insights into the influence of viscoelastic properties of the polymer matrix onto the cyclic behavior of the composite material, we study the dependence of the material behavior on the fiber volume content and the fiber orientation. Goméz et al. [71] studied the influence of the fiber volume content on the dynamic properties of short-fiber reinforced PBT at elevated temperatures by experiments. We use their experimental results to compare trends observed in numerical predictions obtained by the SEC model to experimental observations. Goméz et al. [71] measured the frequency-dependent material properties of PBT at 50.2 • C. In Fig. 12, the dependence of the storage modulus of PBT on the frequency is shown in black. To reproduce this behavior by a linear viscoelastic material model, we chose to use 15 Maxwell elements, about one for each decade in the measured frequency range. To identify the parameters for PBT at 50.2 • C, we proceed with a fit of the material model to the measured loss and storage modulus using an algorithm similar to the one described in Woldekidan [94]. The obtained moduli E j and relaxation time constants τ j = 1/β j are given in Table 3. The Poisson's ratio of all elements is set to ν = 0.4 for all further computations, as commonly identified for PBT. The resulting fit of the material behavior is shown by the red line in Fig. 12. We observe that the material behavior can be represented accurately by the viscoelastic model. The dynamic modulus is mainly governed by the storage modulus as the phase shift in the material behavior is not pronounced (Fig.  13).
Goméz et al. [71] did not report on the fiber orientation inside the experimental specimens. In the present study, we use fiber structures with the eigenvalues of λ 1 = 0.5 and λ 2 = 0.3 in their second-order Advani-Tucker fiber orientation tensor [95]. Higher order momenta are recovered Fig. 13 Loss modulus by the exact closure approximation [96]. We compare fiber structures with 0%, 12% and 29% fiber volume content. For the fiber structures of 12% and 29% fiber volume content, two realizations are generated. The first has a volume edge-length L vol of 1.32 fiber lengths fiber , i.e., we use a volume-fiber length ratio of L vol / fiber = 1.32. In the second structure this ratio is increased to L vol / fiber = 1.77. The number of voxels per fiber diameter is thereby kept constant (6.4 voxels per fiber diameter). The fiber structures with a ratio of L vol / fiber = 1.77 and with 12% and 29% fiber volume content are shown in Fig. 14, visualized with Geo-Dict 1 . The structures were generated using the Sequential Addition and Migration (SAM) algorithm [97].
In Fig. 14, the dependence of the dynamic modulus of the composite material on the frequency is shown. The increase in fiber volume content is reflected in a shift of the dynamic modulus to higher values. The slope of the curve is slightly increased, i.e., higher frequencies show a stronger influence on the composite than on the plain matrix material. These  (Fig. 15).
We further study the influence of the fiber orientation on the viscoelastic material behavior in the periodic orbit. In Fig. 16, the dependence of the composite's dynamic stiffness on the frequency is shown for three fiber structures with different fiber orientations. Each of the fiber structures is subjected to stress-driven extension in x-direction. The reinforcements in the direction of loading lead to a stiffening of the composite material. Thus, the uni-directional fiber structure shows the highest dynamic modulus in x-direction, as all fibers point in loading direction, followed by the planarisotropic and the isotropic fiber structure. Higher frequencies lead to a higher dynamic modulus of the matrix and thus the composite. This effect becomes apparent for all three structures. The increase of the composite dynamic stiffness with increasing frequencies is most pronounced for the isotropic structure, as the matrix behavior is more dominant in this case.

Material model
In polymer-based materials under fatigue loading, we observe a steady degradation of the mechanical behavior over a range of thousands or millions of cycles. This degradation is observable in a decrease of the measured dynamic stiffness. On the other hand, short-term effects of viscoelastic nature can significantly influence the measured stiffness, as discussed for linear viscoelastic materials in the previous Sect. 3. As both of this effects, damage and viscoelasticity, manifest at different time scales, we chose to decouple these effects and model them in two subsequent steps. First, we compute the dynamic stiffness of the viscoelastic material at high cycles, i.e, we evaluate the function C dyn ( f ) using the SEC model discussed in Sect. 3. Secondly, we use the evaluated long-term viscoelastic dynamic stiffness to model the evolution of the dynamic stiffness with a fatigue-damage model. In principle, this approach is independent of the chosen damage model, i.e., the engineer is free to choose a proper damage variable D and to design an appropriate damageevolution law dD/dN = f (ε, D, . . .) to suit the material behavior.
We demonstrate the effects of the proposed ansatz on a recently published compliance-based model [98] in the following. The elastic compliance-based material law is defined based on the free energy density w (ε, D) = 1 2 with the elastic stiffness tensor C 0 , the strain amplitude tensor ε and a scalar damage variable D ∈ [0, ∞), and the dissipation potential The damage evolution, directly formulated in logarithmic cycle space D ≡ d(D)/d log 10 N is governed by the damage speed parameter α. We refer to Magino et al. [98, section 2] for a detailed discussion of the model. We introduce N = log 10 N as short-hand notation for the variable in logarithmic cycle space. To account for the frequency-dependent long-term behavior of the viscoelastic matrix material, we replace the elastic stiffness with the dynamic stiffness C dyn . This yields w (ε, D, f ) = 1 2 16 Influence of fiber orientation on dynamic stiffness for isotropic (iso), planar-isotropic (piso) and uni-directional (ud) fiber structures The modification effects both the stress-strain relation (4.4) where σ is the stress amplitude tensor in the current cycle, and the damage evolution equation i.e., the damage evolution is governed by the load amplitude and the frequency. Equation (4.5) reveals that irreversibility of the damage evolution is ensured if the dynamic stiffness C dyn is positive semidefinite and the inequality D(N 0 ) ≥ −1 holds for some N 0 . Note that the damage evolution is directly formulated in cycle space in the work at hand. For a discussion of the evolution in time space and damage evolution under loading and unloading, we refer to Magino et al. [75]. Experimental studies revealed that damage evolution inside the polymer matrix and the resulting stiffness degradation of the composite material is indeed influenced by the load amplitude [99][100][101]. Moreover, a recent study of Imaddahen et al. [102] reports that the fatigue process in short-glass fiber reinforced polypropylene depends on the strain rate or frequency. They observed an improvement of the fatigue strength at higher frequencies, i.e., a higher number of bearable load cycles at higher frequencies. As a possible explanation for this result they refer to the static experiments of Fitoussi et al. [103] where the authors showed that damage effects at high strain rates are delayed in polymer composites. This observation suggests that the damage evolution in thermoplastic materials is reduced by an increase in frequency.
To uncover the influence of an increased frequency in the proposed model (4.2)-(4.3), we study the influence of the frequency on the stiffness degradation in simulations. Prior to this study, we investigate the necessary spatial resolution to accurately predict the fatigue damage evolution in the composite material. As the reference, we chose a resolution with 6.4 voxels per fiber diameter, as shown in Fig. 17b for an isotropic fiber structure. We compare this resolution to a finer resolution with 12.8 voxels per fiber diameter (h/2) and a coarser resolution with 3.2 voxels per fiber diameter (2 × h) in Fig. 17a. The effective stiffness of the reference and the finer resolutions are close, i.e., the relative deviation at N = 0 is 0.5 %. The deviation becomes larger with increasing mesh width, thus we fix the resolution of the reference structure for further studies. Subsequently, we turn our attention to the stiffness degradation of the composite at different frequencies. In Fig. 18, the degradation of the dynamic Young's modulus in x-direction E x dyn is shown for the structure with isotropic fiber orientation shown in Fig. 17b and stress-driven uniaxial extension of the material with a load amplitude of 50 MPa. The material parameters are listed in Table 2. The stiffness of the material increases with increasing frequency. Thus, the evolution of the dynamic modulus is shifted to higher values for higher frequencies. Plotting the relative evolution of the Young's modulus, we observe that higher frequencies lead to slower degradation. Higher frequencies lead to a slower degradation. This in accordance with the experimental observations of Imaddahen et al. [102].

Experiments versus computational predictions
To validate the viscoelastic damage model discussed in Sect. 4.1, we performed fatigue experiments at different load amplitudes and frequencies. The experiments were done on specimen geometries as shown in Fig. 19a. The geometries were cut from an injection-molded plate. The layered fiber orientation structure inside the plate, which evolves due to the injection process, was characterized by µCT-scans and is shown in Fig. 19b. For the scans, a small volume at the center of the injection molded plate was cut out. Subsequently, the fiber orientation in nine layers subdivided over the thickness of the plate was analyzed, see Hessman et al. [89] for more details on the method. We investigate two types of specimens, 0 • -oriented and 90 • -oriented. By 0 • -oriented, we refer to specimens that were cut from the plate in the direction of the mold injection and in which the resulting fiber orientation is thus mainly oriented in the loading direction. By 90 • -oriented specimens, we mean specimens that were cut transverse to the mold injection direction. The specimens were subjected to cyclic loading with stress ratios R = 0. The different load cases are listed in Table 1. To get an estimate for the scattering of the measurements, each load case was conducted three times. For the computational predictions, we use the (visco)elastic material properties of the matrix and fiber listed in Table 2. For the damage evolution in the matrix material, a straightforward parameter-identification strategy [69] was followed to identify the parameter α. As the damage evolution equation depends linearly on the parameter α, it may be rearranged in the form [69]   Thus, using a prediction of the stiffness computed with any parameter value α can be easily utilized to identify the parameter α by comparison to experimental data. In this manner, we computed the stiffness degradation of the 0 • -oriented spec-imens on the test specimen and identified α by rescaling. This led to the identified parameter of α = 0.18 MPa −1 . Subsequently, the reduced order model was trained for the frequencies of interest. The reduced order model is implemented in the finite element solver Abaqus [104] via a user material subroutine (UMAT). In a primary step, the composite material is characterized computationally via FFT-based precomputations, from which the system matrices of the reduced order model are identified and stored for every of the precomputed microstructures, see Fig. 20. For more details on the database generation process, i.e., the offline phase, we refer to Magino et al. [75, section 5.2]. In the online phase every element of the macroscopic specimen is assigned with a fiber orientation, encoded by the second order Advani Tucker fiber orientation tensor. The sum of its eigenvalues is 1 and the space of possible eigenvalues λ 1 , λ 2 can be represented by the fiber orientation triangle [69], shown in the database box of Fig. 20. In the online computation, each element, according to its assigned orientation, is associated to the system matrices of the three closest precomputed orientations via the fiber orientation interpolation approach [69]. This is done once at the beginning of the computation using the UEX-TERNALDB routine of Abaqus. During the computation, the material law implemented in an user subroutine (UMAT) is evaluated. For the macroscopic specimen, we considered a layered fiber structure resolving nine layers, see Fig. 19b. In Fig. 21, we compare the resulting stiffness evolution in the component simulations to the experimental measurements. Both the measured and the predicted dynamic stiffness increase with increasing frequency. An increase of the frequency by several Hz results in an increase of the dynamic stiffness by several hundred MPa. The computational predictions of this increase match the experimental values quite well, both qualitatively and quantitatively. The scattering of the measurements of 0 • -oriented specimens is significantly higher than for the 90 • -specimens. For 0 • , the initial predicted dynamic stiffness matches the measured mean value and the predicted degradation is within the experimental scattering up to the last decade of the experiment. For the 90 •oriented specimen the dynamic stiffness within the first few hundred is overestimated by the computational predictions. With increasing stress amplitude, the deviation increases.
The model at hand accounts for linear viscoelastic effects, which is typically valid for small strain rates in polymer materials. However, the local strain rates in the matrix material increase at higher stress amplitudes. At large strain rates, non-linear viscoelastic effects may have a significant influence. This fact might explain the increasing deviation of the linear viscoelastic prediction from experimental measurements. Additionally, the current model does not account for possible (visco)plastic effects. Yet, for cycle numbers exceeding 100, the deviation between computational predictions and experimental measurements is acceptable. The applied damage model does not account for localization effects by construction and is therefore unable to represent a negative second derivative in the stiffness degradation. Rather, we focus on the reproduction of the approximately linear trends in the stable stiffness degradation [75,98]. The damage processes leading to the stiffness degradation within the last decade of the fatigue process are most likely already localized. This is responsible for the higher deviations between computational and experimental results with increasing cycle number. To improve predictions in this regime, a localizing damage model needs to be applied.

Conclusion
In the work at hand, the dynamic stiffness of short-fiber reinforced composite materials under high-cycle fatigue loading at different stress amplitudes and frequencies was studied. Experiments reveal a frequency dependence that is not typically accounted for in contemporary existing elastic fatigue damage models. To accurately predict the dynamic stiffness in fatigue experiments, this dependence needs to be accounted for.
In the contribution, we restrict ourselves to linear viscoelastic material model for the thermoplastic matrix. At high stress levels locally occurring in the composite material, this is an approximation. Yet, we are able to reproduce the composite material behavior at a range of strain-rates already with this simple material model, see Sect. 3.2. For cyclic loading of linear viscoelastic materials, a frequency dependent relation between the stress and strain amplitudes at the steady-state is known. We use this relation to approximate the dynamic stiffness of a short-fiber reinforced polymer with equal Poisson's ratios in every Maxwell element valid for small phase shifts. The effective dynamic stiffness may be determined as an elastic effective stiffness (SEC) and is shown to be rather accurate for short-fiber reinforced PBT, see Sect. 3.4. The proposed algorithm is both memory and time efficient.
We validated the approach using experiments of shortfiber reinforced PBT at different load amplitudes, frequencies and fiber orientations. The model is able to predict the frequency dependence in the dynamic stiffness degradation. Additionally, the trend of a slower damage evolution at higher frequencies is reproduced, as the the damage evolution equation is governed by the dynamic stiffness.
To calibrate the model, we need to identify the viscoelastic material parameters of the matrix (see Sect. 3.2) from strainrate experiments on the composite and the damage parameter α from a single fatigue experiment for the matrix. The experimental expense is thus significantly reduced in comparison to a multitude of fatigue measurements. It remains to propose an interpolation approach for the frequency dependent material behavior in the reduced order model. For the computations on the microstructure, it is straightforward to account for a change in frequency by adjusting C dyn ( f ) to the new frequency. However, for the reduced order model, a database for every frequency of interest has to be trained and the macroscopic computation is run at a constant prescribed frequency. To account for a frequency change, the reduced order model thus needs further development.
Additional effects that are of interest in short-fiber reinforced components include temperature dependence. Since polymers are known to show interchangeable properties in frequency and temperature change, it would be interesting to investigate an extension to temperature dependency. Having understood temperature and frequency dependence, the experimental investigation of a wider range of frequencies would be possible. A more thorough experimental investigation of a wider range of frequencies and their effect on the stiffness degradation is a vital step prior to the an accurate prediction of dynamic stiffness.