# Continuum approach for modeling fatigue in amorphous glassy polymers. Applications to the investigation of damageratcheting interaction in polycarbonate

Sami Holopainen, Thierry Barriere, Gang Cheng, Reijo Kouhia

► To cite this version:

Sami Holopainen, Thierry Barriere, Gang Cheng, Reijo Kouhia. Continuum approach for modeling fatigue in amorphous glassy polymers. Applications to the investigation of damageratcheting interaction in polycarbonate. International Journal of Plasticity, 2017, 91, pp.109-133. 10.1016/j.ijplas.2016.12.001. hal-02131311

HAL Id: hal-02131311

https://hal.science/hal-02131311

Submitted on 3 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Continuum approach for modeling fatigue in amorphous
glassy polymers. Applications to the investigation of
damage-ratcheting interaction in polycarbonate.

Sami Holopainen†*, Thierry Barriereᵇ, Gang Chengᵇ, Reijo Kouhia†

† Tampere University of Technology
Department of Mechanical Engineering and Industrial Systems
P.O. Box 589, FI-33101 Tampere, Finland

ᵇFEMTO-ST Institute, Applied Mechanics Department
24 Rue de l'Epitaphe, 25000 Besançon, France

ᵇMechanical Engineering Design Department, INSA Lyon, 85 Rue Henri Becquerel, 01100 Bellignat, France

**Abstract**

In this article, we propose an approach suitable for modeling isothermal fa-
tigue in amorphous polymers. The theory is formulated in a rate form within
continuum mechanics framework without the need to measure damage changes
per loading cycles. Using the approach, contribution of ratcheting to fatigue
of polycarbonate (PC) was investigated and the results were compared to pre-
vious experimental observations. When subjected to uniaxial stress-controlled
cyclic loadings, ratcheting deformation apparently occurs and increases with
mean stress and amplitude. The development of ratcheting deformation shows
an initial growth followed by a decrease to almost a constant growth rate which
occupies majority of the total lifetime. Ratcheting behavior under multiaxial
stress states was also investigated based on finite element analyses of a dogbone-
shaped test specimen. The results show that fatigue damage develops at the sites
following closely the localized plastic deformation and increases with ratcheting
deformation during cyclic loadings. The results indicate that the ratcheting
behavior can be exploited in the evaluation of the entire fatigue lifetime.

***Keywords:*** fatigue damage, amorphous polymers, ratcheting, endurance surface, FEM

**1. Introduction**

Manufacturers of polymer materials are interested in the assurance of com-
ponents lifetime, especially when the lifetime cannot be easily inspected or may

<footnote>*Corresponding author, e-mail: sami.holopainen@tut.fi</footnote>
lead to a catastrophe in service, Maxwell et al. (2005). The total annual loss of engineering components due to fatigue failure has been assessed to be the most dominant causing significant financial losses. Thus, design against fatigue constitutes an integral part of mechanical engineering analysis. Examples of mechanical components that are manufactured from polymers and may experience fatigue during their service life are from fields such as sporting goods, automotive glazing, marine structures, aeronautics, and armour. The design of such constructions could benefit from capable models and the strong computational capability currently available. Moreover, due to the costly and time-consuming testing of long-term fatigue life, which still suffers from its sensitivity to environmental conditions and previous operation histories, there is a need to develop approaches that can simulate tests well.

Knowledge of the ultimate behavior of amorphous glassy polymers has developed rapidly while only a little of this research has been devoted to fatigue behavior, its modeling, and testing, cf. Krairi and Doghri (2014, p. 176). Specific features of these matters are seen e.g. from Lesser (2002); Janssen et al. (2008b); Drozdov (2011); Del Vecchio et al. (2014); Lugo et al. (2014); Xi et al. (2015); Mortazavian and Fatemi (2015). To explore microstructural characteristics of polymers with regard to their fracture toughness, the damage mechanisms ahead of the crack tip have been paid considerable attention, see Li and Chandra (2003); Pijnenburg et al. (2005); Pasta (2011) including references therein. Much research has also devoted to the investigation of fatigue crack propagation in polymers and their compounds, see Fang et al. (2008); Nittur et al. (2013); Ravi Chandran (2016); Kanters et al. (2016) to mention a few. However, these fracture mechanics approaches do not consider the crack initiation stage which may cover over 95% of the total fatigue life of amorphous polymers, Marissen et al. (2001); Janssen et al. (2008b); Lugo et al. (2014).

Fatigue failure of amorphous polymers in their glassy state (termed amorphous glassy polymers) is generally attributed to a two-step process. In the first, initiation step, failure is typically attributed to deficiencies or impurities producing significant stress concentrations which exceed the strength limits of the material, Marissen et al. (2001); Lu and Kim (2007); Drozdov (2011). During repeated loadings, those defects can nucleate and grow during the service life at stress levels well below the nominal yield strength, Kanters et al. (2016). Under such conditions, fatigue is influenced by the localized yield-like deformation process which provides fatigue crack initiation sites and thus controls most of the fatigue life (cycles N to failure) and is of specific interest in the applied fatigue stress S (S-N curve), Lemaitre and Desmorat (2005); Ravi Chandran (2016).

The second, propagation step is characterized by the damage growth through the coalescence of micro-cracks and propagation of small cracks to form large cracks which ultimately lead to component failure, Lesser (2002); Lemaitre and Desmorat (2005); Drozdov (2011); Nittur et al. (2013); Lawrimore et al. (2016). However, the duration of the initiation step is invariably orders of magnitude greater than the propagation period and thus plays the most important role in the fatigue behavior, Janssen et al. (2008b); Lugo et al. (2014). Based on this observation, the influence of crack propagation in the material behavior is often omitted in
the fatigue models.

Although the literature concerning the fatigue damage in amorphous poly-
mers is still scarce, an extensive research has been devoted to investigation
of damage mechanisms of several polymers under monotonic loadings and re-
peated loading loops producing large inelastic deformations. Zaïri et al. (2011);
Holopainen (2014); Seeling and Van der Giessen (2015); Krairi and Doghri (2014);
Lawrimore et al. (2016) represent the recent research literature of this field. This
research also has relevance for fatigue since the failure kinetics of such polymers
shows many similarities under cyclic and static (creep) loadings, Janssen et al.
(2008b); Lugo et al. (2014). However, material subjected to cyclic loadings may
fail at stress levels substantially lower than observed under monotonic loadings,
Abdel-Karim (2005); Drozdov (2011). When the stress is below the (static) yield
strength, the material response is primarily (visco)elastic but may ultimately
fail after a number of load cycles due to the accumulation of plastic strains,
Kang (2008). This behavior is known as ratchetting (cyclic creep) which closely
contributes to the fatigue of polymers under cyclic loadings, i.e. rapidly grow-
ing ratcheting strain (cyclic creep strain) reduces fatigue life, Liu et al. (2008);
Zhang et al. (2015). Xi et al. (2015) conducted experiments on polycarbonate
(PC) polymers and showed that the polymers with large molecular weights are
more prone to ratcheting than those with lower weights. Loading conditions,
such as mean stress, amplitude, and temperature, influence also ratcheting de-
velopment, Kang (2008); Zhang et al. (2015); Lu et al. (2016). Although a lot
of approaches are proposed for the simulation of ratcheting behavior, they are
rarely aimed at amorphous polymer materials or require modeling of hysteresis
loop action leading to a large number of material parameters to be defined,
Liu et al. (2008). Based on a ratcheting model for metals Liu et al. (2008) pro-
posed a simple approach that captures exactly the steady ratcheting strain be-
havior of polymethyl methacrylate (PMMA) under different temperatures and
loading conditions. However, this model is only indended to the description
of ratcheting accumulative rates at steady stage, while the initial ratcheting
behavior and the behavior prior to ultimate failure are neglected. Based also
on metals, Liu et al. (2010) proposed a stress-based, multiaxial fatigue failure
model suitable for the investigation of fatigue-ratcheting interaction. Although
the model provides good predictions, it suffers from its restriction to a reduced
region of both mean stress and shear mean stress.

It is largely acknowledged that ratcheting deformation is triggered by accu-
mulated plastic stretching present in cyclic processes, Liu et al. (2008). How-
ever, the accumulation manner of ratcheting strains or ratcheting vs viscous
(elastic or plastic) deformation are still open questions, at least in amorphous
polymers. Also, through existing literature in the field, Li et al. (1995); Lesser
(2002); Kang (2008); Lugo et al. (2014); Mortazavian and Fatemi (2015); Ravi Chandran
(2016), key morphological or microstructural mechanisms that could explain the
origin of the plastic deformation and subsequent progress of fatigue damage un-
der cyclic loadings are presently unclear.
Research design

This article continues by introducing a model suitable for predicting fatigue in amorphous glassy polymers and by discussing the phenomenon through theory. The governing constitutive model employed is an extension of the celebrated Boyce-Park-Argon (BPA) model for predicting inhomogeneous plastic deformation in amorphous polymers, Boyce et al. (1988). Thereafter, Section 2.4 describes the proposed, multiaxial fatigue model and its numerical treatment in detail. Modeling of fatigue behavior per se is based on an appealing model introduced in Ottosen et al. (2008), which model is formulated in continuum mechanics framework by using evolution equations that make the definition of damage changes per cycle redundant, i.e. cycle-counting techniques do not need to be applied. The approach is calibrated to data taken from both cold drawing experiments and isothermal fatigue tests on dogbone-shaped PC-specimens. Section 3 is dedicated to the evaluation of the proposed approach. The studies involve tangible examples for ratcheting behavior, and the model results are compared with previous research results. The article closes with the conclusions as well as with the limitations and future research avenues.

## 2. The model

### 2.1. Kinematics - basics

The Kröner-Lee decomposition,

$$F = F^e F^p, \qquad (1)$$

of the deformation gradient $F$ is applied, i.e. $F^p$ and $F^e$ define the plastic and elastic contribution, respectively, cf. the original paper Boyce et al. (1988). The symmetric, positive definite stretch tensor $v$ is given by

$$F = vR, \qquad (2)$$

and $R$ is the rotation tensor. Motion of a material body is defined through the velocity field $v$ that defines the spatial velocity gradient $l$, i.e.

$$l := \mathrm{grad}(v) = \dot{F}F^{-1}. \qquad (3)$$

In particular, the accent dot means the material time derivative and the quantities in the intermediate placement will be highlighted by the bar. The split (1) in (3) further provides expressions for the elastic and plastic contributions of $l$, Boyce et al. (1988); Holopainen (2013).

For later use, the elastic and plastic deformation are given in terms of symmetric elastic and plastic deformation tensors defined as

$$\bar{C}^e := F^{e,T} F^e, \quad \bar{C}^p := F^p F^{p,T} \qquad (4)$$

where the tensor transpose is marked by the letter T.
<figure><img src="image_2.png" /><figcaption>Figure 1: Deformation of a solid body and the mappings between different configurations, $\tilde{N}$, $\tilde{N}^p$, $N^R$, and $N^e$. In accordance with the BPA model, Boyce et al. (1988), $R^e$ is chosen to be unity when the elastic intermediate configuration $N^e$ coincides with $\tilde{N}$.</figcaption></figure>

The elastic deformation in the spatial configuration is defined by the elastic Finger tensor,

$$b^e := F^e F^{e,T}. \qquad (5)$$

The polar decomposition of $F^e$ allows definition of the alignment of the elastic relaxed placement by the symmetric elastic stretch $v^e := \sqrt{b^e}$ and the rotation $R^e$, i.e.

$$F^e = v^e R^e. \qquad (6)$$

Similar to (6),

$$F^p = \bar{V}^p R^p \qquad (7)$$

where $R^p$ is the plastic rotation and $\bar{V}^p$ governs the plastic stretching defined in the relaxed intermediate placement, see Fig. 1.

## 2.2. Constitutive model

Polymer materials show the time-dependent behavior under both monotonic and cyclic loadings, Lesser (2002); Anand and Ames (2006). This viscous effect is due to the macromolecular characteristic of polymers, i.e. the polymer chains need a relaxation time to attain their equilibrium state after deformation. A consequence is that the fatigue behavior is different from elastic metals, i.e. a model's ability to follow the viscous (viscoelastic-plastic) behavior of polymers is mandatory.

For convenience, the governing constitutive model for amorphous polymers is briefly reviewed, see Holopainen (2013). The model is an extension of the celebrated non-affine 8-chain model proposed in Boyce et al. (1988) termed here
<figure><img src="image_3.png" /><figcaption>Figure 2: Idealized structure of chain network according to the 8-chain model in a) its initial and b) deformed state. The base vectors $\bar{N}_\alpha$ are the unit eigenvectors of $\bar{V}^p$. The dimension of the undeformed element is $a_0$ and $\bar{\lambda}^p$ denotes the plastic chain stretch which emerges in the direction $\bar{m}$. The unit vector $\bar{m}$ is defined by the angles $\bar{\theta}$ and $\bar{\phi}$, whereas the plastic network stretch $\bar{\lambda}_{ec}^p$ is related to the direction $\bar{N} := 1/\sqrt{3}\sum_\alpha \bar{N}_\alpha$.</figcaption></figure>

the BPA model. The idealized chain network according to this model is depicted in Fig. 2.

The extended model (EBPA model) captures creep and recovery through the Kelvin-Voigt element that involves an elastic spring 2) and a viscous dashpot 3) as shown in Fig. 3. Its configuration with the elastic spring 1) predicts the stress relaxation. The second Kelvin-Voigt-like element comprises a viscoplastic dashpot 4) arranged parallel with a nonlinear spring 5) that results in the anisotropic material behavior at large strains.

The split of $F^e$ into a viscous and an elastic component is applied, i.e.

$$F^e = F_1^e F_2^e \qquad (8)$$

where $F_1^e$ and $F_2^e$ define the elastic stretching in the spring 1) and 2), respectively, see Fig. 3 and Holopainen (2013) for a more detailed account. The polar decomposition of $F_1^e$ and $F_2^e$ allows definition of the orientation of the elastic intermediate placement as

$$F_1^e = v_1^e R_1^e \quad \text{and} \quad F_2^e = v_2^e R_2^e \qquad (9)$$

where $R_1^e$ and $R_2^e$ are the elastic rotations and $v_1^e$ and $v_2^e$ are the elastic stretch tensors, defined in the current deformed and the elastic intermediate placement, respectively. Furthermore, the elastic deformation is given by the two deformation tensors,

$$b_1^e := F_1^e F_1^{e,T} = (v_1^e)^2 \quad \text{and} \quad b_2^e := F_2^e F_2^{e,T} = (v_2^e)^2. \qquad (10)$$

The product decompositions (1) and (8) in (3) yields

$$\mathrm{sym}(l) := d = d^e + d^p = d_1^e + d_2^e + d^p \qquad (11)$$

where

$$d_1^e := \mathrm{sym}(\dot{F}_1^e F_1^{e-1}), \quad d_2^e := \mathrm{sym}(F_1^{e-T} \bar{l}_2^e F_1^{e-1}), \quad d^p := \mathrm{sym}(F^{e-T} \bar{L}^p F^{e-1})$$

and

$$\bar{l}_2^e := \dot{F}_2^e F_2^{e-1}, \quad \bar{L}^p := \dot{F}^p F^{p-1}. \qquad (12)$$
## Thermodynamics

Based on the so-called dissipation inequality thermodynamics of the model is investigated. Two potential functions termed the Helmholtz free energy $\varphi$ and the dissipation potential $\phi$ are needed for the treatment, Lemaitre and Chaboche (1994). Assuming that isothermal conditions prevail, the specific Helmholtz' free energy is given by the following elastic, plastic, and damage contributions

$$ \varphi = \hat{\varphi}(b_1^e, b_2^e, \bar{C}^p, D) = \varphi^{e,1}(b_1^e, D) + \varphi^{e,2}(b_2^e, D) + \varphi^p(\bar{C}^p) \quad (13) $$

where $\varphi^{e,1}$, $\varphi^{e,2}$, and $\varphi^p$ account for the evolving stress in the springs 1), 2), and 5) shown in Fig. 3. The free energy also considers the degradation of the material through a damage variable $D$ evolving from an initial value (typically zero) to one.

In terms of the internal power $\mathfrak{W}^{\text{int}}$ and the rate $\dot{\varphi}$, the power of local dissipation in the current placement is governed by

$$ \mathcal{D} = \mathfrak{W}^{\text{int}} - \dot{\varphi} := \tau : d_1^e + \tau : d_2^e + \tau : d^p - \dot{\varphi} \geq 0 \quad (14) $$

where the operator $:$ is defined by the tensor trace as $A : B := \text{trace}(AB^T)$. The dissipation power can also be defined by a complementary dissipation potential $\phi$ given in terms of the viscoelastic, viscoplastic, and the degradation parts as

$$ \phi(\tau_2, \tilde{\tau}, Y) := \phi^{\text{ve}}(\tau_2) + \phi^{\text{vp}}(\tilde{\tau}) + \phi^{\text{d}}(Y) $$

where $\tau_2$ and $\tilde{\tau}$ are driving stresses in the dashpot 3) and 4), see Fig. 3, and $Y$ stands for a driving force for damage. Using the proposed dissipation potential, the dissipation power is defined as

$$ \mathcal{D} = \frac{\partial \phi^{\text{ve}}(\tau_2)}{\partial \tau_2} : \tau_2 + \frac{\partial \phi^{\text{vp}}(\tilde{\tau})}{\partial \tilde{\tau}} : \tilde{\tau} + \frac{\partial \phi^{\text{d}}(Y)}{\partial Y} Y. \quad (15) $$

The rate of the free energy $\dot{\varphi}$ in (14) becomes

$$ \dot{\varphi} = \frac{\partial \varphi^{e,1}}{\partial b_1^e} : \dot{b}_1^e + \frac{\partial \varphi^{e,2}}{\partial b_2^e} : \dot{b}_2^e + \frac{\partial \varphi^p}{\partial \bar{C}^p} : \dot{\bar{C}}^p + \frac{\partial \varphi^d}{\partial D} : \dot{D}. \quad (16) $$

Taking advantage of (10), and using symmetry of $\partial\varphi^e/\partial b_1^e$ yield the following expression

$$ \frac{\partial \varphi^{e,1}}{\partial b_1^e} : \dot{b}_1^e = 2 \frac{\partial \varphi^{e,1}}{\partial b_1^e} b_1^e : d_1^e. \quad (17) $$

In accordance with (17),

$$ \frac{\partial \varphi^{e,2}}{\partial b_2^e} : \dot{b}_2^e = 2 \frac{\partial \varphi^{e,2}}{\partial b_2^e} b_2^e : d_2^e \quad (18) $$

and

$$ \frac{\partial \varphi^p}{\partial \bar{C}^p} : \dot{\bar{C}}^p = 2 \text{sym} \left( \frac{\partial \varphi^p}{\partial \bar{C}^p} \bar{C}^p \right) : \bar{D}^p =: \bar{B} : \bar{D}^p = \beta : d^p \quad (19) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Demonstration of the model: elastic springs 1) and 2), a viscoelastic dash-pot 3), a viscoplastic dashpot 4), and a Langevin spring 5). Difference between the Kirchhoff stress $\tau$ and the backstress $\beta$ defines the driving stress $\tilde{\tau}$.</figcaption></figure>

where $\bar{B}$ is a backstress representing the plastic work reserved in the material, and $\beta$ is its counterpart given in the current placement, i.e.

$$ \beta = F^e \bar{B} F^{e,T}. \qquad (20) $$

Substituting (17), (18), and (19) into (14) and taking (15) into account, yield

$$ \left( \tau - 2 \frac{\partial \varphi^{e,1}}{\partial b_1^e} b_1^e \right) : d_1^e + \left( \tau - 2 \frac{\partial \varphi^{e,2}}{\partial b_2^e} b_2^e \right) : d_2^e - \frac{\partial \phi^{ve}}{\partial \tau_2} : \tau_2 + (\tau - \beta) : d^p - \frac{\partial \phi^{vp}(\tilde{\tau})}{\partial \tilde{\tau}} : \tilde{\tau} - \left( \frac{\partial \varphi}{\partial D} \dot{D} + \frac{\partial \phi^d}{\partial Y} Y \right) = 0. \qquad (21) $$

Since equation (21) has to be fulfilled with all possible thermodynamically admissible processes, the following constitutive equations

$$ \tau = 2 \frac{\partial \varphi^{e,1}}{\partial b_1^e} b_1^e, \qquad (22) $$

and

$$ \dot{D} = \frac{\partial \phi^d}{\partial Y} \qquad (23) $$

with the additional equations,

$$ \tau : d_2^e = 2 \frac{\partial \varphi^{e,2}}{\partial b_2^e} b_2^e : d_2^e + \frac{\partial \phi^{ve}}{\partial \tau_2} : \tau_2, \qquad (24) $$

and

$$ (\tau - \beta) : d^p = \frac{\partial \phi^{vp}(\tilde{\tau})}{\partial \tilde{\tau}} : \tilde{\tau} =: -\dot{s}, \qquad (25) $$

are obtained. In (23), the formula $Y = -\partial\varphi/\partial D$ for the thermodynamic force $Y$ was defined.
Taking advantage of (22) and (24), the power of dissipation takes the form

$$ \mathcal{D} = \frac{\partial \phi^{\text{ve}}}{\partial \tau_2} : \tau_2 + \tilde{\tau} : d^p + Y \dot{D} \geq 0. \quad (26) $$

It will be shown later that the condition (26) is fulfilled by the model.

### Specific potential functions

The chains in polymer are assumed to be randomly oriented when isotropic strain energy functions

$$ \varphi^{e,1} = (1-D)(\frac{1}{2}\kappa_1(I_1^{e,1})^2 + 2\mu_1 J_2^{e,1}), \quad \varphi^{e,2} = (1-D)(\frac{1}{2}\kappa_2(I_1^{e,2})^2 + 2\mu_2 J_2^{e,2}), \quad (27) $$

where $\kappa_1$, $\kappa_2$, $\mu_1$, and $\mu_2$ are the bulk and the shear moduli, are used. The integrity $1 - D$ describes the degradation of the material. This reduction is in agreement with experimental observations which show a notable stress decay prior to complete failure, cf. e.g. Lemaitre and Desmorat (2005). The logarithmic invariants present in (27) are defined as

$$ I_1^{e,1} := \text{trace}(\ln v_1^e) = \ln J^{e,1}, \quad I_1^{e,2} := \text{trace}(\ln v_2^e) = \ln J^{e,2}, $$

$$ J_2^{e,1} := \frac{1}{2}(\ln v_1^e)^{\text{dev}} : (\ln v_1^e)^{\text{dev}}, \quad J_2^{e,2} := \frac{1}{2}(\ln v_2^e)^{\text{dev}} : (\ln v_2^e)^{\text{dev}} $$

where $J^{e,1} := \det(v_1^e)$ and $J^{e,2} := \det(v_2^e)$. The notation, dev, refers to the deviatoric component, defined as $(\cdot)^{\text{dev}} := (\cdot) - \frac{1}{3}\text{trace}((\cdot))\mathbf{i}$ for all second order tensors $(\cdot)$ where $\mathbf{i}$ is the identity.

Making use of functions (27) reveals that $-Y = \partial\varphi/\partial D$ is always positive in (26), i.e. the free energy needed for a damage increase constantly reduces.

The first part of the dissipation potential is defined to be as

$$ \phi^{\text{ve}} := 1/2\tau_2 : \eta^{-1} : \tau_2 \quad (28) $$

where $\eta$ is a constant, symmetric isotropic fourth order tensor defined subsequently. Defining also that $\tau_2 := \eta : d_2^e$, it follows from (22), (24), (27), and (28) that the elastic constitutive law is based on the equilibrium

$$ \tau = (1 - D)\mathcal{L}_1^e : \ln v_1^e = \eta : d_2^e + (1 - D)\mathcal{L}_2^e : \ln v_2^e, \quad (29) $$

where $\mathcal{L}_1^e$ and $\mathcal{L}_2^e$ stand for the standard fourth order elastic stiffness tensor, i.e. generally $\mathcal{L}^e := 2\mu(\mathfrak{J} + (3\kappa - 2\mu)/(6\mu)\mathbf{i} \otimes \mathbf{i})$ in which the tensor $\mathfrak{J}$ is the fourth order identity. In subsequent treatments, the bulk modulus $\kappa$ and the shear modulus $\mu$ will be replaced by the Young's modulus $E$ and the Poisson's ratio $\nu$.

The viscous damper in (29) is given by

$$ \eta = \eta_1(\mathfrak{J} - 1/3\mathbf{i} \otimes \mathbf{i}) + \eta_2\mathbf{i} \otimes \mathbf{i} \quad (30) $$

where $\eta_1$ and $\eta_2$ are the viscosities that govern elastic shear and volumetric deformation, respectively, Reese and Govinjee (1998). The notation $\otimes$ here means
the standard tensor product, i.e. the components of $\eta$ in an orthonormal coordinate system take the form

$$ \eta_{ijkl} = \frac{\eta_1}{2} (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk} - 2/3\delta_{ij}\delta_{kl}) + \eta_2\delta_{ij}\delta_{kl}. \quad (31) $$

For the sake of simplicity, a usual assumption that the latter term in (30) and (31) vanishes is done, i.e. $\eta_2 = 0$, cf. Boyce et al. (1988); Reese and Govinjee (1998). Then, only a few new material parameters, $E_2$, $\nu_2$, and $\eta_1 := \eta$, enter the proposed model.

The rate of the plastic deformation $d^p$ is taken to align with the spatial normalized direction of $\tilde{\tau}^{\text{dev}}$, i.e.

$$ d^p := \dot{\gamma}^p n, \quad n := \frac{\tilde{\tau}^{\text{dev}}}{\sqrt{2\tau}}, \quad \tilde{\tau}^{\text{dev}} := \tau^{\text{dev}} - \beta^{\text{dev}}, \quad \tau := \sqrt{\frac{1}{2}\tilde{\tau}^{\text{dev}} : \tilde{\tau}^{\text{dev}}} \quad (32) $$

in which the use of the deviatoric components reveals that $J^p := \det(F^p) = 1$, i.e. the plastic deformation is incompressible. It then follows that $J = J^e := \det(F^e)$. The development of the plastic deformation is governed by a positive plastic shear strain rate $\dot{\gamma}^p$, see Boyce et al. (1988); Holopainen (2013) for a more detailed account. From (32) follows then that $\tilde{\tau} : d^p = \sqrt{2\tau}\dot{\gamma}^p \ge 0$.

Function of the viscoplastic dashpot is defined by the following evolution equation

$$ \dot{s}_1 = h_1(1 - s_1/s_{ss})\dot{\gamma}^p, \quad s_1(0) = s_0, \quad (33) $$

where $h_1$ and $s_{ss}$ are positive material parameters. In (33), $s_1 > s_{ss}$ when $\dot{s}_1 < 0$. Although $-\dot{s}_1 \neq -\dot{s} = \sqrt{2\tau}\dot{\gamma}^p \ge 0$ in (25), $\dot{s}_1$ also fulfills the condition (26) for the dissipation inequality.

In terms of a non-affine plastic network stretch,

$$ \lambda_{\text{ec}}^{\text{p}} := \sqrt{\text{trace}(\bar{C}^{\text{p}})/3}, \quad (34) $$

the orientational hardening is governed by the backstress

$$ \bar{B} = \frac{C^R}{3\lambda_{\text{ec}}^p} \sqrt{N} \mathcal{L}^{-1} \left( \frac{\lambda_{\text{ec}}^p}{\sqrt{N}} \right) \bar{C}^p \quad (35) $$

where $C_R$ and $N$ are material parameters. When evaluating the inverse of the Langevin function $\mathcal{L}$ present in (35), the Pade approximation is employed.

By substitution of (28) and (32) to (26), the power of local dissipation becomes

$$ \mathcal{D} = d_2^e : \eta : d_2^e + \sqrt{2\tau}\dot{\gamma}^p + Y\dot{D} \ge 0, \quad (36) $$

i.e. the dissipation inequality is always fulfilled since $Y \ge 0$ and $\dot{D} \ge 0$.

## 2.3. Numerical treatment of the constitutive model

Since long-term periods are investigated, an Euler backward method, which enables large time steps, is applied. This method is reported to be efficient for
the modeling of rate-dependent solids, Weber and Anand (1990). To compress the notation, the index n+1 for the updated state is neglected, and the quantities solely at the known state t_n are symbolized by n.

In analogy with the BPA model, the elastic deformation is presumed to be irrotational and therefore the elastic rotation $R^e$ is identity. Moreover, $\tilde{W}^P$ is nonzero and is calculated by introducing an algorithmic plastic spin $\tilde{W}^P$, Wallin and Holopainen (2012). Due to the imposed symmetry of $F^e$, $\tilde{W}^P$ is skew-symmetric ('skew'), and the priori assumption that $F_1^e$ in the decomposition (8) is symmetric ('sym') determines the orientation of the elastic relaxed placement. Referring to Wallin and Holopainen (2012), the formula

$$F^e = F F_n^{P-1} \exp(-\Delta t \tilde{L}^P), \quad (37)$$

in which $\tilde{L}^P = \bar{D}^P + \tilde{W}^P$, is employed for the update of $F^e$. The tensorial exponent is numerically computed by exploiting the Pade approximation.

According to the stress equilibrium (29), the viscoelastic stretching evolves as

$$d_2^e = \frac{1-D}{\eta}(\mathcal{L}_1^e : \ln v_1^e - \mathcal{L}_2^e : \ln v_2^e). \quad (38)$$

It is assumed that the total deformation (deformation gradient $F$) is priori known, which is the situation in standard finite element codes. Then, the system of nonlinear equations are completed by the integration of the internal rate variables defined by (33), (37), and (38), respectively, i.e.

$$ 
\begin{aligned}
\mathbf{R}_1 &:= \mathbf{F}^{\mathrm{e},\mathrm{T}} - \mathbf{F} \mathbf{F}_n^{\mathrm{P}-1} \exp(-\Delta t \tilde{\mathbf{L}}^{\mathrm{P}}), \\
\mathbf{R}_2 &:= \mathrm{skew}(\mathbf{F}^{\mathrm{e}}), \\
\mathbf{R}_3 &:= \mathrm{sym}(\tilde{\mathbf{W}}^{\mathrm{P}}), \\
\mathbf{R}_4 &:= \mathrm{skew}(\mathbf{F}_1^e), \\
\mathbf{R}_5 &:= \ln v_2^e - \ln v_{2,n}^e - \mathrm{d}(\ln v_2^e)/\mathrm{d}t, \\
\mathbf{R}_6 &:= (s_1 - s_{1,n} - \dot{s}_1 \Delta t)/s_{ss}
\end{aligned} 
\quad (39) $$

where the relation between $d \ln(v_2^e)/dt$ and $d_2^e$ by (38) is determined e.g. in Hoger (1987).

The residuals $R_2$ and $R_4$ consist of only upper (lower) non-diagonal components and the residuals $R_3$ and $R_5$ consist of a maximum of six equations. Solution of (39) is discussed in the Appendix.

## 2.4. Fatigue model

Fatigue failure of polymers is generally due to either mechanically or thermally dominated mechanisms, Lesser (2002); Janssen et al. (2008b). The interest here is placed on the mechanically dominated fatigue, where the material ultimately fractures brittly. Mechanical modes that occur at relatively low stresses and frequencies are characterized by the two-step process described in Introduction.
When dealing with fatigue under complex variable loadings, a suitable dam-
age rule constitutes an integral part of the analysis. Various approaches for
fatigue analysis under such conditions exist and three representative groups
can be identified: energy approach, strain approach, and stress approach, cf.
Wang and Yao (2004); Ottosen et al. (2008); Nittur et al. (2013). The stress
approach, which has generally been used for mechanically dominated, ulti-
mately brittle and high-cycle fatigue, is considered as a basis of the proposed
model. Many of these approaches represents fatigue-limit criteria in which
the fatigue limits are described under an infinite number of identical cycles.
For finite life predictions these models are equipped with cumulative damage
theories, which determine the damage growth per cycle and thus require that
the loading is comprised of well-defined cycles, Lemaitre and Desmorat (2005);
Ottosen et al. (2008). To define equivalent, representative cycles for load his-
tories, cycle-counting methods need to be applied, Fatemi and Yang (1998);
Mortazavian and Fatemi (2015). However, it is often challenging to extract a
standard cycle from a complex load spectrum, which makes the cycle-counting
approaches difficult for demanding practical applications. Another way is to
formulate the model within the continuum damage mechanics (CDM) frame-
work exploiting an incremental formalism without the need to measure dam-
age changes per loading cycles, Lemaitre and Desmorat (2005); Ottosen et al.
(2008); Krairi and Doghri (2014).

An appealing model suitable for describing a mechanically dominated fa-
tigue behavior was proposed by Ottosen et al. (2008). According to this model,
uniaxial and multiaxial stress states for arbitrary loading histories are treated in
a unified manner by using incremental expressions, thus avoiding cycle-counting
methods. Only few macroscopical quantities and a single parameter set are used,
which characteristic makes the model simple and suitable for practical applica-
tions. By exploiting this evolution equations based fatigue modeling concept, a
model for predicting fatigue of amorphous glassy polymers is proposed.

The model relies on the concept of an evolving damage variable and an en-
durance surface which can shift in the space of stress, see Ottosen et al. (2008).
Although the original approach is aimed at the metal fatigue, its formulation is
also suitable for polymers that show similar fatigue characteristics:

- the process of fatigue failure is comprised of three stages ( 1) nucleation and progress of voids and micro-cracks around local inhomogeneities, 2) stable crack growth, 3) unstable cracking leading to failure), cf. Sec. 1;

- defects can nucleate and grow during the service life at stress levels well below the nominal yield strength;

- most of the fatigue life is spent in the 1st stage;

- crack growth during stages 2) and 3) shows a marked resemblance between polymers and metals;

- asymptotic behavior at the extremes of lifetime, i.e. the endurance limit exists and fatigue under this limit is suppressed,

Lugo et al. (2014); Ravi Chandran (2016); Kanters et al. (2016). The model
should also predict a decreasing fatigue life with an increasing particle size, which property is typical for both metals and polymers, Lugo et al. (2014). However, this microstructure-based fatigue modeling is out of the scope of this work. Given fatigue characteristics are the important prerequisites for the proposed model extension.

For many amorphous glassy polymers such an endurance surface can be determined, i.e. the cyclic lifetime increases with a decreasing accumulation of applied stress suggesting a well-defined plateau when ultimate failure finally takes place at finite numbers of cycles just above the endurance limit, Lu and Kim (2007); Lugo et al. (2014). Wöhler curves are commonly used to illustrate these characteristics and to identify polymers' endurance limits. Examples of polymers showing an endurance limit are PC, PMMA, and amorphous polyamide 6,6 (PA-6,6), Lesser (2002); Lu and Kim (2007); Ravi Chandran (2016).

The endurance surface is considered as a function of the stress history and its movement in the stress space is modeled with a reduced deviatoric stress measure that determines the center of the endurance surface. In contrast to the plasticity theories for metals where the endurance surface may lie inside the yield surface, amorphous polymers do not show an explicitly defined yield surface and the fatigue damage development is always induced by the propagation of plastic deformation. Also, since the polymer chains start to align with the loading direction at relatively low stresses and plastic strains, polymer materials show an anisotropic response that is in the model described by the backstress $\beta$, (35).

When reducing loadings, the effect of the backstress decreases and under high-cycle-fatigue (HCF) loadings, the presence of another backstress quantity, $\alpha$ (defined subsequently), is mandatory for governing the fatigue of virtually elastic responses (cf. e.g. steels). Furthermore, many polymers are susceptible to fatigue only once the yield strength is almost reached. In the beginning of such cyclic loadings, plastic deformations and backstress $\beta$ are low, and backstress $\alpha$ dominates fatigue and defines the endurance surface. During continued cyclic loadings, however, plastic deformation may accumulate (ratcheting) and then increases the influence of the backstress $\beta$ in relation to $\alpha$. Therefore, the endurance function,

$$ \beta = \frac{1}{\sigma_0}(\bar{\tau} + aI_1 - \sigma_0), \quad (40) $$

where

$$ \bar{\tau} = \sqrt{3J_2(\tau^{\text{dev}} - \beta^{\text{dev}} - \alpha)} = \sqrt{\frac{3}{2}(\tau^{\text{dev}} - \beta^{\text{dev}} - \alpha) : (\tau^{\text{dev}} - \beta^{\text{dev}} - \alpha)} \quad (41) $$

is proposed. The endurance function also includes the first stress invariant $I_1 = \text{tr} \, \tau$ that reflects the effect of hydrostatic stress $1/3\text{tr} \, \tau$, i.e. the hydrostatic tension enhances the fatigue development while fatigue is suppressed under hydrostatic compression. The parameter $a$ in (40) is considered as positive and dimensionless, and it defines in uniaxial cyclic loading, the slope of the Haigh-diagram. The last parameter $\sigma_0$ corresponds to the endurance limit for infinite fatigue life as mean stress is zero.
<figure><img src="image_5.png" /><figcaption>Figure 4: Alternating stress state for (nonproportional) loading (left). The endurance surface tracks the current stress point (cf. (42)) and then moves between the states A and B (generally not fixed) (right). Peripheries of the endurance surfaces in the initial and final state are highlighted by the dashed and solid lines, respectively. Notation d denotes the increment.</figcaption></figure>

The endurance surface in the deviatoric plane is spherical as illustrated in Fig. 4. As shown, it is the $\alpha + \beta^{\text{dev}}$ tensor that defines the center of the endurance surface. Once an alternating loading is applied, the endurance surface tracks the current stress since the movement of $\alpha + \beta^{\text{dev}}$ always is in the direction of $\tau^{\text{dev}} - \beta^{\text{dev}} - \alpha$, i.e. the $\alpha + \beta^{\text{dev}}$ tensor stores the load history and results in the movement of the endurance surface. It then follows that the evolution of $\alpha$ is governed by a rule

$$ \alpha^{\nabla} = C(\tau^{\text{dev}} - \beta^{\text{dev}} - \alpha)\dot{\beta} - \beta^{\text{dev},\nabla} \quad (42) $$

where

$$ \alpha^{\nabla} = \dot{\alpha} - l\alpha - \alpha l^{\text{T}} $$

is an objective stress rate consistent with the rate of the backstress $\beta^{\text{dev},\nabla}$ and the rate of the deviatoric Kirchhoff stress $\tau^{\text{dev},\nabla}$ defined subsequently in Sec. 2.5.

In (42), a single material parameter C was introduced. Hydrostatic stress influences also the evolution of $\alpha$, which effect is included to (42) through the rate $\dot{\beta}$ of the endurance surface (40).

Since $\alpha$ is considered an overall driving force for fatigue damage, it only evolves once the current stress is crossed the endurance surface and recedes from it, i.e. the condition

$$ \beta \ge 0, \quad \dot{\beta} \ge 0 \quad \Rightarrow \quad \dot{D} \ge 0, \quad \alpha^{\nabla} \neq 0 \quad (43) $$

holds. This situation is demonstrated in Fig. 4.

## Demonstration of the model

Considering uniaxial stress states when the only nonzero stress component is $\tau := \tau_{11}$, while the relation between the backstress components are $\beta^{\text{dev}} :=$
$\beta_{11}^{\text{dev}} = -2\beta_{22}^{\text{dev}} = -2\beta_{33}^{\text{dev}}$ and $\alpha := \alpha_{11} = -2\alpha_{22} = -2\alpha_{33}$. Then, the effective stress (41) takes a simple form as

$$ \bar{\tau} = \kappa(\tau - \frac{3}{2}(\beta^{\text{dev}} + \alpha)) \qquad (44) $$

where the parameter $\kappa$ controls the sign, i.e. $\kappa = 1$ as $(\tau - \frac{3}{2}(\beta^{\text{dev}} + \alpha)) > 0$ and $\kappa = -1$ as $(\tau - \frac{3}{2}(\beta^{\text{dev}} + \alpha)) < 0$. Moreover, the fatigue limit criterion of the endurance surface (40) reduces to a form given as

$$ \beta = \frac{1}{\sigma_0} \left\{ \kappa \left( \tau - \frac{3}{2}(\beta^{\text{dev}} + \alpha) \right) + a\tau - \sigma_0 \right\} = 0, \qquad (45) $$

which form can be expressed as

$$ \tau = \frac{1}{\kappa + a} \left\{ \kappa \frac{3}{2}(\beta^{\text{dev}} + \alpha) + \sigma_0 \right\}, \qquad (46) $$

representing the limit (stress limit) of the endurance surface below which fatigue damage does not accumulate, Lesser (2002).

Under a proportional loading where the stress varies linearly between arbitrary stress states $\tau(A)$ and $\tau(B)$ (stress path pierces the center of the endurance surface, see Fig. 4), one can show that the endurance surface (40) in a high cycle regime ($\beta \to 0$) will reduce to the form given in Sines (1959), see Ottosen et al. (2008, Eqs. 8-11) for a more detailed account. Moreover, when the stress state is uniaxial, the fatigue limit criterion of (40) further reduces to

$$ \tau_a + a\tilde{I}_1 - \sigma_0 = 0 \qquad (47) $$

where $\tau_a$ is the stress amplitude and $\tilde{I}_1 := \frac{1}{2}(I_1(\tau(A)) + I_1(\tau(B)))$ is the mean value of the traces.

Assuming also that the stress ranges periodically between $\tau(A) = \tau_m + \tau_a$ and $\tau(B) = \tau_m - \tau_a$. In such a situation, $\tilde{I}_1 = \tau_m$ in (47), which indicates that the variation between the fatigue stress amplitude and mean stress is linear. Considering further that $\tilde{I}_1$ does not alter, (47) then implicates that the mean shear stress does not effect on the fatigue development. As long as the stress is under the static yield limit and the mean stress is positive or virtually zero, previous observations for many metals, amorphous polymers, and their compounds indicate that the linear relation (47) is in good agreement with experiments, cf. Janssen et al. (2008a, Figure 11a) and Ottosen et al. (2008) with the references therein.

Referring to Fig. 4 and Fig. 5 for an alternating uniaxial loading, the concept is further demonstrated. During loading from state 1 to 2, the stress point recedes from the endurance surface when damage evolves, i.e. $\beta > 0$ and $\dot{\beta} > 0$, cf. (43). Between the states 2 and 3, the stress path has crossed the surface and is in the space within the surface. Then, $\dot{\beta} < 0$ when damage and the backstress do not evolve, i.e. $\dot{D} = 0$ and $\alpha^\nabla = 0$ until the path again pierces the endurance surface at state 3. It then follows that $\alpha_3 = \alpha_2$. From state 3 to 4, damage again evolves since $\beta > 0$ and $\dot{\beta} > 0$. When loading from state 4, the damage development is inhibited until state 5 is reached, i.e. $\alpha_5 = \alpha_4$.
<figure><img src="image_6.png" /><figcaption>Figure 5: Alternating uniaxial loading and demonstration of a damage development during cycling is indicated by a double curve. The time periods for damage developments depend on the applied loading and the endurance surface inherent to the material at hand.</figcaption></figure>

## Damage evolution

Despite possible strain hardening and subsequent directional damage fields especially in large deformations, the fatigue behavior is described by a scalar valued quantity so as to keep the model simple. Assuming damage increases only if $\beta > 0$ and $\dot{\beta} > 0$, an exponential form

$$ \dot{D} = K \exp(f(\beta; L_1, L_2, \vartheta)) \dot{\beta} \geq 0, \quad (48) $$

with the values $0 \le D \le 1$, is valid for the damage evolution law, i.e. it satisfies the dissipation inequality (26). In (48), $K, L_1, L_2$, and $\vartheta$ are positive material parameters and thus, the damage evolution equation (48) shows that damage increases with the distance from the endurance surface.

Many amorphous glassy polymers (such as PC) show only a moderate accumulation of the applied stress as the cyclic lifetime reduces. To capture such a behavior, a function $f$ having two linear asymptotes for positive $\beta$ is defined, i.e.

$$ f(\beta; L_1, L_2, \vartheta) = L_1\beta - L_2 \left[ \beta + \frac{L_2}{\vartheta} \left( \exp(-\vartheta\beta/L_2) - 1 \right) \right], $$

which has the asymptote $L_1\beta$ when $\beta \to 0$ (HCF-regime) and $(L_1 - L_2)\beta$ when $\beta$ is large (LCF-regime). The curvature $\vartheta$ determines how rapidly the second asymptote is reached.

Let $x, t_0$, and $t_c$ be the material placement, initial time instant, and the final, critical time instant, respectively. Then, $D = D(x, t_0) := D_0 \ge 0$ is an initial value, while $D(x, t_c) = D_c = 1$ represents a local fatigue failure right prior to a rapid progress of small cracks to a form of large cracks (rupture).

## 2.5. Numerical treatment of the fatigue model

To integrate the evolution equations (42) and (48), the function (40) for $\beta$ needs to be differentiated. Assuming that the current deformation gradient $F$
and an increment for the strain tensor $\Delta\epsilon$ can be extracted during incremental calculations. The increment of the strain tensor relies on the relation (3), i.e.

$$ \Delta\epsilon = \frac{1}{2}(\Delta l + \Delta l^T) \qquad (49) $$

which defines the increment $\Delta F := \Delta t \dot{F}$ as

$$ \Delta l := \Delta F F^{-1}. \qquad (50) $$

Then, based on (1) and (3) one obtains

$$ \Delta l := \Delta l^e + F^{e-T} \Delta \bar{L}^p F^{e-1} \qquad (51) $$

where

$$ \Delta l^e := \Delta F^e F^{e-1}, \quad \text{and} \quad \Delta \bar{L}^p := \Delta t (\bar{D}^p + \tilde{W}^p) \qquad (52) $$

are the incremental presentations for the elastic and plastic velocity gradient, respectively. Owing to the symmetry of $\Delta F^e$ finally at the integration interval,

$$ \Delta F^e = \frac{1}{2}((\Delta\epsilon - F^{e-T}\Delta\bar{L}^p F^{e-1})F^e + F^e(\Delta\epsilon - F^{e-T}\Delta\bar{L}^{p,T}F^{e-1})). \qquad (53) $$

Substituting (53) into (52) and then the obtained result into (51) gives the expressions for $\Delta l^e$ and $\Delta l$, respectively.

The differentiation of $\beta$ given by (40) yields

$$ \Delta\beta = \frac{1}{\sigma_0 + C\bar{\tau}} \left[ \frac{3}{2} \frac{(\tau^{\text{dev}} - \beta^{\text{dev}} - \alpha)}{\bar{\tau}} : \Delta\tau^{\text{dev},\nabla} + a \text{trace}(\Delta\tau^{\nabla}) \right] \qquad (54) $$

where $\Delta\beta := \Delta t \dot{\beta}$ and $\Delta\tau^{\nabla} := \Delta t \tau^{\nabla}$ is an objective stress increment defined below.

According to (42),

$$ \Delta\alpha^{\nabla} = C(\tau^{\text{dev}} - \beta^{\text{dev}} - \alpha)\Delta\beta - \Delta\beta^{\text{dev},\nabla} \qquad (55) $$

where $\Delta\alpha^{\nabla} = \Delta t \alpha^{\nabla}$,

$$ \Delta\beta^{\nabla} := \Delta t \beta^{\nabla} = F \Delta \left( F^{-1} \beta F^{-T} \right) F^T = \Delta\beta - \Delta l \beta - \beta \Delta l^T, \qquad (56) $$

and

$$ \Delta\beta := F^e \Delta \bar{B}^{\nabla} F^{e,T} + \Delta l^e \beta + \beta \Delta l^{e,T}. \qquad (57) $$

The increment $\Delta \bar{B}^{\nabla}$ of the backstress in the relaxed placement is obtained by taking the increment of (35). The inverse Langevin function included in $\Delta \bar{B}^{\nabla}$ is approximated by the Pade-approximation when its derivative with respect to its argument is given by

$$ \xi_{ec} := \frac{d}{dX} \mathcal{L}^{-1}(X) \bigg|_{X=\frac{\lambda_{ec}^p}{\sqrt{N}}} = \frac{\left(\frac{\lambda_{ec}^p}{\sqrt{N}}\right)^4 + 3}{\left(1 - \left(\frac{\lambda_{ec}^p}{\sqrt{N}}\right)^2\right)^2}. \qquad (58) $$
Differentiation of (34) yields

$$ \frac{\partial \lambda_{ec}^{p}}{\partial \bar{C}^{p}} = \frac{1}{6 \lambda_{ec}^{p}} \bar{I} \qquad (59) $$

and its substitution into the incremental expression of (35) leads to the form

$$ \Delta \bar{B}^{\nabla} = \frac{C_R \sqrt{N}}{3 \lambda_{ec}^{p}} \mathcal{L}^{-1} \left( \frac{\lambda_{ec}^{p}}{\sqrt{N}} \right) \left[ \left( \frac{\chi_{ec} \xi_{ec}}{3 \sqrt{N} \lambda_{ec}^{p}} - \frac{1}{3 (\lambda_{ec}^{p})^2} \right) \text{trace}(\Delta \bar{L}^{p} \bar{C}^{p}) \bar{C}^{p} + 2 \text{sym}(\Delta \bar{L}^{p} \bar{C}^{p}) \right] \qquad (60) $$

where $\chi_{ec} := 1/\mathcal{L}^{-1}(\lambda_{ec}^{p}/\sqrt{N})$.

The non-linear system of equations for fatigue are completed by the integration of the internal rate variables defined by (42) and (48), i.e. $R_{\alpha} = [R_{\alpha,1} \quad R_{\alpha,2}]$ where

$$ \begin{aligned} R_{\alpha,1} &:= \alpha - \alpha_n - \Delta t \alpha^{\nabla}, \\ R_{\alpha,2} &:= D - D_n - \Delta t \dot{D}. \end{aligned} \qquad (61) $$

The treatment of the nonlinear systems given in (39) and (61) is discussed in the Appendix.

### 2.5.1. Stress increment

A glance at (49) reveals that $\Delta l = \Delta \epsilon + \tilde{\omega}$ when the stress increment becomes

$$ \Delta \boldsymbol{\tau}^{\nabla} = \boldsymbol{F} \Delta (\boldsymbol{F}^{-1} \boldsymbol{\tau} \boldsymbol{F}^{-T}) \boldsymbol{F}^{T} = \Delta \boldsymbol{\tau} - \tilde{\omega} \boldsymbol{\tau} - \boldsymbol{\tau} \tilde{\omega}^{\mathrm{T}} = \boldsymbol{c} : \Delta \boldsymbol{\epsilon} \qquad (62) $$

where

$$ \boldsymbol{c} := \boldsymbol{c}^{\epsilon} - \boldsymbol{c}^{\omega} := \frac{\mathrm{d}\boldsymbol{\tau}}{\mathrm{d}\boldsymbol{F}} : (\boldsymbol{i} \odot \boldsymbol{F}) - (\boldsymbol{\tau} \odot \boldsymbol{i} + \boldsymbol{i} \odot \boldsymbol{\tau}) - \boldsymbol{c}^{\omega} \qquad (63) $$

is the algorithmic tangent stiffness (ATS) tensor in which $\boldsymbol{c}^{\epsilon} : \Delta \boldsymbol{\epsilon} := \Delta \boldsymbol{\tau}$ and $\boldsymbol{c}^{\omega} : \Delta \boldsymbol{\epsilon} := \tilde{\omega} \boldsymbol{\tau} + \boldsymbol{\tau} \tilde{\omega}^{\mathrm{T}}$. The non-standard tensor product in (63) is defined as

$$ (\boldsymbol{A} \odot \boldsymbol{B})_{ijkl} = 1/2 (\boldsymbol{A}_{il} \boldsymbol{B}_{jk} + \boldsymbol{A}_{ik} \boldsymbol{B}_{jl}). $$

Due to the imposed symmetry of $\boldsymbol{F}^{\epsilon}$, the elastic spin vanishes, and thus, $\tilde{\omega} := \boldsymbol{F}^{\mathrm{e}-\mathrm{T}} \tilde{\boldsymbol{W}}^{\mathrm{p}} \boldsymbol{F}^{\mathrm{e}-1}$ in (62). In view of (29),

$$ \frac{\mathrm{d}\boldsymbol{\tau}}{\mathrm{d}\boldsymbol{F}} = \mathcal{L}_{1}^{\mathrm{e}} : \left( (1-D) \frac{\partial(\ln \boldsymbol{v}_{1}^{\mathrm{e}})}{\partial \boldsymbol{v}_{1}^{\mathrm{e}}} : \frac{\mathrm{d}\boldsymbol{v}_{1}^{\mathrm{e}}}{\mathrm{d}\boldsymbol{F}} - \ln \boldsymbol{v}_{1}^{\mathrm{e}} \otimes \frac{\mathrm{d}D}{\mathrm{d}\boldsymbol{F}} \right) \qquad (64) $$

where the derivative of the tensor logarithm is determined following Miehe (1998). The details for obtaining the terms $\mathrm{d}\boldsymbol{v}_{1}^{\mathrm{e}}/\mathrm{d}\boldsymbol{F}$ and $\mathrm{d}D/\mathrm{d}\boldsymbol{F}$ set out in the Appendix.

The stress increment (62) represents an appealing objective time derivative of the Kirchhoff stress termed the convected rate. In user-defined material model subroutines of commercial finite-element packages, the objective stress rate and the corresponding stress increment required in the equilibrium iterations may differ from that presented in (62). A usual one is the Jaumann rate, $\boldsymbol{\sigma}^{\nabla J}$, and its incremental representation based on (62) is given by

$$ \Delta \boldsymbol{\sigma}^{\nabla J} := \Delta t \boldsymbol{\sigma}^{\nabla J} = J^{\mathrm{e}-1} \Delta \boldsymbol{\tau}^{\nabla} + \Delta \boldsymbol{\epsilon} \odot \boldsymbol{\sigma} + \boldsymbol{\sigma} \odot \Delta \boldsymbol{\epsilon} - \text{trace}(\Delta \boldsymbol{\epsilon}) \boldsymbol{\sigma}. \qquad (65) $$
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Testing arrangement (right) and shots of a dogbone-shaped specimen at $u = 10$ mm, 20 mm, 30 mm, and 40 mm (left). The displacement at $y = 0$ is restrained. The details of the testing arrangement is defined in ASTM-D638 and ISO 527-2.</figcaption></figure>

# 3. Applications and discussion

## 3.1. Constitutive model parameters

The model introduced in Section 2.2 was first calibrated to data obtained from cold drawing tests on dogbone-shaped PC-specimens. Fig. 6 shows the specimens' geometry. The un-notched test specimens comply with ISO 527 and were fabricated from Lexan® 223R PC-polymer by injection-molding. The tests were conducted by the Instron® 5967 machine controlled by Blue Hill 3 software. During drawing, the applied load $f$ and the elongation $u$ were monitored. To prevent a temperature rise owing to the dissipative heating, a low crosshead speed $\dot{u} = 2$ mm per a minute was applied. Repeated unloadings to varying stress levels were conducted at an equable rate 1 kN/min. For the implementation, a finite-element method is applied, and the specimen was meshed with linear solid elements. Resulting responses and the observed data for overall $f-u$-responses are shown in Fig. 7. A more detailed account for the test program involving repeated loading cycles is found from Holopainen (2013). The calibrated parameters based on the least squares fitting are presented in Table 1.

### Uniaxial compression - non-monotonic loading

As an introduction to fatigue, let us consider a simple, uniaxial non-monotonic compression case involving few cycles. Anand and Ames (2006) conducted such cyclic as well as monotonic tests on PMMA. The proposed model was also calibrated to those tests under monotonic loadings. Figure 8 shows the model predictions for uniaxial cyclic compression loading in which the loading cycles comprise the straining up until the logarithmic strain of 20% is reached and a subsequent unloading to zero strain.

The development of the backstress is also shown. The model predicts increasing backstress during both the loading and unloading, which then smoothly reduces during the subsequent period to zero strain. Therefore, in addition to the
<figure><img src="image_9.png" /><figcaption>Figure 7: Overall *f* – *u* responses predicted by the model with calibrated material parameters and measured from the experiment. The unloadings are performed to *f* = 60 N and *f* = 1000 N followed by the dwell period of *t* = 120 s.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: The compressive true stress vs logarithmic strain response under cyclic loading for PMMA (left). The solid line and dash-and-dot line curve with markers are the model response and experimental data taken from Anand and Ames (2006). The undermost loop shows the stabilized behavior. The plot in right shows the evolution of the backstress $\beta^{dev}$. The calibrated material parameters for PMMA are $E = 2300$ MPa, $s_0 = 121$ MPa, $s_{ss} = 63.3$ MPa, $h_1 = 184$ MPa, and $\dot{\gamma}_0 = 5.6 \cdot 10^{15} s^{-1}$. The rest are equal with those given in Table 1.</figcaption></figure>

viscoelastic constitutive treatment, the model predicts nonlinear stress-strain relationship during both the loading and unloading followed by a continued elongation at a constant stress termed dwelling, Wallin and Holopainen (2012). In contrast to a purely elastic deformation, viscous deformation results in a loss (dissipation) of mechanical energy which is equal to the area of the hysteresis loops through loading cycles. Fig. 8 also shows that the model predicts the Bauchinger effect, which stabilizes during few cycles leading to the saturated
state of hardening when the predicted stress-strain response matches well with the experimentally measured response.

## 3.2. Fatigue model parameters

Once the constitutive model parameters were found, the parameters for modeling fatigue were determined from in situ measurements taken from Janssen et al. (2008a) and Janssen et al. (2008b). The material employed in the tensile fatigue test is a quenched PC (Lexan® 101R and 161R), and the specimen geometry used for fatigue studies is a common dogbone-shaped, injection molded tensile specimen (ASTM D638-IV), Foster (2015). Isothermal test conditions for mechanically induced fatigue have included uniaxial stress submitted to load control at the room temperature of 23°C. Referring to Janssen et al. (2008b), the fatigue failure occurs in the gauge section. The test specimen is also termed a uniaxial tensile bar, and it will be shown later that the stress state in that region is virtually uniaxial, i.e. the only nonzero stress component employed in the calibration is given by

$$ \sigma = \sigma_m + \sigma_a \sin(\omega t) \qquad (66) $$

representing an average true stress somewhere in the gauge section. In (66), $\omega$ denotes the frequency, $t$ is time, $\sigma_m$ is mean stress, and $\sigma_a$ is the stress amplitude.

Since fatigue tests are both costly and time consuming, an accelerated life testing with a high mean stress has been applied, Lesser (2002); Maxwell et al. (2005); Lu and Kim (2007); Zhang et al. (2015). In the referred test, see Janssen et al. (2008b), the stress waveform is sinusoidal, and the maximum stress is allowed

<table><caption>Table 1: Constitutive model parameters for PC fitted to data taken from the cold drawing experiments on a dogbone-shaped test specimens. The remaining viscoelastic constitutive parameters are $E = 2000$ MPa, $E_1 = 1000$ MPa, $\eta = 1.5 \cdot 10^4$ MPas, and $\nu = 0.37$.</caption><thead><tr><th>Parameter</th><th>s<sub>0</sub></th><th>s<sub>ss</sub></th><th>h<sub>1</sub></th><th>γ̇<sub>0</sub></th><th>A</th><th>C<sup>R</sup></th><th>N</th><th>α</th></tr></thead><tbody><tr><td>Unit ........</td><td>MPa</td><td>MPa</td><td>MPa</td><td>s<sup>-1</sup></td><td>MPa<sup>-1</sup>K</td><td>MPa</td><td></td><td></td></tr><tr><td>Value ........</td><td>96</td><td>76</td><td>720</td><td>5.6 · 10<sup>15</sup></td><td>240</td><td>14</td><td>2.2</td><td>0.08</td></tr></tbody></table>

<table><caption>Table 2: Fatigue model parameters for PC. The fit was performed to experimentally observed fatigue data given in Janssen et al. (2008a) and Janssen et al. (2008b).</caption><thead><tr><th>Parameter</th><th>σ<sub>0</sub></th><th>a</th><th>C</th><th>K · 10<sup>5</sup></th><th>L<sub>1</sub></th><th>L<sub>2</sub></th><th>ϑ</th></tr></thead><tbody><tr><td>Unit ........</td><td>MPa</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Value ........</td><td>34.0</td><td>0.46</td><td>0.05</td><td>8.50</td><td>7.60</td><td>3.80</td><td>2.63</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 9: Fatigue strengths of PC under the application of the high mean stress $\sigma_m = \sigma_a + 2.2$ MPa (left). The solid line denotes the model result, and the marker ▲ refers to data points taken from Janssen et al. (2008b) (Fig. 6b). The upper and lower horizontal dashed lines refer to the static tensile yield strength $\sigma_y$ and an estimated endurance limit, respectively. Haigh-diagram for the endurance limit according to the model (right). The parameter $a$ determines the slope in the diagram, and the dashed curve represents the expected experimental behavior. A tangent line for the loading $\sigma_{max} = \sigma_y$ at $\sigma_m = 35$ MPa is also represented.</figcaption></figure>

to vary, while the minimum stress is kept at 2.2 MPa. To guarantee isothermal testing conditions, a low frequency of 2 Hz has been used.

First, the parameter $a$, which is the slope of the Haigh-diagram near zero mean stress, was determined, see the endurance function (40). Unfortunately, data which account for this region of the amplitude-mean stress relationship are not available, and the value of $a$ is fitted to an expected experimental curve shown in Fig. 9 (right). For many materials, however, the slope of the Haigh-diagram is not constant, but reduces towards the zero mean stress, Ottosen et al. (2008). The slope of the curve at high mean stresses is extracted from data given in Janssen et al. (2008a) (Fig. 11a therein). The data show observed mean stresses 55-70 MPa for a wide range of fatigue lives as two low stress amplitudes of 5 MPa and 10 MPa have been applied. In this amplitude region, the relationship between the mean stress and amplitude appears to be unaltered, i.e. the addition of 5 MPa in the stress amplitude reduces the mean stress virtually the same degree for all observed fatigue lives. Thus, $a$ is almost unity. The predicted Haigh-diagram is also shown in Fig. 9 (right). According to (47), the variation between mean stress and the fatigue strength is linear, and the slope of the line is given in terms of the parameter $a$.

The fatigue limit $\sigma_0$ for zero mean stress and the remaining parameters $C$, $K$, $L_1$, $L_2$, and $\vartheta$ were calibrated to data shown in Fig. 9 by using the concept analogous to that given in Ottosen et al. (2008) (section 5 therein). According to this concept, the equations for the backstress and damage are first integrated in order to define the damage accumulation per a loading cycle and the number of cycles $n$ needed to reach fatigue failure. This expression for $n$ is then used in a single object function to be minimized. The least squares fitting, which utilizes the Nelder-Mead simplex algorithm, is employed. The least squares error of number of cycles $n$ that leads to fatigue failure is chosen as the object function,
<figure><img src="image_12.png" /><figcaption>Figure 10: Lifetime of PC using R-values of 0.7, 0.4, and 0.1 (left). The solid, dashed, and dotted curves represent the corresponding responses, specified by the thick black and thin red curves for the model results and data fittings taken from Kanters et al. (2016) (Lexan® 143R). Sinusoidal waveform with the frequency of 1 Hz has been applied. Fatigue lifetimes under uniaxial compression with $\sigma_{max} = 2.2$ MPa (right). The model result and data fitting taken from Janssen et al. (2008b) are highlighted by the solid and dashed curve, respectively.</figcaption></figure>

i.e.

$$ \min F := \sum_{i=1}^{\hat{I}} (w^{(i)})^2 \left( 1 - \frac{\ln(n^{(i)})}{\ln(n_{em}^{(i)})} \right)^2 \quad (67) $$

where $\hat{I}$ denotes the count of data points, $w^{(i)}$ are the weights, and $n^{(i)}$ and $n_{em}^{(i)}$ denote the predicted and recorded number of cycles, respectively. All data points shown in Fig. 9 are used in the calibration.

The calibrated parameters and model responses are depicted in Table 2 and Fig. 9, respectively. Although the model governs only the crack initiation step, inaccuracy of the model response is small since this step typically covers over 95% of the total fatigue life, Marissen et al. (2001); Janssen et al. (2008a). The results indicate that the cyclic lifetime smoothly decreases with increasing accumulation of applied stress. Reducing the stress level, a transition in the failure mode takes place from ductile to brittle, and ultimate failure can finally be reached after finite numbers of cycles at the stress level right above the expected endurance limit.

### 3.3. Results - uniaxial stress state

Based on the calibrated parameters in Tables 1 and 2, time-to-failure versus maximum stress is shown in Fig. 10. Due to a slightly different PC used in the tests, a small distinction between the model responses and data can be observed at low fatigue lifetimes, especially under the low stress ratio $R$ (the ratio between the smallest and largest stresses). The bigger $R$ is, the lower is the maximum
<figure><img src="image_13.png" /><figcaption>Figure 11: The cyclic true stress (solid line) and the endurance surface (limits are highlighted by the dashed and dash-and-dot lines) during the first few stress cycles (left). The damage development during the first two cycles is indicated by the double curve. Development of damage variable according to (48) (right).</figcaption></figure>

stress for a short time-to-failure whereas curves show opposite trend at high fatigue lifetimes. Moreover, the difference between the curves seems to reduce with increasing lifetimes, cf. Kanters et al. (2016), i.e. the responses show an asymptotic behavior at the extremes of lifetime and almost equal endurance limits.

Some data for the fatigue strengths in compression are also available, cf. Janssen et al. (2008b). Under such conditions (a notable negative mean stress), the linear relationship between the effective stress and the hydrostatic stress in the expression of the endurance surface (40) is unsuitable, and the slope parameter *a* of the Haigh-diagram needs to be modified, see Fig. 9. Here, trivially a negative value *a* = -0.17 is used. Despite the calibration process only for tension, this elementary modification yields the model result that reflects well the available experimental data for a heavy compression, see Fig. 10(right). The results show that the fatigue strengths for the compression loading are substantially greater than under tension.

The predicted endurance surface (46) under a sinusoidal cyclic tension is shown in Fig. 11. Once the periodic loading is applied, the endurance surface tracks the current stress point and reaches its periodic state since the movement of $\alpha$ always is in the direction of $\tau^{\text{dev}} - \beta^{\text{dev}} - \alpha$, cf. equation (42) along with the plots in Fig. 4. During continued cycling, damage always develops as the stress state recedes from the endurance surface, i.e the upper limit of the endurance surface in Fig. 11 shows an increase and the stress is greater than this limit. A strong initiation of damage is due to the high mean stress applied, see Fig. 11(right).

Figure 12 further shows a periodically increasing development of the logarithmic strain in the beginning of cyclic tension. Since the stress always remains below the yield strength of the present polymer material, plastic deformation is initially low and the total strain and its elastic plus viscoelastic component are virtually overlapping at this phase. In tensile test specimens, the growth of
<figure><img src="image_14.png" /><figcaption>Figure 12: Periodic development of the logarithmic strain in the beginning of cyclic tension shown in Fig. 11.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Lifetime of PC under static (creep) loading $\sigma = \sigma_m$. The solid line is the model result, and the markers △ are data points taken from Janssen et al. (2008b) (Fig. 3a).</figcaption></figure>

viscoelastic and plastic strains leads to a notable elongation that may cause a plastic instability termed necking and is finally followed by a brittle rupture as shown in Janssen et al. (2008b).

The proposed approach also provides time-to-failure predictions under fixed stress levels ($\sigma = \sigma_m$) in a reasonable accuracy, see. Fig. 13. At room temperature, aging kinetics has not influence on lifetimes in this regime, and the experimental results show almost a linear response. The model tends to overestimate higher lifetimes which characteristic is due to a rise of the endurance surface at mean stresses far above the ones applied in the calibration.

Almost a fully reversed uniaxial loading case consisting of cycling in the frequency of 2 Hz is also studied. The maximum stress up almost to the tensile yield strength is applied. Under this loading, damage increases rapidly leading to an ultra-low fatigue life as depicted in Fig. 15, cf. Drozdov (2011); Del Vecchio et al. (2014). Once the material response has stabilized, hysteresis
<figure><img src="image_16.png" /><figcaption>Figure 14: The cyclic true stress (solid line) and the periodic movement of the endurance surface (dashed and dash-and-dot lines) (left). Model prediction for the hysteresis loops (right).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15: Periodic development of the compressive strain (dashed line) and its plastic component $\epsilon^P$ (solid line) with a low mean stress and a high stress amplitude (left). Damage development up to the final failure (right).</figcaption></figure>

loops demonstrated in Fig. 14 show a shift and Bauchinger effect similar to that observed in metals. A slightly reducing area of the loops indicates a progressive increase in energy-dissipation and a suppression of the irreversible work during evolving fatigue, cf. Lesser (2002); Janssen et al. (2008b); Drozdov (2011). This result is a consequence of neglected temperature effects and physical aging, both of which being relevant assumptions under the applied low frequency loadings.

The shift of hysteresis loops and the development of the plastic strain shown in Fig. 15 refer to ratcheting (cyclic creep), i.e. constantly accumulated plastic strain develops as the cyclic loading continues, cf. Lemaitre and Desmorat (2005); Kang (2008); Del Vecchio et al. (2014); Xi et al. (2015). The model is able to simulate ratcheting since the inverse Langevin function (35) accounts for the extensibility of polymer chains generating an asymptotically growing backstress and strain hardening as the plastic stretches increase.

To further demonstrate the minute ratcheting behavior of the PC applied, similar to that used in Xiong et al. (2014); Xi et al. (2015), the uniaxial ratch-
<figure><img src="image_18.png" /><figcaption>Figure 16: Predictions of ratcheting strain (absolute values, only the response with the high mean stress is positive) under the stresses presented in Figs. 11 and 14. Data curve (blue color) for a PC-polymer is taken from Lu et al. (2016) (left). The numerical values denote the applied mean stress and amplitude. Creep as the stress is kept as a constant after 30 seconds (at σ = 10 MPa or σ = 2.2 MPa) and after 800 seconds (σ = 30.6 MPa or σ = 2.2 MPa) (right).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17: Prediction of the ratcheting strain responses (black lines) with different mean stresses at the frequency of ω = 1 Hz and at room temperature (left). Corresponding data curves (blue color) for PC-7030PJ polymer are taken from Xi et al. (2015). Development of fatigue damage (right).</figcaption></figure>

eting strain is defined as a mean strain,

$$ \epsilon_r = \frac{\epsilon_{\max} + \epsilon_{\min}}{2} \qquad (68) $$

where $\epsilon_{\max}$ and $\epsilon_{\min}$ denote the maximum and minimum logarithmic strain in each cycle, respectively.

The characteristics of ratcheting strain development under different loadings are shown in Figs. 16, 17, and 18. Three different stages can be observed. In the first stage, ratcheting strain accumulates rapidly, but then the growth rate reduces to a nearly constant value. The accumulation of ratcheting strain and the duration of this stage depend strongly on the applied stress level and appear more evidently under low stresses. In the second stage, ratcheting strain shows a steady accumulative region which occupies most of the fatigue life, i.e. the
<figure><img src="image_20.png" /><figcaption>Figure 18: Prediction of the ratcheting strain responses with different amplitudes at the frequency of $\omega = 1$ Hz and at room temperature (left). Development of the fatigue damage (right).</figcaption></figure>

duration of this stage reflects the entire fatigue life. The third stage occurs solely at high stress levels during which the cyclic creep growth rate shows a marked increase prior to ultimate failure. The predicted accumulative behavior of ratcheting strain well corresponds with experimentally observed behavior of amorphous polymers, cf. Kang (2008); Liu et al. (2008); Xi et al. (2015).

Fig. 16(right) then shows the strain behavior as the stress is kept as constant after 30 seconds or 800 seconds. Under the long hold (dwell) time, the response shows evidently (static) creep while neither recovery nor creep occur after 800 seconds. This behavior is in line with experimental results, i.e. the appearance of creep progressively vanishes during loading cycles, see Dreistadt et al. (2009) for a more detailed account. Since the applied constant stress levels are low, fatigue damage becomes suppressed and lifetime explodes when compared to the original cyclic loading.

The impact of mean stress on the ratcheting behavior is shown in Fig. 17. The cyclic creep accumulation rate is initially high but rapidly reduces to almost zero. The greater the mean stress is, the more rapidly the ratcheting strain increases and finally reaches asymptotically its constant state. Experimental results show similar characteristics, see Figs. 16, 17, and Xi et al. (2015). When also compared to the observed lifetimes under the high creep loadings present in Fig. 13, it is not a surprise, that the fatigue development is suppressed under the applied low stress amplitudes, see Fig. 17(right).

Numerical investigations also indicated that the stress amplitudes ranging between $\sigma_a$ = 10...30 MPa have not a marked influence on ratcheting under the applied mean stresses shown in Fig. 17. It means that ratcheting is more prone to the applied mean stress than to the amplitude, cf. Xi et al. (2015). However, when subjected to low mean stresses, substantially below the stress amplitudes, the development of the ratcheting strain is greatly influenced by the stress amplitude, see Figs. 16 and 18.

An interesting observation is that the damage evolution resembles closely the development of ratcheting strain, see Figs 17 and 18, i.e. ratcheting strain can be exploited in the evaluation of the entire fatigue lifetime.
<figure><img src="image_21.png" /><figcaption>Figure 19: Damage development in a test specimen right after 140 cycles (left). Prediction of the ratcheting strain $\epsilon_r$ in the middle of the specimen (thin black curve) (top right). For a comparison, also the predicted ratcheting strain under homogeneous deformation (uniaxial stress) is shown (thick black curve). Corresponding damage developments represented on a logarithmic scale (bottom right).</figcaption></figure>

### 3.4. Multi-axial stress state - fatigue of a dogbone-shaped specimen

Due to the localization of the plastic deformation, damage grows unevenly, leading to a reduced fatigue life somewhere in the material. Thus, it is of interest to investigate the fatigue damage development of the entire test specimen by using a finite-element method. The same test specimen's geometry and finite element mesh are applied as previously introduced in Section 3.1, see Fig. 6. In accordance with the fatigue testing discussed in Section 3.2, the displacement at $y = 0$ is fastened by a clamp and the evenly distributed loading at $y = L$ is controlled by a sinusoidal waveform of the frequency 2 Hz and the mean and amplitude stress 14.85 MPa and 13.75 MPa, respectively. The corresponding average stress state in the gauge section is twofold. Under fatigue loadings, the advantage of adaptive concepts cannot be benefited, and time steps of 0.01 s are used to accurately capture the sinusoidal fatigue signal. Obviously, large computation times can be expected already in LCF-regime.

In the analysis, the specimen does not suffer from impurities or imperfections, which is the case in practise, when the fatigue damage progresses most intensively between the gauge section and clamps, see Fig. 19. Since experiments show final failure somewhere in the gauge section, a small notch or the like would be required to simulate tests more realistically. Fatigue damage also starts with the formation of an oblate region in the middle of the specimen. Li et al. (1995) reported agglomeration of similar crazing zones in PC under
cyclic tension. Although not shown here, the damage development in the specimen closely follows the fields of increasing localized plastic deformation. The result implicates that the localized plastic deformation is a precursor that provides crazing initiation sites which then control the fatigue propagation in the material, cf. Li et al. (1995); Lesser (2002); Lugo et al. (2014).

The predicted development of ratcheting strain in the middle of the specimen is also shown in Fig. 19. Ratcheting strain has been calculated approximatively by dividing the overall displacement *u* by the grip-to-grip distance *H*, see Fig. 6, and then applying (68). It appears that the multiaxial stress path causes an additional hardening that makes ratchetting strain lower and its development more linear compared to the development under the uniaxial stress state. This characteristic is also observed in experimental testing, see Lu et al. (2016). In addition, the development of the plastic deformation under cyclic loadings (ratcheting) is susceptible to the specimens geometry.

In the gauge section, the axial (longitudinal) stress component is orders of magnitude greater than others, i.e. a uniaxial stress state is representative as it was presumed in the calibration. Depending on the load level, specimen will show a rapid rupture after neck initiation (ultra low cycle fatigue) or brittle fracture after hundreds or tens of thousands of cycles, (Janssen et al., 2008b, Fig. 4). Since ratcheting strain becomes suppressed during the applied loading, the damage distribution shown in Fig. 19 represents the latter case. The difference between the damage responses seems to be a consequence of the ratcheting strain accumulation which initially is substantially greater under the uniaxial stress state. However, once ratcheting strain has increased much enough, damage growth under the multiaxial stress shows a tendency towards the lifetime predicted under the uniaxial stress state.

4. Conclusions

High cycle fatigue approach for metals present in Ottosen et al. (2008) was extended for the modeling of fatigue in amorphous polymers. The proposed approach covers mechanically dominated fatigue behavior intrinsic to polymers of amorphous classes, and is formulated in a continuum mechanics framework by using evolution equations that make the definition of damage changes per cycle redundant, i.e. cycle-counting techniques do not need to be applied. The approach takes into account of the influence of plastic deformations on fatigue, and the development of plastic deformations is based on the celebrated 8-chain Boyce-Parks-Argon (BPA) model where the description of the transformation from a single chain deformation to the continuum deformation relies on chain stretching due to the plastic deformation. Since the relaxation time of the polymer chains to attain their equilibrium state after deformation is also affected by the viscoelastic effects, the BPA model was enriched by a viscoelastic element. The approach was calibrated to data taken from both non-monotonic cold drawing experiments and accelerated isothermal fatigue tests on dogbone-shaped PC-specimens.
Ratcheting behavior of polycarbonate (PC) was investigated and a comparison was made with fatigue tests at various mean stresses and amplitudes. Three regions of ratcheting strain could be observed. In the first stage, ratcheting strain accumulates rapidly, but then the growth rate reduces to a nearly constant value. In the second stage, ratcheting strain shows a steady accumulative region which occupies most of the fatigue life, i.e. ratcheting behavior can be exploited in the evaluation of the fatigue lifetime. The third stage occurs solely at high stress levels during which the cyclic creep growth rate shows a marked increase prior to ultimate failure. In each stage, ratcheting strain increases with both the applied mean stress and amplitude, although it is more prone to mean stress.

Also, finite element studies of a dogbone-shaped test specimen were performed for analyzing a compromise between inhomogeneous plastic deformations and fatigue in PC. The results show that

- fatigue damage increases at the sites of localized plastic deformation, which behavior is microscopically manifested by crazing;

- an apparent ratcheting deformation is observed during the stress-controlled cyclic loading;

- multiaxial stress path results in an additional hardening that makes the ratchetting strain lower and its development more linear compared to the development under uniaxial stress;

- damage development obeys the increase of the ratcheting strain closely, i.e. the ratcheting strain is an important indicator of fatigue.

Validated by the experiments, the proposed approach was able to predict fatigue characteristics of PC under different loadings. The results indicate that ratcheting obviously contributes to fatigue and it can be exploited in evaluation of the entire fatigue life of amorphous polymers. To ensure practical significance, the approach needs a low number of material parameters and thus, is well suitable for numerical solution methods.

However, the proposed approach could be further extended. For instance, in low cycle regime, the influence of energy dissipation and subsequent rise in temperature increase when a temperature-dependent material behavior should be considered in the modeling. In the literature, an accelerated (uniaxial) fatigue testing of a high mean stress has been paid considerable attention, and more data that also account for directional fatigue properties, multiaxial ratcheting, and rate dependence of amorphous polymers under various loading situations are needed. One limitation with regard to the fatigue damage modeling on the micro-level is the presumption of material homogeneity, while real polymers show a diffuse microstructure. Moreover, several details for fatigue damage development are quite ambiguous at the moment. Important supplements to be involved in the modeling are certain morphological and microstructural changes, such as crazing and disentanglement of a chain network.
References

Abdel-Karim, M., 2005. Shakedown of complex structures according to various hardening rules. International Journal of Pressure Vessels and Piping 82, 427-458.

Anand, L., Ames, N.M., 2006. On modeling the micro-indentation response of an amorphous polymer. Int. J. Plasticity 22, 1123-1170.

Boyce, M.C., Parks, D.M., Argon, A.S., 1988. Large inelastic deformation of glassy polymers, Part I: Rate-dependent constitutive model. Mechanics of Materials 7, 15-33.

Dreistadt, C., Bonnet, A.S., Chevrier, P., Lipinski, P., 2009. Experimental study of the polycarbonate behaviour during complex loadings and comparison with the Boyce, Parks and Argon model predictions. Materials and Design 30, 3126-3140.

Drozdov, A.D., 2011. Cyclic viscoelastoplasticity and low-cycle fatigue of polymer composites. Int. J. Solids and Structures 48, 2026-2040.

Fang, Q.Z., Wang, T.J., Beom, H.G., Li, H.M., 2008. Effect of cyclic loading on tensile properties of PC and PC/ABS. Polymer Degradation and Stability 93, 1422-32.

Fatemi, A., Yang, L., 1998. Cumulative fatigue damage and life prediction theories: a survey of the state of the art for homogeneous materials. International Journal of Fatigue 20, 9-34.

Foster, A.M., 2015. Materials Testing Standards for Additive Manufacturing of Polymer Materials: State of the Art and Standards Applicability. Technical Report. Materials and Structural Systems Division Engineering Laboratory, NISTIR 8059.

Hoger, A., 1987. The stress conjugate to logarithmic strain. Int. J. Solids and Structures 12, 1645-56.

Holopainen, S., 2013. Modeling of the mechanical behavior of amorphous glassy polymers under variable loadings and comparison with state-of-the-art model predictions. Mechanics of Materials 66, 35-58.

Holopainen, S., 2014. Influence of damage on inhomogeneous deformation behavior of amorphous glassy polymers. Modeling and algorithmic implementation in a finite element setting. Engng. Fract. Mech. 117, 28-50.

Janssen, R.P.M., Govaert, L.E., Meijer, H.E.H., 2008a. An analytical method to predict fatigue life of thermoplastics in uniaxial loading: Sensitivity to wave type, frequency, and stress amplitude. Macromolecules 41, 2531-40.

Janssen, R.P.M., Kanter, D.K., Govaert, L.E., Meijer, H.E.H., 2008b. Fatigue life predictions for glassy polymers: A constitutive approach. Macromolecules 41, 2520-30.

Kang, Q., 2008. Ratchetting: Recent progresses in phenomenon observation, constitutive modeling and application. International Journal of Fatigue 30, 1448-72.

Kanters, M.J.W., Kurokawa, T., Govaert, L.E., 2016. Competition between plasticity-controlled and crack-growth controlled failure in static and cyclic fatigue of thermoplastic polymer systems. Polymer Testing 50, 101-110.

Krairi, A., Doghri, I., 2014. A thermodynamically-based constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage. Int. J. Plasticity 60, 163-181.

Lawrimore, W.B., Francis, D.K., Bouvard, J.L., Hammi, Y., Horstemeyer, M.F., 2016. A mesomechanics parametric finite element study of damage growth and coalescence in polymers using an Elastoviscoelastic-Viscoplastic internal state variable model. Mechanics of Materials 96, 83-95.
Lemaitre, J., Chaboche, J.L., 1994. Mechanics of Solid Materials. Cambridge University Press, Cambridge.

Lemaitre, J., Desmorat, R., 2005. Engineering Damage Mechanics. Ductile, Creep, Fatigue and Brittle Failures. Springer-Verlag, Berlin.

Lesser, A.J., 2002. Fatigue behavior of polymers. In: Encyclopedia of Polymer Science and Technology 6, 197–251.

Li, H., Chandra, N., 2003. Analysis of crack growth and crack-tip plasticity in ductile materials using cohesive zone models. Int. J. Plasticity 19, 849–882.

Li, X., Hristov, H.A., Yee, A.F., Gidley, D.W., 1995. Influence of cyclic fatigue on the mechanical properties of amorphous polycarbonate. Polymer 36, 759–765.

Liu, W., Gao, Z., Yue, Z., 2008. Steady ratcheting strains accumulation in varying temperature fatigue tests of PMMA. Materials Science and Engineering A 492, 102–109.

Liu, Y., Kang, G., Gao, Q., 2010. A multiaxial stress-based fatigue failure model considering ratchetting-fatigue interaction. International Journal of Fatigue 32, 678–684.

Lu, F., Kang, G., Zhu, Y., Xi, C., Jiang, H., 2016. Experimental observation on multiaxial ratchetting of polycarbonate polymer at room temperature. Polymer Testing 50, 135–144.

Lu, H., Kim, G., 2007. Characteristics of accelerated lifetime behavior of polycarbonate under athermal and high loading frequency conditions. Polymer Testing 26, 839–45.

Lugo, M., Fountain, J.E., Hughes, J.M., Bouvard, J.L., Horstemayer, M.F., 2014. Microstructure-based fatigue modeling of an acrylonitrile butadiene styrene (abs) copolymer. J. Appl. Polym. Sci. 131, 1–12.

Marissen, R., Schudy, D., Kemp, A.V.J.M., Coolen, S.M.H., Duijzings, W.G., der Pol, A.V., Gulick, A.J.V., 2001. The effect of material defects on the fatigue behaviour and the fracture strain of ABS. Journal of Material Sciences 36, 4167–80.

Maxwell, A.S., Broughton, W.R., Dean, G., Sims, G.D., 2005. Review of accelerated ageing methods and lifetime prediction techniques for polymeric materials. Technical Report. National Physical Laboratory Hampton Road, Teddington, Middlesex, TW11 0LW.

Miehe, C., 1998. Comparison of two algorithms for the computation of fourth-order isotropic tensor functions. Comp. Struct. 66, 37–43.

Mortazavian, S., Fatemi, A., 2015. Fatigue behavior and modeling of short fiber reinforced polymer composites: A literature review. International Journal of Fatigue 70, 297–321.

Nittur, P.G., Karlsson, A.M., Carlsson, L.A., 2013. Implementation of a plastically dissipated energy criterion for three dimensional modeling of fatigue crack growth. International Journal of Fatigue 54, 47–55.

Ottosen, N., Stenström, R., Ristinmaa, M., 2008. Continuum approach to high-cycle fatigue modeling. International Journal of Fatigue 30, 996–1006.
doi:doi:10.1016/j.ijfatigue.2007.08.009.

Pasta, S., 2011. Fatigue crack growth through particulate clusters in polycarbonate material. Engng. Fract. Mech. 78, 397–411.

Pijnenburg, K.G.W., Seeling, T., Van der Giessen, E., 2005. Successively refined models for crack tip plasticity in polymer blends. European Journal of Mechanics A/Solids 24, 740–756.
Ravi Chandran, K.S., 2016. Mechanical fatigue of polymers: A new approach to characterize the S-N behavior on the basis of macroscopic crack growth mechanism. Polymer 91, 222–238.

Reese, S., Govinjee, S., 1998. A theory of finite viscoelasticity and numerical aspects. Int. J. Solids and Structures 35, 3455–82.

Seeling, T., Van der Giessen, E., 2015. Computational modeling of rubber-toughening in amorphous thermoplastic polymers: a review. Int. J. Fracture 196, 207–222.
doi:doi:10.1007/s10704-015-0066-6.

Sines, G., 1959. Behavior of metals under complex static and alternating stresses. McGraw-Hill. pp. 145–169.

Del Vecchio, F.J.C., Reis, J.M.L., da Costa Mattos, H.S., 2014. Elasto-viscoplastic behaviour of polyester polymer mortars under monotonic and cyclic compression. Polymer Testing 35, 62–72.

Wallin, M., Holopainen, S., 2012. Modeling of long-term behavior of amorphous glassy polymers. ASME J. Engng. Materials Technol. 135, 1–11.

Wang, Y.Y., Yao, W.X., 2004. Evaluation and comparison of several multiaxial fatigue criteria. International Journal of Fatigue 26, 17 – 25. doi:doi:10.1016/S0142-1123(03)00110-5.

Weber, G., Anand, L., 1990. Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids. Comput. Methods Appl. Mech. Engng. 79, 173–202.

Xi, C., Kang, G., Lu, F., Zhang, J., Jiang, H., 2015. An experimental study on uniaxial ratcheting of polycarbonate polymers with different molecular weights. Materials and Design 67, 644–648.

Xiong, Y., Yu, Q., Jiang, Y., 2014. An experimental study of cyclic plastic deformation of extruded zk60 magnesium alloy under uniaxial loading at room temperature. Int. J. Plasticity 53, 107–124.

Zaïri, F., Naït-Abdelaziz, M., Gloaguen, J.M., Lefebvre, J.M., 2011. A physically-based constitutive model for anisotropic damage in rubber-toughened glassy polymers during finite deformation. Int. J. Plasticity 27, 25–51.

Zhang, J., Jiang, H., Jiang, C., Kang, G., Lu, F., 2015. Accelerated ratcheting testing of polycarbonate using the time-temperature-stress equivalence method. Polymer Testing 44, 8–14.

## Appendix A.

### *Solution of the nonlinear systems*

The proposed algorithm for a solution of the systems (39) and (61) is based on an assumption that those two systems are considered uncoupled, i.e. they are solved in succession. A benefit of the proposed concept is that the integration of the viscoelastic-plastic constitutive model and the fatigue model can be implemented by using suitable numerical schemes, not necessarily one and the same. Concerning e.g. an increasing fatigue damage prior to final failure, an explicit integration of the fatigue model may provide an efficient solution without the dependency on the step size which is a weakness of implicit solution procedures. Since both the models are rate-dependent and long-term periods
are investigated, the algorithmic update is demonstrated by using an implicit (backward Euler) scheme.

When solving the system (39), the state variables are given in the vectorized form, i.e.

$$ Y := [F^e \tilde{W}^p F_1^e s_1]. \qquad (A.1) $$

Using the Newton-Raphson procedure, the solution is $Y^{i+1} = Y^i + \Delta Y$ where the increment is defined as

$$ \Delta Y = -J^{-1}R, \quad \text{and} \quad J := \frac{\partial R}{\partial Y} \qquad (A.2) $$

is the Jacobian. The vectorized form $R := [R_1 \ R_2 \ R_3 \ R_4 \ R_5 \ R_6]$ consists of the residuals introduced in (39).

For a solution of the system (61), fatigue state variables are also expressed by exploiting the vectorized form, i.e.

$$ Y_{\alpha} := [\alpha \ D]. \qquad (A.3) $$

The solution is $Y_{\alpha}^{i+1} = Y_{\alpha}^{i} + \Delta Y_{\alpha}$ where the increment $\Delta Y_{\alpha}$ is calculated in common with (A.2).

### ATS-tensor

Based on the fact that the residuals $R$ in (39) vanish finally at the integration interval for all $F$, i.e. $R(Y(F), F) = 0$, one obtains

$$ \frac{dY}{dF} = - \left( \frac{\partial R}{\partial Y} \right)^{-1} \frac{\partial R}{\partial F} =: -J^{-1} \frac{\partial R}{\partial F}. \qquad (A.4) $$

Using (A.4), one can extract $d\mathbf{v}_1^e/dF$ that is consistent with the applied integration method. In accordance with (A.4), one obtains $dY_{\alpha}/dF$ from (61), and thus all components for the calculation of the ATS-tensor (63) are available. The overall implementation of the scheme is given in Table 3.
<table><caption>Table 3: Solution procedure based on the Kröner-Lee decomposition and the evolution equations based fatigue model.</caption>
  <tbody>
    <tr>
      <td>1. Load data: <b>F</b>, &Delta;&epsilon;, <b>Y</b><sub>n</sub> := [<b>F</b><sub>n</sub><sup>e</sup> <b>W&#x0303;</b><sub>n</sub><sup>p</sup> <b>F</b><sub>1,n</sub><sup>e</sup> s<sub>1,n</sub>], <b>Y</b><sub>&alpha;,n</sub> := [&alpha;<sub>n</sub> <b>D</b><sub>n</sub>].</td>
    </tr>
    <tr>
      <td>2. Initialize: k = 0, <b>F</b><sup>p</sup>|<sub>k=0</sub> = <b>F</b><sub>n</sub><sup>e-1</sup><b>F</b>, <b>F</b><sub>2</sub><sup>e</sup>|<sub>k=0</sub> = <b>F</b><sub>1,n</sub><sup>e-1</sup><b>F</b><sub>n</sub><sup>e</sup>.</td>
    </tr>
    <tr>
      <td>3. Update the internal variables:<br/>(i) Compute &sigma;<sub>n</sub> from (29).<br/>WHILE ||<b>R</b>|| &gt; tol<br/>(ii) Compute <b>D&#x0303;</b><sup>p</sup> by (32).<br/>(iii) Compute s&#x0307;<sub>1</sub> by (33).<br/>(iv) Compute the residuals <b>R</b> according to (39), and the Jacobian J := &part;<b>R</b>/&part;<b>Y</b>.<br/>(v) Update:, <b>Y</b> &larr; <b>Y</b><sub>n</sub> + &Delta;<b>Y</b> by (A.2) and set k &larr; k + 1.<br/>END WHILE LOOP<br/>(vi) Compute &beta;<sub>n</sub> by (35) and (20).<br/>(vii) Compute &beta; by (40) and &Delta;&beta; by (54) taking (62) into account.<br/>IF &beta; &ge; 0 and &Delta;&beta; &ge; 0 THEN fatigue evolves<br/>WHILE ||<b>R</b><sub>&alpha;</sub>|| &gt; tol<br/>(viii) Compute <b>D&#x0303;</b> from (48) and &alpha;<sup>&nabla;</sup> from (42) taking (56), (57), and (60) into account.<br/>(ix) Compute the residuals <b>R</b><sub>&alpha;</sub> according to (61), and the Jacobian J<sub>&alpha;</sub> := &part;<b>R</b><sub>&alpha;</sub>/&part;<b>Y</b><sub>&alpha;</sub>.<br/>(x) Update internal variables, <b>Y</b><sub>&alpha;</sub> &larr; <b>Y</b><sub>&alpha;,n</sub> + &Delta;<b>Y</b><sub>&alpha;</sub> by analogy to (A.2).<br/>END WHILE LOOP<br/>ELSE<br/>(xi) &alpha;<sup>&nabla;</sup> = 0, <b>D&#x0303;</b> = 0.<br/>END IF</td>
    </tr>
    <tr>
      <td>4. Store variables: <b>Y</b> := [<b>F</b><sup>e</sup> <b>W&#x0303;</b><sup>p</sup> <b>F</b><sub>1</sub><sup>e</sup> s<sub>1</sub>], <b>Y</b><sub>&alpha;</sub> := [&alpha; <b>D</b>] and proceed to the equilibrium iteration for <b>F</b> and &Delta;&epsilon;.</td>
    </tr>
  </tbody>
</table>