# A multi-temporal scale model reduction approach for the computation of fatigue damage

Mainak Bhattacharyya, Amélie Fau, Udo Nackenhorst, David Néron, Pierre Ladevèze

## ► To cite this version:

Mainak Bhattacharyya, Amélie Fau, Udo Nackenhorst, David Néron, Pierre Ladevèze. A multi-temporal scale model reduction approach for the computation of fatigue damage. Computer Methods in Applied Mechanics and Engineering, 2018, 340, pp.630-656. 10.1016/j.cma.2018.06.004. hal-01832346

HAL Id: hal-01832346

https://hal.science/hal-01832346

Submitted on 26 Dec 2019

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# A multi-temporal scale model reduction approach for the computation of fatigue damage

Mainak Bhattacharyya<sup>a,b</sup>, Amélie Fau<sup>a</sup>, Udo Nackenhorst<sup>a</sup>, David Néron<sup>b</sup>, Pierre Ladevèze<sup>b</sup>

<sup>a</sup>IBNM, Leibniz Universität Hannover, Appelstraße 9a, 30167 Hannover, Germany

<sup>b</sup>LMT, ENS Paris-Saclay, CNRS, Université Paris Saclay, 61 avenue du Président Wilson, F-94235, Cachan, France

## Abstract

One of the challenges of fatigue simulation using continuum damage mechanics framework over the years has been reduction of numerical cost while maintaining acceptable accuracy. The extremely high numerical expense is due to the temporal part of the quantities of interest which must reflect the state of a structure that is subjected to exorbitant number of load cycles. A novel attempt here is to present a non-incremental LATIN-PGD framework incorporating temporal model order reduction. LATIN-PGD method is based on separation of spatial and temporal parts of the mechanical variables, thereby allowing for separate treatment of the temporal problem. The internal variables, especially damage, although extraneous to the variable separation, must also be treated in a tactical way to reduce numerical expense. A temporal multi-scale approach is proposed that is based on the idea that the quantities of interest show a slow evolution along the cycles and a rapid evolution within the cycles. This assumption boils down to a finite element like discretisation of the temporal domain using a set of “nodal cycles” defined on the slow time scale. Within them, the quantities of interest must satisfy the global admissibility conditions and constitutive relations with respect to the fast time scale. Thereafter, information of the “nodal cycles” can be interpolated to simulate the behaviour on the whole temporal domain. This numerical strategy is tested on different academic examples and leads to an extreme reduction in numerical expense.

## Keywords:

Model reduction, Two-time scale, Fatigue, Damage mechanics, (Visco)plasticity, Cyclic loading, LATIN method

## 1. Introduction

Risk of failure due to fatigue, i.e. due to repetitive loading, is an essential factor of failure for most mechanical components and civil structures (see Sobczyk and Spencer, 1992). Continuum damage mechanics, being a branch of the classical continuum mechanics, introduces an extra internal variable, that quantifies in a thermodynamically consistent framework the loss of load carrying capacity of a structure (see Lemaitre and Desmorat, 2005; Lemaitre, 1996). From the first usage of continuum damage mechanics to predict fatigue (see Chaboche and Lesne, 1988), significant developments have taken place for modelling different types of fatigue processes (see Lemaitre and Desmorat, 2005; Lemaitre, 1996). One of the most important developments, proposed first in Ladevèze and Lemaitre (1984), was the idea of micro-defects closure effects, where the difference in damage behaviour during tension and compression was introduced. This idea has later been used extensively to model fluctuating loads (Desmorat et al., 2007, 2006; Desmorat and Cantournet, 2008; Lemaitre et al., 1999) especially for fatigue applications.

Continuum damage mechanics offers a flexible platform for simulating damage evolution for a large variety of applications, however, this leads to computations, which may be extremely costly. More specifically, the numerical challenge for fatigue simulation in a continuum mechanics framework lies in the number of cycles to be computed, for instance, 10⁴ cycles including 100 time steps per cycle lead to one million of time steps. For fatigue computation applied to a large number of cycles, advanced numerical strategies are then required to circumvent the computational cost due to the large temporal domain.

Email address: mainak.bhattacharyya@ibnm.uni-hannover.de (Mainak Bhattacharyya)
One of the classical approaches in this area is the usage of the so called “jump-in-cycle” technique, where full
blocks of cycles are skipped and only certain cycles are calculated (see Lesne and Savalle, 1989). The detailed
computations of some sets of cycles are used to extrapolate the further evolution of the non-cyclic quantities of
interest, such as damage and accumulated plastic strain. The length of the cycle jumps may be monitored using a
control function to guarantee a certain accuracy (see Cojocaru and Karlsson, 2006). Mostly developed for quasi-linear
systems, some extensions to non-linear behaviour have been proposed using the control function to automatically
shorten or even cancel the cycle jumps (see Nesnas and Saanouni, 2000).

On the other hand, Bhamare et al. (2014) proposes not to skip any cycles but to reduce and optimise the number
of temporal degrees of freedom by using enrichment functions in time for high-cycle fatigue applications.

For combined cycle fatigue in which a combination of large-amplitude, low-frequency and small-amplitude, high-frequency loads is applied to the structure (see Suresh, 2001), homogenisation schemes for time have been proposed (see Guennouni, 1988; Puel and Aubry, 2014). They lie on a prior assumption of a scale separation between a small scale due to the high-frequency load and a large scale due to the low-frequency load. As long as the ratio between both of them allows to consider the two scales as independent, a homogenisation scheme in time can be used similarly to space homogenisation approaches (see Guennouni and Aubry, 1986). This numerical approach, which was initially established for plasticity and quasi-static computations (see Fish et al., 1997), has been more recently extended to compute structures under damaging loads (see Oskay and Fish, 2004) and dynamic effects (see Aubry and Puel, 2010).

For considering effectively large temporal domain, time domain may be decomposed in several subdomains, which
are computed independently in parallel (see Gander, 2015). A comparison of different parallelised strategies has been
proposed in De Vuyst (2016) in the context of thermal applications.

The LATIN method, where the approximation of the solution is looked for at every iteration on a time-space
domain (see e.g. Ladevèze, 1999; Ladevèze, 1989) is used as a solution strategy in this article. Each LATIN iteration
is composed of two steps. On one hand the global equilibrium is solved for the time-space domain of interest as a
linearised problem. The non-linear evolution of the material properties is then evaluated locally also for the whole
time-space domain of interest (see Cognard et al., 1999). This method defined as non-incremental in time is based as
usual on a discretisation scheme in space and in time, but it does not prescribe the accurate evaluation of the problem
at time tᵢ before investigating the subsequent time tᵢ₊₁ as commonly required by incremental schemes. To reduce
the numerical cost, a separation of variables according to the Proper Generalised Decomposition (PGD) method is
employed to handle conveniently the global problem. Using this model order reduction (MOR) technique, the solution
is obtained by solving the Galerkin problem within reduced-order bases, the dimensions of which are much smaller
than the original problem. The bases are generated on the fly by a greedy algorithm. Interested readers are referred to
Chinesta et al. (2011) for a comprehensive review of MOR. The efficiency and robustness of the LATIN-PGD strategy
has been shown for non-linear problems (see Ladevèze, 1999; Ladevèze, 2016) as well as for parametric problems
(see González et al., 2015; Ammar et al., 2015; Aguado et al., 2015) or for quasi-brittle failure (see Vandoren et al.,
2013). Although, PGD-based model reduction in LATIN method provides a drastic decrease in numerical expense, it
is not enough to simulate structures subjected to extremely large number of load cycles. The temporal ODEs in these
cases are very large and sophisticated strategies are necessary to solve them.

Therefore, a multi-temporal scale approach has been proposed to compute efficiently cyclic (visco)plastic prob-
lems in LATIN-PGD framework (see Cognard and Ladevèze, 1993; Ladevèze, 1999). Only a small number of cycles
of interest, called “nodal cycles” are computed. These cycles, which can be interpreted as temporal nodes, delimit
the time elements. Within the time element, the evolution of the quantities of interest is interpolated using shape
functions as for conventional finite element method. This approach results in a two-time scale discretisation with a
rough discretisation between the “nodal cycles”, and a fine discretisation within each “nodal cycle”.

The LATIN-PGD technique has been used to simulate unilateral damage in Bhattacharyya et al. (2017). It has been
proven to be a flexible MOR framework for cyclic damage computations, feasible however, for very few number of
cycles (about 10-20). Hence, the multi-temporal discretisation approach must be used, which may reduce drastically
the numerical expense. The idea of separated temporal scales applied for the calculation of unilateral damage is
introduced in the present study. The Galerkin-based PGD approach used in Cognard and Ladevèze (1993) although
has its benefits in the resolution of the temporal problems, there is no guaranteed convergence of the problem. PGD
based on minimisation of residual, however, provides guaranteed convergence despite being numerically expensive
(Nouy, 2010). A hybrid PGD strategy proposed in Bhattacharyya et al. (2017) benefiting from both is probably a better
option. Therefore the multi-temporal scale methodology used in Cognard and Ladevèze (1993) has to be modified for the usage of hybrid PGD formulation. The use of damage as an internal variable, has its own numerical complexity when incorporated in the multi-temporal LATIN-PGD algorithm, for instance calculation of the initial conditions at each “nodal cycle”.

The structure of the article is as follows. In section 2, the general problem in time-space domain with corresponding admissibility conditions and constitutive formulations for the state and internal variables are presented. Then, the LATIN algorithm incorporating damage and model reduction based on the separation of time and space variables is introduced in section 3. The numerical treatment of the time problem, which allows to reduce the numerical cost, is detailed in section 4, which is the core of this contribution. Finally, several numerical examples are explored in section 5 to analyse the performance of the new algorithm, along with an example depicting the comparison with a reference solution.

## 2. Damage description in continuum mechanics framework

Consider a quasi-static and iso-thermal analysis within time domain $[0, T]$ of a general continuous structure defined by its spatial domain $\Omega$ having a boundary $\partial\Omega$. This structure may be subjected to volumetric forces $f_d$, prescribed displacements $\underline{u}_d$ over a part $\partial_1\Omega$ of the boundary $\partial\Omega$ and to prescribed traction forces $F_d$ over the complementary part $\partial_2\Omega$ of the boundary $\partial\Omega$. The field of internal stresses thus generated is permissible with the applied external forces if the static admissibility condition defined as $\forall \underline{u}^* \in \mathcal{U}^*$

$$ \int_{[0,T]\times\Omega} \sigma : \varepsilon(\underline{u}^*) d\Omega dt = \int_{[0,T]\times\Omega} \underline{f}_d \cdot \underline{u}^* d\Omega dt + \int_{[0,T]\times\partial_2\Omega} \underline{F}_d \cdot \underline{u}^* dS dt, \quad (1) $$

is satisfied. Here $\sigma$ denotes the statically admissible stress tensor, and $\mathcal{U}^*$ is the vector space associated with the space $\mathcal{U}$ of kinematically admissible displacement fields $\underline{u}$. The compatibility of the strain generated within the structure with respect to the prescribed displacement is guaranteed by the satisfaction of the kinematic admissibility condition written as $\forall \sigma^* \in \mathcal{F}^*$

$$ \int_{[0,T]\times\Omega} \sigma^* : \varepsilon \, d\Omega \, dt = \int_{[0,T]\times\partial_1\Omega} \sigma^* \underline{n} \cdot \underline{u}_d \, dS \, dt, \quad (2) $$

where $\mathcal{F}^*$ is the vector space associated with the space $\mathcal{F}$ of statically admissible fields $\sigma$. Under the assumption of infinitesimal strain theory, the kinematically admissible total strain tensor field $\varepsilon$ can be additively decomposed into an elastic part $\varepsilon^e$ and a plastic part $\varepsilon^p$.

A set of constitutive relations characterises the mechanical behaviour of the material, and its specific properties. Here an elasto-(visco)plastic material with kinematic hardening and unilateral damage is considered. From the free energy function, the equations of state for elasto-(visco)plastic materials subjected to unilateral damage are established as (see Lemaitre and Desmorat, 2005)

$$ \varepsilon_{ij}^e = \frac{1 + \nu}{E} \left[ \frac{\langle \sigma_{ij} \rangle}{1 - D} + \frac{\langle -\sigma_{ij} \rangle}{1 - hD} \right] - \frac{\nu}{E} \left[ \frac{\langle \sigma_{kk} \rangle}{1 - D} + \frac{\langle -\sigma_{kk} \rangle}{1 - hD} \right] \delta_{ij}, \quad (3a) $$

$$ \beta_{ij} = C\alpha_{ij}, \quad (3b) $$

$$ Y = \frac{1 + \nu}{2E} \left[ \frac{\langle \sigma_{ij} \rangle \langle \sigma_{ij} \rangle}{(1 - D)^2} + h \frac{\langle -\sigma_{ij} \rangle \langle -\sigma_{ij} \rangle}{(1 - hD)^2} \right] - \frac{\nu}{2E} \left[ \frac{\langle \sigma_{kk} \rangle^2}{(1 - D)^2} + h \frac{\langle -\sigma_{kk} \rangle^2}{(1 - hD)^2} \right], \quad (3c) $$

where the Macauley brackets $\langle\cdot\rangle$ represents the non-negative part of $\bullet$. Equation (3a) is the elastic state law where $E$ is the modulus of elasticity, $\nu$ the Poisson ratio and $D$ the isotropic damage variable. The elastic state equation also shows the effect of the closure of the micro-defects during compression. This phenomenon is governed by the closure parameter $h$ and represents greater reduction of the stiffness in tension as compared to compression. Equation (3b) provides the relation between the back strain tensor $\alpha_{ij}$ and its corresponding thermodynamic force $\beta_{ij}$ through the
material parameter C. The relation between the strain energy release rate Y, which is conjugate to the damage variable D is given by eq. (3c).

For classical standard inelastic materials, the history dependence of the quantities of interest, represented by the evolution equations can be written as (see Bellenger and Bussy, 1998),

$$ \dot{D} = \lambda_d \frac{\partial f_d}{\partial Y}, \qquad (4a) $$

$$ \dot{\varepsilon}_{ij}^{p} = \lambda_{p} \frac{\partial f_{p}}{\partial \sigma_{ij}}, \qquad (4b) $$

$$ \dot{\alpha}_{ij} = \dot{\lambda}_p \frac{\partial f_p}{\partial \beta_{ij}}. \qquad (4c) $$

From the damage function $f_d$ and the damage multiplier $\lambda_d$, the evolution of damage with respect to time is given by eq. (4a). Equation (4b) describes similarly the evolution of the plastic strain with respect to time, where $f_p$ is the yield function and $\lambda_p$ is the plastic multiplier. Finally the evolution of the internal variable for kinematic hardening is defined by eq. (4c).

In this case, the yield function is defined from the Marquis-Chaboche visco-plastic model (see Lemaitre, 1996) as

$$ f_p = \sqrt{\frac{3}{2} \varrho_{ij} \varrho_{ij} + \frac{a}{2C} [\beta_{ij} \beta_{ij}] - \sigma_y}, \qquad (5) $$

where $\varrho_{ij} = \frac{\sigma_{ij}^p}{1-D} - \beta_{ij}$, with $\sigma_{ij}^p$ being the deviatoric part of the stress tensor, $\sigma_y$ being the yield stress, and $a$ is a material parameter for kinematic hardening. Contrary to plasticity, $\lambda_p$ has a closed form solution for visco-plasticity, which is defined by the Norton's law as

$$ \dot{\lambda}_p = k \langle f_p \rangle^n, \qquad (6) $$

with $k$ and $n$ being the viscous coefficient and exponent respectively, which are material dependent (see Lemaitre and Desmorat, 2005).

By analogy with visco-plasticity, the damage multiplier is defined as

$$ \dot{\lambda}_d = k_d \langle f_d \rangle^{n_d}, \qquad (7) $$

where $k_d$ and $n_d$ are material-dependent parameters (see Lemaitre and Desmorat, 2005). This type of damage law is similar to what is used in Bellenger and Bussy (1998). The damage threshold is expressed in terms of a critical energy release rate $Y_0$ through

$$ f_d = Y - Y_0. \qquad (8) $$

Damage increases according to the evolution equation previously defined until reaching the critical damage $D_c$ corresponding to the failure of the material due to the initiation of macro-cracks.

The global admissibility conditions and the local constitutive relations are solved to obtain the quantities of interest using the LATIN-PGD framework.

### 3. LATIN-PGD approach for damage simulation

The numerical scheme is based on the LATIN strategy, which considers the problem on the whole time-space domain or a subset of the time-space domain at every iteration (see Ladevèze, 1999). It does not compute consequently every time increment after reaching convergence for the previous instant, but the solution is improved iteratively for all space-time points of the considered domain. The LATIN iterations to solve the coupled damage elasto-viscoplastic problem lie on the successive computations of the global admissibilities on one hand, and of the non-linear evolution equations and non-linear elastic law on another hand. The desired quantities of interest that indicate the state of a structure is given by the solution set $s = \{\dot{\varepsilon}^p, \varepsilon^e, \dot{\alpha}, \dot{D}, \sigma, \beta, Y\}$.
LATIN numerical scheme begins by the resolution of the problem in the time-space domain of interest considering the material remains in an undamaged elastic zone for the exact boundary conditions of the problem. Thereafter, at each LATIN iteration, plastic and damage corrections are successively added to the linear elastic solution. Each LATIN iteration tackles two successive steps. The first step deals with the local non-linear problem. It is solved in the space $\Gamma$ defined from the evolution equations and the elastic state law, which is not linearisable due to the presence of unilateral damage. The subsequent step is to build the solution set based on linear equations which might be global. The problem is solved in the space $A_d$ which is defined from the admissibility conditions and the state laws for damage and kinematic hardening. The strain energy release rate is updated from the knowledge of the stress tensor and the damage variable at the end of each iteration as a post-process.

The transfer of information between the solution sets of the two manifolds takes place through certain linear operators, called search direction operators. This two-step algorithm seeking the approximate solution at the global stage $s_i \in A_d$ and the approximation at the local stage $\hat{s}_{i+1/2} \in \Gamma$ continues till convergence is reached. The convergence, theoretically, is obtained when the exact solution $s_{ex} \in A_d \cap \Gamma$ can be calculated. However, in a numerical point of view, the convergence criterion is based on the distance between $s_{i+1}$ and $\hat{s}_{i+1/2}$, which should be lower than a pre-defined tolerance to achieve convergence.

### 3.1. Local stage

The evolution equations as well as the elastic state law, which are local in space and non-linear, are solved during the local stage. From the solution set $s_i \in A_d$, the estimation of $\hat{s}_{i+1/2} \in \Gamma$ satisfies the search direction equation

$$ \left[ \begin{array}{c} \hat{D}_{i+1/2} - \dot{D}_i \\ \hat{\boldsymbol{\varepsilon}}_{i+1/2}^p - \dot{\boldsymbol{\varepsilon}}_i^p \\ - (\hat{\boldsymbol{\alpha}}_{i+1/2} - \dot{\boldsymbol{\alpha}}_i) \\ \hat{\boldsymbol{\varepsilon}}_{i+1/2}^e - \boldsymbol{\varepsilon}_i^e \end{array} \right] + \mathbf{B}^+ \left[ \begin{array}{c} \hat{Y}_{i+1/2} - Y_i \\ \hat{\boldsymbol{\sigma}}_{i+1/2} - \boldsymbol{\sigma}_i \\ \hat{\boldsymbol{\beta}}_{i+1/2} - \boldsymbol{\beta}_i \\ \hat{\boldsymbol{\sigma}}_{i+1/2} - \boldsymbol{\sigma}_i \end{array} \right] = 0, \quad (9) $$

where the direction of ascent $\mathbf{B}^+$ is assumed to be vertical, i.e. $(\mathbf{B}^+)^{-1} = 0$ (see Ladevèze, 1999). The approximation $\hat{s}_{i+1/2}$ is defined from the search direction equation along with the non-linear elastic law and the evolution equations. Once the solution set $\hat{s}_{i+1/2}$ is calculated, it is used thereafter to solve the global stage to obtain $s_{i+1}$.

### 3.2. Global stage including PGD

The solution set $s_{i+1} \in A_d$ is defined from the admissibility conditions, the state laws for damage and kinematic hardening, and the search direction equations, which are given as

$$ \left[ \begin{array}{c} \dot{D}_{i+1} - \hat{D}_{i+1/2} \\ \hat{\boldsymbol{\varepsilon}}_{i+1}^p - \hat{\boldsymbol{\varepsilon}}_{i+1/2}^p \\ - (\dot{\boldsymbol{\alpha}}_{i+1} - \hat{\boldsymbol{\alpha}}_{i+1/2}) \\ \boldsymbol{\varepsilon}_{i+1}^e - \hat{\boldsymbol{\varepsilon}}_{i+1/2}^e \end{array} \right] - \mathbf{B}^{-} \left[ \begin{array}{c} Y_i - \hat{Y}_{i+1/2} \\ \boldsymbol{\sigma}_{i+1} - \hat{\boldsymbol{\sigma}}_{i+1/2} \\ \boldsymbol{\beta}_{i+1} - \hat{\boldsymbol{\beta}}_{i+1/2} \\ \boldsymbol{\sigma}_{i+1} - \hat{\boldsymbol{\sigma}}_{i+1/2} \end{array} \right] = 0 \quad (10) $$

with $\mathbf{B}^{-} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & \mathbf{H}^{-} & 0 \\ 0 & 0 & \mathbf{C}^{-1} \end{bmatrix}$, which means $\dot{D}_{i+1} = \hat{D}_{i+1/2}$. The operator $\mathbf{H}^{-}$ lies in the tangent space to the solution set $\hat{s}_{i+1/2}$ in the manifold $\Gamma$, $\mathbf{C}$ denotes the Hooke tensor for the undamaged material. It should be noted that $\mathbf{H}^{-}$ is of the form $\mathbf{H}^{-} = \begin{bmatrix} \mathbf{H}_{\sigma} & \mathbf{H}_{\sigma\beta} \\ \mathbf{H}_{\beta\sigma} & \mathbf{H}_{\beta} \end{bmatrix}$. For the sake of simplicity, the off-diagonal terms are set to zero (see Bhattacharyya

et al., 2017; Relun et al., 2011), which results into $\mathbf{B}^{-} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & \mathbf{H}_{\sigma} & 0 & 0 \\ 0 & 0 & \mathbf{H}_{\beta} & 0 \\ 0 & 0 & 0 & \mathbf{C}^{-1} \end{bmatrix}$.

The search direction equation for hardening variables along with the state equation for kinematic hardening can be expressed as

$$ -(\dot{\boldsymbol{\alpha}}_{i+1} - \hat{\boldsymbol{\alpha}}_{i+1/2}) = \mathbf{H}_{\beta}(\boldsymbol{\Lambda}\boldsymbol{\alpha}_{i+1} - \hat{\boldsymbol{\beta}}_{i+1/2}), \quad (11) $$
where $\Lambda$ is a tensor similar to the Hooke tensor which contains the material parameter $C$. Once the search direction is known, the hardening variables can be computed from eq. (10).

The quantities of interest are represented in a corrective form at iteration $i+1$ as $\Delta\sigma_{i+1} = \sigma_{i+1} - \sigma_i$, $\Delta\varepsilon_{i+1}^e = \varepsilon_{i+1}^e - \varepsilon_i^e$ and $\Delta\dot{\varepsilon}_{i+1}^p = \dot{\varepsilon}_{i+1}^p - \dot{\varepsilon}_i^p$. To begin with the developments, consider the search direction equation relating elastic strain and stress written in terms of corrections,

$$ \Delta\sigma_{i+1} = C\Delta\varepsilon_{i+1}^e - \Delta\mathbf{R}_{i+1}, \qquad (12) $$

with the residual stress term $\Delta\mathbf{R}_{i+1}$ is given by

$$ \Delta\mathbf{R}_{i+1} = (\sigma_i - \hat{\sigma}_{i+1/2}) - C(\varepsilon_i^e - \hat{\varepsilon}_{i+1/2}^e). \qquad (13) $$

Introducing a residual strain term $\Delta\varepsilon_{i+1}^R$ given by

$$ \Delta\varepsilon_{i+1}^R = C^{-1}\Delta\mathbf{R}_{i+1}, \qquad (14) $$

eq. (12) can be represented as,

$$ \Delta\sigma_{i+1} = C(\Delta\varepsilon_{i+1}^e - \Delta\varepsilon_{i+1}^R). \qquad (15) $$

Using the additive strain decomposition relationship (infinitesimal strain theory), i.e.

$$ \Delta\varepsilon_{i+1}^e = \Delta\varepsilon_{i+1} - \Delta\varepsilon_{i+1}^p, \qquad (16) $$

eq. (15) is transformed into

$$ \Delta\sigma_{i+1} = C(\Delta\varepsilon_{i+1} - \Delta\varepsilon_{i+1}^p) - C\Delta\varepsilon_{i+1}^R. \qquad (17) $$

The idea herein is to represent the plastic strain in PGD form such that the static admissibility condition is satisfied. If there is no damage corrections (for instance visco-plasticity), the residual strain term is zero and the satisfaction of the static admissibility is straightforward, and the total strain, plastic strain, and stress can be represented as separated variable form (see Relun et al., 2011). However, for non-zero residual strain, it is difficult to have a separable representation for all the aforementioned quantities and satisfying the static admissibility condition. Therefore an additive split of the total strain is proposed, such that

$$ \Delta\varepsilon_{i+1} = \Delta\varepsilon'_{i+1} + \Delta\tilde{\varepsilon}_{i+1}, \qquad (18) $$

where both the quantities have to be kinematically admissible to zero. Plugging this into eq. (17)

$$ \Delta\sigma_{i+1} = C(\Delta\varepsilon'_{i+1} - \Delta\varepsilon_{i+1}^p) + C(\Delta\tilde{\varepsilon}_{i+1} - \Delta\varepsilon_{i+1}^R), \qquad (19) $$

such that

$$ \Delta\sigma'_{i+1} = C(\Delta\varepsilon'_{i+1} - \Delta\varepsilon_{i+1}^p), \quad \Delta\tilde{\sigma}_{i+1} = C(\Delta\tilde{\varepsilon}_{i+1} - \Delta\varepsilon_{i+1}^R), \qquad (20) $$

with the complete description of stress is given by,

$$ \Delta\sigma_{i+1} = \Delta\sigma'_{i+1} + \Delta\tilde{\sigma}_{i+1}. \qquad (21) $$

The quantities $\Delta\sigma'_{i+1}$, $\Delta\varepsilon'_{i+1}$, $\Delta\varepsilon_{i+1}^p$ will satisfy the static admissibility through PGD representation, and the quantities $\Delta\tilde{\sigma}_{i+1}$, $\Delta\tilde{\varepsilon}_{i+1}$, will satisfy the static admissibility bereft of PGD as $\Delta\varepsilon_{i+1}^R$ is a known quantity.

If the plastic strain is obtained through the integration of eq. (4b) and not as rate form, the additive strain decomposition relationship can directly be used in local stage to obtain the total strain. The global stage in that case will involve solving the static admissibility in terms of the total strain with the right-hand side of the equation already obtained from the local stage. This will give a direct PGD representation of the total strain, and plastic strain will not come into the scenario of the global stage. To be more general, if all the constitutive relations i.e. eqs. (3) and (4) are solved in the local stage, which involves temporal integration of eq. (4), the total strain is the only quantity that needs to be calculated in the global stage, and the stress can be obtained from a search direction equation relating stress and total strain. In that case, no additive split as mentioned before is necessary (see Allix et al., 1989, in case of damage).
The problem in that case however is the temporal integration at each Gauss point for the internal variables, which can render the local stage extremely expensive, hence, the quantities namely the plastic strain are obtained in rate form, and integration is performed only in the global stage through PGD-based model reduction. For viscoplasticity, this approach is straightforward, as the Hooke's tensor is constant (see Relun et al., 2011), however for unilateral damage, especially for non-proportional loading, the aforementioned additive split becomes a necessity. The idea behind this whole ordeal is to have the global stage linear for PGD based variable separation.

Now, the search direction relating to the plastic strain is written through corrective terms as

$$ \Delta \dot{\boldsymbol{\epsilon}}_{i+1}^{p} - \mathbf{H}_{\sigma} \Delta \boldsymbol{\sigma}_{i+1} + \bar{\boldsymbol{\Lambda}}_{i+1} = 0, \qquad (22) $$

where $\bar{\boldsymbol{\Lambda}}_{i+1}$ provides information from the local stage and is given as

$$ \bar{\boldsymbol{\Lambda}}_{i+1} = \mathbf{H}_{\sigma}(\hat{\boldsymbol{\sigma}}_{i+1/2} - \boldsymbol{\sigma}_{i}) - (\hat{\boldsymbol{\epsilon}}_{i+1/2}^{p} - \dot{\boldsymbol{\epsilon}}_{i}^{p}). \qquad (23) $$

To be consistent with the additive split, introduced before, eq. (22) should be expressed in terms of those quantities where separated variable representations will be used. Following this logic, eq. (22) is transformed into

$$ \Delta \dot{\boldsymbol{\epsilon}}_{i+1}^{p} - \mathbf{H}_{\sigma} \Delta \boldsymbol{\sigma}'_{i+1} + \bar{\boldsymbol{\Lambda}}_{i+1} \approx 0. \qquad (24) $$

The mathematical validity of this equation is acceptable in the presence of low damage value, which is indeed the case for metals (see Lemaitre and Desmorat, 2005, for instance), investigated in this article. The search direction operators $\mathbf{H}_{\sigma}$ and $\mathbf{H}_{\beta}$ are given in Bhattacharyya et al. (2017) as,

$$ \mathbf{H}_{\sigma} = kn \langle f^p \rangle^{n-1} \frac{1}{(1-D)^2} \frac{3/2 \boldsymbol{\rho}}{\sqrt{J_2}} \otimes \frac{3/2 \boldsymbol{\rho}}{\sqrt{J_2}} + k \langle f^p \rangle^n \frac{1}{(1-D)^2} \frac{3/2 \left( \sqrt{J_2} \mathbb{I} - 3/2 \frac{\boldsymbol{\rho} \otimes \boldsymbol{\rho}}{\sqrt{J_2}} \right)}{J_2} \qquad (25a) $$

$$ \mathbf{H}_{\beta} = kn \langle f^p \rangle^{n-1} \left( -\frac{3/2\boldsymbol{\rho}}{\sqrt{J_2}} + \frac{a}{C}\mathbb{I}\boldsymbol{\beta} \right) \otimes \left( -\frac{3/2\boldsymbol{\rho}}{\sqrt{J_2}} + \frac{a}{C}\mathbb{I}\boldsymbol{\beta} \right) + k \langle f^p \rangle^n \frac{3/2 \left( \sqrt{J_2} \mathbb{I} - 3/2 \frac{\boldsymbol{\rho} \otimes \boldsymbol{\rho}}{\sqrt{J_2}} \right)}{J_2} \frac{a}{C} \mathbb{I} \qquad (25b) $$

with $J_2 = \frac{3}{2}\boldsymbol{\rho} \cdot \boldsymbol{\rho}$ and $\mathbb{I}$ being the identity matrix.

PGD is then used to estimate the correction terms corresponding to the plastic behaviour $\Delta \boldsymbol{\sigma}'_{i+1}$ and $\Delta \dot{\boldsymbol{\epsilon}}_{i+1}^p$ in a separable form.

### 3.2.1. Separation of variables

As the LATIN global stage is the most computationally demanding part of the algorithm, it commonly incorporates a model reduction method, namely PGD (see Ladevèze, 1999, 2014). Considering any separable quantity of interest defined over the complete space-time domain represented by $\chi(\underline{x}, t)$, it is possible to approximate it as a finite sum of products of two separate functions, one depending on time and the other on space. If $\iota$ is the number of terms chosen to approximate the quantity of interest, $\chi$ can be written as

$$ \chi(\underline{x}, t) \approx \sum_{j=1}^{\iota} \chi_j^{\iota}(\underline{x}) \chi_j^{\iota}(t). \qquad (26) $$

This separation allows to rewrite the governing equation into separated space and time problems, which simplifies drastically the numerical complexity of the problem (see Chinesta et al., 2010b). For the current problem of interest, $\Delta \boldsymbol{\sigma}'$ and $\Delta \dot{\boldsymbol{\epsilon}}^p$ are written in separable variable forms and at each LATIN iteration the objective is to obtain the approximation either by updating the temporal basis or if needed by enriching the spatial and temporal bases by adding space-time PGD pairs.
### 3.2.2. Updating stage

Considering $\mu$ space-time modes have already been generated to estimate the stress and plastic strain rate at the previous iteration $i$, the decomposition at iteration $i + 1$ is written by re-using the space functions and re-computing the time functions (see Bhattacharyya et al., 2017, for details). The separation of the quantities of interest is of the form

$$ \begin{aligned} \Delta \dot{\boldsymbol{\varepsilon}}_{i+1}^{p}(\underline{x}, t) &= \sum_{j=1}^{\mu} \Delta \dot{\lambda}_{j}(t) \bar{\varepsilon}_{j}^{p}(\underline{x}), \\ \Delta \boldsymbol{\sigma}'_{i+1}(\underline{x}, t) &= \sum_{j=1}^{\mu} \Delta \lambda_{j}(t) \bar{\sigma}_{j}(\underline{x}), \end{aligned} \qquad (27) $$

with $\{\Delta\lambda_j\}_{j=1}^\mu$ being corrections of the time functions $\{\lambda_j\}_{j=1}^\mu$. If this updating stage does not improve the accuracy of the approximation significantly, the greedy algorithm allows to enlarge the space and time bases.

### 3.2.3. Enrichment of space-time bases

If the approximation obtained from the update phase is not adequate enough, the reduced order basis is enriched using the PGD-based approximation. The numerical strategy chosen is to add one new space-time pair at a given LATIN iteration. Although, the method, in principle, is not restricted to the addition of single space-time mode per LATIN iteration, it appears as an optimal strategy to avoid expansion of the basis at the initial iterations whereas the information from the local stage is quite distant from accuracy, for highly non-linear problems. In case of enrichment of the bases, the correction terms for the plastic strain and stress at LATIN iteration $i + 1$ are respectively

$$ \begin{aligned} \Delta \dot{\boldsymbol{\varepsilon}}_{i+1}^{p}(\underline{x}, t) &= \dot{\lambda}_{\mu+1}(t) \bar{\varepsilon}_{\mu+1}^{p}(\underline{x}), \\ \Delta \boldsymbol{\sigma}'_{i+1}(\underline{x}, t) &= \lambda_{\mu+1}(t) \bar{\sigma}_{\mu+1}(\underline{x}), \end{aligned} \qquad (28) $$

where $\dot{\lambda}_{\mu+1}$, $\bar{\varepsilon}_{\mu+1}^p$, and $\bar{\sigma}_{\mu+1}$ are newly calculated. Contrary to the previous versions of LATIN-PGD techniques, a hybrid strategy is used to calculate the new PGD-pair. It has been tested for different academic applications and has shown robustness (see Bhattacharyya et al., 2017). However, the gain in numerical cost and memory comparing to a full computation is not large enough to enable the computation of large number of cycles that a structure can sustain before reaching the critical damage. In the following section, an extension of this method is proposed to reduce the numerical cost of the time problem to handle fatigue damage loads while maintaining low CPU cost and reasonable computational time.

## 4. A two-time discretisation approach

By decoupling the time and space dependencies, the LATIN-PGD strategy reduces drastically the size of the linear system, and thus allows to tackle nowadays sophisticated industrial problems (see Ladevèze, 2016). The spatial dependency is usually the most expensive part of the mechanical problem, and its numerical cost can be efficiently handled by using domain decomposition strategies (see Ladevèze and Lorong, 1991) or multi-scale approaches (see Chinesta et al., 2010a; Néron and Ladevèze, 2010). Considering fatigue, the number of cycles may be very large, and so the cost of the time-dependent problem may become prohibitive for the computation. For instance, if $10^5$ cycles are simulated with 100 time steps per cycle, it will result in $10^7$ total time steps. The fundamental goal of this contribution is to propose a numerical strategy to tackle this difficulty and thereby to allow the estimation of the fatigue lifetime of structures subjected to cyclic loading using continuum damage mechanics.

### 4.1. Multi-time scale discretisation

Consider structures under periodic loads with constant time periods, but with amplitudes that may vary with time. When the temporal domain $[0, T]$ includes several thousands of cycles, a multi-temporal discretisation scheme has to be used. Originating from the propositions of Cognard and Ladevèze (1993) for elasto-viscoplasticity, this idea is based on the slow evolution of the quantities of interests along the cycles and their fast evolution within a cycle. Based
on these assumptions, a coarse temporal discretisation $\tilde{\theta}$ is defined on the domain $[0, T]$ representing the beginning of every cycle, and a fine temporal discretisation $\tau$ defined within a cycle. Any instant of interest $t$ can be discretised by the pair $(\tilde{\theta}, \tau)$ as

$$t = \tilde{\theta}_i + \tau_i, \text{ with } i \in [0, N-1] \text{ and } \tau_i \in [0, \Delta T], \forall t \in [0, T], \quad (29)$$

where $N$ is the number of cycles included in the time domain.

## 4.2. Finite element like time interpolation scheme

The two-time scale allows to define certain cycles where the mechanical problem is solved using the LATIN-PGD technique, which are called the “nodal cycles” (see fig. 1). The temporal domain is divided into “time elements”, which are delimited by the “nodal cycles”. Once the “nodal cycles” are computed, one-dimensional interpolations are used to estimate the quantities of interest over the “time elements” (see Cognard and Ladevèze, 1993). For a particular “time element” $[\Theta_m, \Theta_{m+1}]$ containing $\tilde{p}$ cycles and demarcated by the “nodal cycles” $m$ and $m+1$, a local variable $\theta_k, \forall k \in [0, \tilde{p}-1]$ is introduced which represents the coarse discretisation for the given “time element”. This representation provides the time $t$ to be defined as

$$t = \theta_k + \tau_k, \text{ with } k \in [0, \tilde{p}-1] \text{ and } \tau_k \in [0, \Delta T], \forall t \in [\Theta_m, \Theta_{m+1} + \Delta T]. \quad (30)$$

For the given “time element” $[\Theta_m, \Theta_{m+1}]$ temporal shape functions $v_m$ and $v_{m+1}$ are defined as

$$v_m(\theta_k) = \frac{\Theta_{m+1} - \theta_k}{\Theta_{m+1} - \Theta_m}, \quad v_{m+1}(\theta_k) = \frac{\theta_k - \Theta_m}{\Theta_{m+1} - \Theta_m}, \quad \forall t \in [\Theta_m, \Theta_{m+1} + \Delta T]. \quad (31)$$

<figure><img src="image_2.png" /><figcaption>Figure 1: Discretisation of the time domain using two time scale discretisations: $\theta$ along the cycles and $\tau$ within a cycle, $\tilde{p}$ being the number of cycles for the time element $[\Theta_m, \Theta_{m+1}]$</figcaption></figure>

If any temporal quantity of interest $\chi$ is known at the “nodal cycles” $m$ and $m+1$, the approximation on the whole “time element” $[\Theta_m, \Theta_{m+1}]$, containing $\tilde{p}$ cycles is given by

$$\chi(t = \theta_k + \tau_k) = v_m(\theta_k)\chi(\Theta_m + \tau_m) + v_{m+1}(\theta_k)\chi(\Theta_{m+1} + \tau_{m+1}), \forall t \in [\Theta_m, \Theta_{m+1} + \Delta T], \quad (32)$$

with $\tau_m = \tau_{m+1} = \tau_k$ belonging to $[0, \Delta T]$ as the time period is constant and $k = \{0, 1, \dots, \tilde{p}-1\}$. The discretisation is defined such that $\theta_k = \Theta_m + k\Delta T$, $\theta_0 = \Theta_m$ and $\theta_{\tilde{p}-1} = \Theta_{m+1}$.

It can be noticed that compared to classical finite element interpolations, the interpolations are based not on certain points (nodes) but on whole sets of points (“nodal cycles”). The finite element like description in time for the quantities of interest lies on a set of time elements which spans the temporal domain and an *a posteriori* time interpolation scheme is used with the knowledge of the “nodal cycles”.
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic representation of the multi-temporal scale LATIN-PGD algorithm</figcaption></figure>

After computing the first few cycles using the LATIN-PGD approach as presented in Bhattacharyya et al. (2017), the two-time scheme is used. This calculation of the initial cycles can be interpreted as “training stage”, and the reduced order basis (ROB) obtained from it is reused for the subsequent “nodal cycles”. For each “nodal cycle”, the two-step LATIN-PGD algorithm is used till a converged solution is found. In details, the first element is built from the “nodal cycle” 0, which is the last cycle of the “training stage”, defined over [Θ₀, Θ₀ + ΔT]. Then, from “nodal cycle” 0, the goal is to calculate “nodal cycle” 1 defined over [Θ₁, Θ₁ + ΔT], and successively from knowing “nodal cycle” 1, computing “nodal cycle” 2 and so on (see fig. 2). The computation is pursued for the different “nodal cycles” in succession until the ultimate “nodal cycle” concluding the time domain of interest is calculated, or until the critical damage has been reached.

### 4.3. Computation of one “nodal cycle”

The two-step LATIN-PGD algorithm is used independently at a particular “nodal cycle” till a converged solution is obtained, and only after that the next “nodal cycle” is computed. After initialisation, the local stage and global stage calculations work as explained in section 3.

#### 4.3.1. Initialisation

Unlike classical LATIN strategy (see Ladevèze, 1999), the initialisation for a “nodal cycle” is not based on linear elasticity but benefits from the knowledge of the previous “nodal cycle”. The quantities of interest can be divided into two groups: cyclic quantities, and non-decreasing quantities. The stresses, elastic strains, variables for kinematic hardening, and the time functions $\{\lambda_j\}_{j=1}^\mu$ belong to the former group, and quantities like damage belong to the latter.

The initialisation of the cyclic quantities at “nodal cycle” $m + 1$ is based on the duplication of their estimation at “nodal cycle” $m$. Now, to maintain continuity with respect to the successive “nodal cycles”, it is necessary that the final condition of the previous nodal cycle becomes the initial condition of the current nodal cycle. This idea is similar to the initial conditions forced in jump cycle methods (see Lemaitre and Desmorat, 2005) and is meant to avoid any non-physical discontinuity. Therefore the duplicated quantity must be transformed to take this logic into account.

Consider any cyclic quantity of interest $\vartheta(\Theta_m + \tau_m)$ defined $\forall \tau_m \in [0, \Delta T]$ and at $\Theta_m$. This quantity should be transformed into $\tilde{\vartheta}(\Theta_m + \tau_m)$, such that

$$ \tilde{\vartheta}(\Theta_m) = \tilde{\vartheta}(\Theta_m + \Delta T) = \vartheta(\Theta_m + \Delta T). \quad (33) $$

A linear transformation rule can be assumed which is of the form

$$ \tilde{\vartheta}(\Theta_m + \tau_m) = \vartheta(\Theta_m + \tau_m) + a_{\vartheta} \cdot (\Theta_m + \tau_m) + b_{\vartheta}, \quad (34) $$

where $a_{\vartheta}$ and $b_{\vartheta}$ are parameters obtained from the periodicity condition of eq. (33). The initialisation $\vartheta^0(\Theta_{m+1} + \tau_{m+1})$ for the “nodal cycle” $m + 1$ is thus based on the transformed quantity computed for the “nodal cycle” $m$ as

$$ \vartheta^0(\Theta_{m+1} + \tau_{m+1}) = \tilde{\vartheta}(\Theta_m + \tau_m) = \vartheta(\Theta_m + \tau_m) + \frac{\vartheta(\Theta_m) - \vartheta(\Theta_m + \Delta T)}{\Delta T} \cdot \tau_m + \vartheta(\Theta_m + \Delta T) - \vartheta(\Theta_m). \quad (35) $$
Damage, on the other hand, being non-cyclic and non-decreasing function of time, is initialised as constant over the whole “nodal cycle” $m + 1$, such that the constant magnitude is the same as that obtained at $\Theta_m + \Delta T$. It should be noted that, as “nodal cycle” 0 is the last cycle of the “training stage”, it contains the boundary condition and elastic part of the problem through the initialisation of the “training stage”. Thereafter when this “nodal cycle” is duplicated to “nodal cycle” 1, the boundary condition and the elastic contributions are already taken into account, and holds true for the initialisation of all the subsequent “nodal cycles”. Of course, this is true only for fixed boundary conditions with constant amplitude loading. For moving loads, or for variable amplitudes, separate elastic initialisation has to be calculated for each “nodal cycle”, in this initialisation phase.

Thereafter, the quantities of interest are estimated iteratively using the two-step LATIN algorithm over “nodal cycle” $m + 1$.

### 4.3.2. Local stage

The integration of eq. (3a) to obtain the damage variable at each nodal cycle is the main difficulty in the local stage. To integrate eq. (3a) over “nodal cycle” $m + 1$ the initial condition, i.e. $D(\Theta_{m+1})$ must be known. A different approach compared to the jump cycle method, where the initial condition is obtained through linear extrapolation, is proposed (see Bhattacharyya et al., 2018).

If the damage increments for all $\bar{p} - 1$ cycles are known (see fig. 3), $D(\Theta_m = \theta_0)$ can be used to calculate $D(\Theta_{m+1} = \theta_{\bar{p}-1})$ as

$$D(\theta_{\bar{p}-1}) = D(\theta_0) + \Delta D_{ini} + \sum_{k=1}^{\bar{p}-2} \Delta D_k, \qquad (36)$$

where $\Delta Ds$ are the damage increment at each cycle as shown in fig. 3. However, calculating all the damage increments individually can be expensive so the idea is to calculate only the increments at “nodal cycle” $m$ and $m + 1$, i.e. $\Delta D_{ini}$ and $\Delta D_{fin}$, and then to interpolate linearly the increments for the intermediate cycles from them using eq. (31). Using this assumption, the initial value of damage for the “nodal cycle” $m + 1$ can be written as

$$D(\Theta_{m+1}) = D(\Theta_m) + \sum_{k=0}^{\bar{p}-2} \left( \frac{\Theta_{m+1} - \theta_k}{\Theta_{m+1} - \Theta_m} \Delta D_{ini} + \frac{\theta_k - \Theta_m}{\Theta_{m+1} - \Theta_m} \Delta D_{fin} \right), \qquad (37)$$

where $\Delta D_{ini}$ and $\Delta D_{fin}$ are the damage increments for “nodal cycles” $m$ and $m + 1$ respectively. These increments can easily be calculated by solving the evolution eq. (4a) in the corresponding temporal domains using zero initial conditions.

The other evolution equations, i.e. eqs. (3b) and (3c) do not require any integration and plastic strain and back strain are obtained in rate form and integrations are performed only in the global stage.

### 4.3.3. Global stage

The boundary conditions of the structural problem are already taken into account in the initialisation part of a given “nodal cycle”. The quantities of interest are therefore obtained in terms of corrections such that they are kinematically admissible to zero. Following the idea introduced in section 3.2, the static admissibility conditions for a given nodal cycle $m + 1$ and for a given LATIN iteration $i + 1$ can be written $\forall \underline{u}^* \in \mathcal{U}_0$ as

$$\int_{[\Theta_{m+1}, \Theta_{m+1}+\Delta T]\times\Omega} \Delta\sigma'_{i+1} : \varepsilon(\underline{u}^*) d\Omega dt = 0, \qquad (38a)$$

$$\int_{[\Theta_{m+1}, \Theta_{m+1}+\Delta T]\times\Omega} \Delta\tilde{\sigma}_{i+1} : \varepsilon(\underline{u}^*) d\Omega dt = 0, \qquad (38b)$$

with

$$\Delta\sigma'_{i+1} = C(\Delta\varepsilon'_{i+1} - \Delta\varepsilon^p_{i+1}), \qquad (39a)$$

$$\Delta\tilde{\sigma}_{i+1} = C(\Delta\tilde{\varepsilon}_{i+1} - \Delta\tilde{\varepsilon}^R_{i+1}). \qquad (39b)$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Calculating $D(\Theta_{m+1})$, knowing $D(\Theta_m)$</figcaption></figure>

The correction to the plastic strain rate and the corresponding part of the total strain as well as the displacement can be represented through separated variable formulation as

$$ \Delta \boldsymbol{\varepsilon}_{i+1}^{p} = \lambda^{p}(t) \bar{\varepsilon}^{p}(\underline{x}) \Rightarrow \Delta \dot{\boldsymbol{\varepsilon}}_{i+1}^{p} = \dot{\lambda}^{p}(t) \bar{\varepsilon}^{p}(\underline{x}), \quad (40a) $$

$$ \Delta \underline{u}_{i+1} = \lambda^u(t) \bar{u}(\underline{x}) \Rightarrow \Delta \varepsilon'_{i+1} = \lambda^u(t) \bar{\varepsilon}(\underline{x}). \quad (40b) $$

The test field in eq. (38a) when represented in a separated form, i.e. $\underline{u}^* = \lambda^{u*}\bar{u} + \lambda^u\bar{u}^*$ with $\bar{u}^* \in \mathcal{U}_0$ and $\lambda^{u*}$ without any specific condition, the problem (eq. (38a) along with eq. (39a)) can be separated into sovereign spatial and temporal problems. The temporal problem aims at finding $\lambda^u$ so that $\forall \lambda^{u*}$,

$$ \int_{[\Theta_{m+1}, \Theta_{m+1} + \Delta T]} \lambda^u \lambda^{u*} dt \int_{\Omega} \mathbf{C} \bar{\varepsilon} : \varepsilon(\bar{u}) d\Omega = \int_{[\Theta_{m+1}, \Theta_{m+1} + \Delta T]} \lambda^p \lambda^{u*} dt \int_{\Omega} \mathbf{C} \bar{\varepsilon}^p : \varepsilon(\bar{u}) d\Omega. \quad (41) $$

The simplistic assumption from eq. (41) is $\lambda^u = \lambda^p = \lambda$ (see Relun et al., 2011, for instance), which results in the corrective terms of stress $\Delta\sigma'_{i+1}$ be written as

$$ \Delta\sigma'_{i+1} = \lambda(t) \bar{\sigma}(\underline{x}). \quad (42) $$

The spatial problem thereafter attempts to find $\bar{u} \in \mathcal{U}_0$ such that

$$ \int_{\Omega} \mathbf{C} \bar{\varepsilon} : \varepsilon(\bar{u}^*) d\Omega = \int_{\Omega} \mathbf{C} \bar{\varepsilon}^p : \varepsilon(\bar{u}^*) d\Omega, \quad \forall \bar{u}^* \in \mathcal{U}_0. \quad (43) $$

This allows to define operators $\mathbb{E}$ and $\mathbb{C}$ (see Appendix A for details) in the finite element space, that links $\bar{\varepsilon}^p$ with the other space functions as

$$ \bar{\varepsilon} = \mathbb{E}\bar{\varepsilon}^p, \quad (44a) $$

$$ \bar{\sigma} = \mathbb{C}\bar{\varepsilon}^p \quad (44b) $$

The quantities of interest depending on the plastic strain thereby become

$$ \Delta \dot{\varepsilon}_{i+1}^{p} = \dot{\lambda}(t) \bar{\varepsilon}^{p}(\underline{x}), \quad (45a) $$

$$ \Delta \sigma'_{i+1} = \lambda(t) \mathbb{C} \bar{\varepsilon}^p(\underline{x}). \quad (45b) $$
The corrective terms $\Delta \tilde{\sigma}_{i+1}$ and $\Delta \tilde{\varepsilon}_{i+1}$ are obtained from the weak form of eq. (38b) along with eq. (39b), i.e.

$$ \int_{[\Theta_{m+1}, \Theta_{m+1}+\Delta T] \times \Omega} \mathbf{C}(\Delta \tilde{\varepsilon}_{i+1} - \Delta \varepsilon_{i+1}^R) : \tilde{\varepsilon}(\underline{u}^*) \, d\Omega \, dt = 0, \quad \forall \underline{u}^* \in \mathcal{U}_0. \qquad (46) $$

Knowing $\Delta \varepsilon_{i+1}^R$ from the local stage, the corrective terms of strain and the corresponding stress tensor can be estimated using the operators $\mathbb{E}$ and $\mathbb{C}$ as

$$ \Delta \tilde{\varepsilon}_{i+1} = \mathbb{E} \Delta \varepsilon^R, \qquad (47a) $$

$$ \Delta \tilde{\sigma}_{i+1} = \mathbb{C} \Delta \varepsilon^R. \qquad (47b) $$

The initialisation phase has already introduced several PGD modes at “nodal cycle” $m+1$ from “nodal cycle” $m$, these space functions may first be re-used and only the time functions are updated. In case the reduced order model (ROM) hence obtained, is not accurate enough, an additional PGD pair can be added.

*Update of the time functions.* Therefore, first, considering that $\mu$ pairs have been used to approximate the quantities of interest, at the end of the global stage of LATIN iteration $i$, at LATIN iteration $i+1$ the objective is to re-use the spatial basis. The quantities of interest depending on plastic deformation are approximated as

$$ \begin{aligned} \Delta \tilde{\varepsilon}_{i+1}^p &= \sum_{j=1}^{\mu} \Delta \lambda_j(t) \tilde{\varepsilon}_j^p(\underline{x}), \\ \Delta \sigma'_{i+1} &= \sum_{j=1}^{\mu} \Delta \lambda_j(t) \mathbb{C} \tilde{\varepsilon}_j^p(\underline{x}). \end{aligned} \qquad (48) $$

The updates of the time functions $\{\lambda_j\}_{j=1}^\mu$ are evaluated from the minimisation of the constitutive relation error $e_{CE}^2$ defined for the given “nodal cycle” as the norm of the search direction equation given by eq. (22) as

$$ \begin{aligned} e_{CE}^2 &= \| \Delta \tilde{\varepsilon}_{i+1}^p - \mathbf{H}_{\sigma} \Delta \sigma'_{i+1} + \bar{\mathbf{\Lambda}}_{i+1} \|_{\mathbf{H}_{\sigma}^{-1}}^2, \\ &= \left\| \sum_{j=1}^{\mu} \Delta \lambda_j \tilde{\varepsilon}_j^p - \mathbf{H}_{\sigma} \sum_{j=1}^{\mu} \Delta \lambda_j \mathbb{C} \tilde{\varepsilon}_j^p + \bar{\mathbf{\Lambda}}_{i+1} \right\|_{\mathbf{H}_{\sigma}^{-1}}^2, \end{aligned} \qquad (49) $$

which leads to

$$ \{\Delta \lambda_j\}_{j=1}^\mu = \arg \min_{\{\Delta \lambda_i\}_{i=1}^\mu} \left\| \sum_{j=1}^{\mu} \Delta \lambda_j \tilde{\varepsilon}_j^p - \mathbf{H}_{\sigma} \sum_{j=1}^{\mu} \Delta \lambda_j \mathbb{C} \tilde{\varepsilon}_j^p + \bar{\mathbf{\Lambda}}_{i+1} \right\|_{\mathbf{H}_{\sigma}^{-1}}^2. \qquad (50) $$

This minimisation problem is solved by a discontinuous Galerkin method of order zero (see Relun et al., 2011). In case the correction provided by the update of the time functions is not satisfactory enough, a supplementary space-time pair will be added. The quality of the ROM is evaluated through an error saturation indicator, which is given by,

$$ \zeta = \frac{\xi_i - \xi_{i+1}}{\xi_i + \xi_{i+1}}. \qquad (51) $$

The approximation for this LATIN iteration is considered acceptable if the indicator is larger than a particular pre-defined value $\zeta^{\text{tol}}$. If the indicator is lower than $\zeta^{\text{tol}}$, the reduced-order basis is enriched with a new space-time PGD pair.
Addition of space-time modes. The enrichment of the reduced-order basis, performed through a hybrid strategy, involves Galerkin based temporal problem and minimal-residual based spatial problem.

The quantities of interest in the enrichment stage are approximated as

$$ 
\begin{aligned}
\Delta \dot{\boldsymbol{\varepsilon}}_{i+1}^p &= \lambda_{\mu+1}(t) \bar{\boldsymbol{\varepsilon}}_{\mu+1}^p(\underline{x}), \\
\Delta \boldsymbol{\sigma}'_{i+1} &= \lambda_{\mu+1}(t) \mathbb{C} \bar{\boldsymbol{\varepsilon}}_{\mu+1}^p(\underline{x}),
\end{aligned} 
\quad (52) $$

where $\{\dot{\lambda}_{\mu+1}(t), \bar{\boldsymbol{\varepsilon}}_{\mu+1}^p(\underline{x})\}$ is the new space-time pair. The search direction equation is now written as

$$ \Delta \dot{\boldsymbol{\varepsilon}}_{i+1}^p - \mathbf{H}_{\sigma} \Delta \boldsymbol{\sigma}'_{i+1} + \bar{\boldsymbol{\Lambda}}_{i+1} = 0, \quad (53) $$

with

$$ \bar{\boldsymbol{\Lambda}}_{i+1} = \mathbf{H}_{\sigma}(\hat{\boldsymbol{\sigma}}_{i+1/2} - \boldsymbol{\sigma}_{i+1}^{\text{up}}) - (\hat{\boldsymbol{\varepsilon}}_{i+1/2}^p - \dot{\boldsymbol{\varepsilon}}_{i+1}^{p, \text{up}}), \quad (54) $$

where $\dot{\boldsymbol{\varepsilon}}_{i+1}^{p, \text{up}}$ and $\boldsymbol{\sigma}_{i+1}^{\text{up}}$ are known from the update stage. The constitutive relation error can now be defined as

$$ 
\begin{aligned}
e_{CE}^2 &= \| \Delta \dot{\boldsymbol{\varepsilon}}_{i+1}^p - \mathbf{H}_{\sigma} \Delta \boldsymbol{\sigma}'_{i+1} + \bar{\boldsymbol{\Lambda}}_{i+1} \|^2_{\mathbf{H}_{\sigma}^{-1}}, \\
&= \| \dot{\lambda}_{\mu+1} \bar{\boldsymbol{\varepsilon}}_{\mu+1}^p - \mathbf{H}_{\sigma} \lambda_{\mu+1} \mathbb{C} \bar{\boldsymbol{\varepsilon}}_{\mu+1}^p + \bar{\boldsymbol{\Lambda}}_{i+1} \|^2_{\mathbf{H}_{\sigma}^{-1}}.
\end{aligned} 
\quad (55) $$

The space-time pair is solved using a fixed-point algorithm through successive iterations between the spatial and temporal problems till convergence is achieved. Knowing the spatial function $\bar{\boldsymbol{\varepsilon}}_{\mu+1}^p$, the temporal function $\lambda_{\mu+1}$ is calculated through the minimisation of the constitutive relation error $e_{CE}^2$

$$ \lambda_{\mu+1} = \arg \min_{\dot{\lambda}_{\mu+1}} \left\| \dot{\lambda}_{\mu+1} \bar{\boldsymbol{\varepsilon}}_{\mu+1}^p - \mathbf{H}_{\sigma} \lambda_{\mu+1} \mathbb{C} \bar{\boldsymbol{\varepsilon}}_{\mu+1}^p + \bar{\boldsymbol{\Lambda}}_{i+1} \right\|_{\mathbf{H}_{\sigma}^{-1}}^2. \quad (56) $$

This minimisation problem is also solved using zero-order discontinuous Galerkin technique.

For the spatial problem, with the knowledge of the calculated temporal function, a Galerkin-based formulation is proposed. All the superscripts and subscripts that have been used previously are dropped for simplicity for the following development. The strain partition relation combined with the state equation is written as

$$ \Delta \dot{\boldsymbol{\varepsilon}}' = \Delta \dot{\boldsymbol{\varepsilon}}^p + \mathbf{C}^{-1} \Delta \bar{\boldsymbol{\sigma}}', \quad (57) $$

where $\Delta \dot{\boldsymbol{\varepsilon}}^p$ is obtained from eq. (53) and $\Delta \boldsymbol{\sigma}' = \lambda \bar{\boldsymbol{\sigma}}$ with $\bar{\boldsymbol{\sigma}} = \lambda \mathbb{C} \bar{\boldsymbol{\varepsilon}}^p$. The kinematic admissibility condition in rate form is given by

$$ \int_{[\Theta_{m+1}, \Theta_{m+1}+\Delta T] \times \Omega} \Delta \dot{\boldsymbol{\varepsilon}}' : \boldsymbol{\sigma}^* \, d\Omega \, dt = 0, \quad \forall \boldsymbol{\sigma}^* \in \mathcal{S}_0, \quad (58) $$

which leads to

$$ \int_{\Omega} [\langle \mathbf{H}_{\sigma} \lambda^2 \rangle \bar{\boldsymbol{\sigma}} + \langle \dot{\lambda} \lambda \rangle \mathbf{C}^{-1} \bar{\boldsymbol{\sigma}} - \langle \bar{\boldsymbol{\Lambda}} \lambda \rangle] : \bar{\boldsymbol{\sigma}}^* \, d\Omega = 0, \quad \forall \bar{\boldsymbol{\sigma}}^* \in \mathcal{S}_0, \quad (59) $$

with $\langle \cdot \rangle = \int_{[\Theta_{m+1}, \Theta_{m+1}+\Delta T]} \cdot dt$ and $\boldsymbol{\sigma}^* = \lambda \bar{\boldsymbol{\sigma}}$. Introducing $\bar{\boldsymbol{\varepsilon}}$ such that

$$ \bar{\boldsymbol{\varepsilon}} = \mathbf{W}^{-1} \bar{\boldsymbol{\sigma}} - \bar{\boldsymbol{\delta}}, \quad (60) $$

with $\mathbf{W}^{-1} = \langle \mathbf{H}_{\sigma} \lambda^2 \rangle + \langle \dot{\lambda} \lambda \rangle \mathbf{C}^{-1}$ and $\bar{\boldsymbol{\delta}} = \langle \bar{\boldsymbol{\Lambda}} \lambda \rangle$ leads to

$$ \int_{\Omega} \bar{\boldsymbol{\varepsilon}} : \bar{\boldsymbol{\sigma}}^* \, d\Omega = 0, \quad \forall \bar{\boldsymbol{\sigma}}^* \in \mathcal{S}_0. \quad (61) $$
The static admissibility to zero of $\bar{\sigma}$ is written as

$$ \int_{\Omega} \bar{\sigma} : \varepsilon(\bar{\underline{u}}^*) d\Omega = 0, \quad \forall \bar{\underline{u}}^* \in \mathcal{U}_0, \qquad (62) $$

which may be reformulated as

$$ \int_{\Omega} W\varepsilon(\bar{\underline{u}}): \varepsilon(\bar{\underline{u}}^*) d\Omega = - \int_{\Omega} W\bar{\delta}: \varepsilon(\bar{\underline{u}}^*) d\Omega = 0, \quad \forall \bar{\underline{u}}^* \in \mathcal{U}_0, \qquad (63) $$

using eq. (60) and introducing displacement $\bar{\underline{u}} \in \mathcal{U}_0$ such that $\varepsilon(\bar{\underline{u}}) = \bar{\varepsilon}$. The pseudo-displacement $\bar{\underline{u}}$ along with the associated strain $\bar{\varepsilon}$ is obtained by solving the finite element problem as usual. Thereby the space function $\bar{\varepsilon}_{\mu+1}^p$ is estimated as,

$$ \bar{\varepsilon}_{\mu+1}^p = \frac{1}{\langle \lambda \lambda \rangle} [\bar{\varepsilon} - \langle \lambda \lambda \rangle C^{-1} W (\bar{\varepsilon} + \bar{\delta})]. \qquad (64) $$

The spatial basis vector thus obtained is orthonormalised with respect to the pre-existing basis vectors using Gram-Schmidt algorithm (see Ruhe, 1983).

The variables representing kinematic hardening, which are also cyclic, can be treated in a similar manner. With zero initial conditions, eq. (11) is solved over the “nodal cycle” $m + 1$ in a corrective form as

$$ \Delta \dot{\alpha}_{i+1} + \mathbf{H}_{\beta} \Lambda \Delta \alpha_{i+1} + \bar{\mathbf{\Pi}}_{i+1} = 0, \qquad (65) $$

with

$$ \bar{\mathbf{\Pi}}_{i+1} = \mathbf{H}_{\beta} (\Lambda \alpha_i - \hat{\beta}_{i+1/2}) + \dot{\alpha}_i - \hat{\alpha}_{i+1/2}, \qquad (66) $$

and

$$ \Delta \alpha_{i+1} = \alpha_{i+1} - \alpha_i, \quad \Delta \dot{\alpha}_{i+1} = \dot{\alpha}_{i+1} - \dot{\alpha}_i. \qquad (67) $$

It must be noted that the calculation of $\mathbf{H}_{\beta}$ and $\mathbf{H}_{\sigma}$ is computationally expensive. The idea herein is to calculate these operators only once at the first LATIN iteration for a given “nodal cycle” $m + 1$, and re-use them for the rest of the iterations.

The convergence between local and global solutions is evaluated independently for every “nodal cycle” by a LATIN error indicator defined as

$$ \xi = \frac{\| \hat{s}_{i+1/2}^p - s_{i+1}^p \|}{\| \hat{s}_{i+1/2}^p \| + \| s_{i+1}^p \|}, \qquad (68) $$

with

$$ \|s^p\|^2 = \int_{[\Theta_m, \Theta_{m+1} + \Delta T] \times \Omega} (\sigma : \mathbf{H}_{\sigma} \sigma + \beta : \mathbf{H}_{\beta} \beta + \dot{\varepsilon}^p : \mathbf{H}_{\sigma}^{-1} \dot{\varepsilon}^p + \varepsilon^e : \mathbf{C} \varepsilon^e + \dot{\alpha} : \mathbf{H}_{\beta}^{-1} \dot{\alpha}) d\Omega dt. \qquad (69) $$

This approach has been developed to reduce drastically the computational cost compared to classical LATIN technique, and hence enables to compute damage evolution for fatigue damage cases. The performance of the method is experienced in the next section on academic examples.

## 5. Numerical examples

The robustness and abilities of the innovative algorithm are tested on some two-dimensional mechanical problems. For the two-time scale computation, the search direction operators are calculated only once for every “nodal cycle”. For simple geometry and boundary conditions, the ROB obtained from the “training phase” is generally adequate, and the calculation of the “nodal cycles” involves only updating the ROB. For cases involving damage localisations, only updating the basis is not sufficient, therefore the ROB should be enriched by a maximum of one PGD mode per “nodal cycle” if $\zeta < 0.01$. All the two dimensional structures mentioned herein are discretised in space using linear four-noded isoparametric quadrilateral elements with four Gauss points per element.
5.1. Verification with mono-scale LATIN method

Verification analyses are carried out for qualitative estimation of the gain in CPU time and loss in accuracy compared to mono-scale LATIN-PGD framework. For the first set of analysis, a material behaviour characterised by viscoplasticity and kinematic hardening without any damage is considered. The structure used for this verification is a rectangular plate which is constrained on three sides as shown in fig. 4. The geometry of the plate is described by its length L = 50 mm and width W = 100 mm. A sinusoidal displacement of the form $U_d = U_0 \sin \frac{2\pi t}{\Delta T}$, with time period $\Delta T = 1$ s is applied uniformly on the unconstrained side of the plate for 1000 cycles.

The choice of simple structure for verification analysis is to ensure that there are no localisations, which may generate a large number of PGD modes that would lead to unreasonable CPU time for the mono-scale reference solution. The PGD modes obtained after the “training stage” are considered to be optimal as the structure does not incite any localisation and the space functions obtained after the “training stage” are enough to describe the spatial dependency of the quantities of interest.

<figure><img src="image_5.png" /><figcaption>Figure 4: Rectangular plate with different yield stresses</figcaption></figure>

The structure is assumed to be made of Cr-Mo steel at 25°C with properties given in table 1. The material at the central part of the plate is considered to be of yield stress $\sigma_y = 188$ MPa and the peripheral material is of $\sigma_y = 189$ MPa.

<table><caption>Table 1: Material properties for Cr-Mo steel at 20°C (see Lemaitre and Desmorat, 2005)</caption><tr><td>E</td><td>199740 MPa</td></tr><tr><td>ν</td><td>0.3</td></tr><tr><td>C</td><td>20602 MPa</td></tr><tr><td>a</td><td>141.4 MPa</td></tr><tr><td>K</td><td>195 MPa s<sup>1/n</sup></td></tr><tr><td>n</td><td>12.5</td></tr><tr><td>k</td><td>K<sup>-n</sup></td></tr><tr><td>h</td><td>0.2</td></tr></table>

The problem such defined is initially solved for $U_0 = 0.073$ mm using classical mono-scale LATIN-PGD method. Using the two-scale computation, four temporal discretisations are examined, with a minimum of 50 cycles to a maximum of 500 cycles per time element (see table 2).

As plasticity is the defining phenomenon in the process, the relative error is considered to be based on the plastic strain rate, i.e.

$$er^2 = \frac{\int_{\Omega \times [0,T]} (\dot{\varepsilon}_{ms}^p - \dot{\varepsilon}_{ts}^p) : (\dot{\varepsilon}_{ms}^p - \dot{\varepsilon}_{ts}^p) \, d\Omega \, dt}{\int_{\Omega \times [0,T]} \dot{\varepsilon}_{ms}^p : \dot{\varepsilon}_{ms}^p \, d\Omega \, dt + \int_{\Omega \times [0,T]} \dot{\varepsilon}_{ts}^p : \dot{\varepsilon}_{ts}^p \, d\Omega \, dt} \quad (70)$$
<table><caption>Table 2: Different temporal discretisations</caption><thead><tr><th>Temporal discretisation</th><th>Number of cycles per time element</th><th>Total number of time elements</th></tr></thead><tbody><tr><td>I</td><td>50</td><td>20</td></tr><tr><td>II</td><td>100</td><td>10</td></tr><tr><td>III</td><td>200</td><td>5</td></tr><tr><td>IV</td><td>500</td><td>2</td></tr></tbody></table>

where $\dot{\epsilon}_{ms}^p$ is the mono-scale reference solution and $\dot{\epsilon}_{ts}^p$ is the two-scale solution. Expected result is seen in fig. 5, i.e. enlarging the size of time elements decreases the accuracy of the approximation and hence increases $er$. Also, there is a very minimal loss of accuracy for the two-scale computation even if 500 cycles per time element is used. The relative CPU time is computed as $t_{ts}^{cal}/t_{ms}^{cal}$, where $t_{ms}^{cal}$ is the time needed for the mono-scale computation and $t_{ts}^{cal}$ is the calculation time for the two-scale computation. As expected, this quantity decreases with increase in number of cycles per time element. It should also be noted from fig. 5 that there is a drastic drop in the calculation time for two-scale computation, compared to the mono-scale computational time, even if 50 cycles per time element are used.

<figure><img src="image_6.png" /><figcaption>Figure 5: Relative accuracy and saved time for different temporal meshes</figcaption></figure>

For visco-plasticity, bereft of damage, the two-scale solutions provide very good accuracy with extremely low run-time as compared to mono-scale reference solution.

To investigate the influence of damage on the two-time scale method, similar tests are performed for material that is susceptible to damage along with visco-plasticity and kinematic hardening. The damage threshold considered is $(1.4\sigma_y)^2/2E$. Parameters $k_D = 0.128 \text{ MPa}^{-n_D} \text{s}^{-1}$ and $n_D = 2$ are chosen to describe the damage evolution process. The rest of the material properties are considered to be the same as before, i.e. given in table 1.

The problem is solved for $U_0 = 0.070$ mm, initially using mono-scale LATIN-PGD method to obtain a reference solution, and then using the four temporal discretisations as mentioned in the previous case. The relative error estimation with respect to the reference mono-scale solution is considered to be dictated by damage and given by

$$er_d^2 = \frac{\int_{\Omega \times [0, T]} (D_{ms} - D_{ts}) \cdot (D_{ms} - D_{ts}) \, d\Omega \, dt}{\int_{\Omega \times [0, T]} D_{ms} \cdot D_{ms} \, d\Omega \, dt + \int_{\Omega \times [0, T]} D_{ts} \cdot D_{ts} \, d\Omega \, dt}, \quad (71)$$

where $D_{ms}$ and $D_{ts}$ are obtained from a mono-scale and a two-scale LATIN-PGD algorithms respectively. A strong increase in the relative error is seen with respect to the increase in the size of the “time elements” (see fig. 6). A minimum of 2.6% and a maximum of 4.2% are obtained, which can be accepted as reasonable. The reduction in CPU time is however much more drastic compared to the previous case.
<figure><img src="image_7.png" /><figcaption>Figure 6: Accuracy and saved time for different temporal meshes with respect to a mono-scale reference solution</figcaption></figure>

Further tests are conducted to analyse the effect of the damage level on the accuracy of different temporal discretisations. The two-scale LATIN-PGD algorithm is used to calculate the solution for four different load amplitudes $U_0$. For each load level, the solution is estimated using all the four discretisations given in table 2. Instead of calculating mono-scale solutions, the solutions obtained using discretisation *I* (see table 2) are considered to be the references. It has to be mentioned that all the two-scale LATIN analyses included “training stages” consisting of 10 cycles and the ROBs obtained from the “training stages” are not enriched, only the time functions were updated at each “nodal cycle”.

It is evident from fig. 7, where the damage evolution obtained using discretisation *I* for different load levels is plotted, that the increase in load amplitude induces higher damage on both parts of the structure.

<figure><img src="image_8.png" /><figcaption>Figure 7: Damage evolution in the plate for different load levels using discretisation *I*</figcaption></figure>

From fig. 8a it can be inferred that with increase in load level the coarser time discretisation results in higher inaccuracy. The reason behind is the fact that the damage behaviour becomes non-linear with higher load level (see fig. 7). For lower load level, however, the damage level is low and the behaviour is found relatively linear, and coarser temporal mesh can also be used to obtain acceptable accuracy.

For a particular load level, if finer temporal mesh is used, the CPU time for the computation is higher (see fig. 8b). Also, for a given temporal discretisation, higher load level induces higher material non-linearities and more LATIN iterations are needed at each “nodal cycle”, resulting in higher CPU time.

These observations lead to a conclusion that the magnitude of damage is a defining factor in the choice of the
<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison with the size of time elements relative to discretisation *I* solution</figcaption></figure>

type of temporal discretisation. For low damage level, coarse temporal mesh can be used for reasonable accuracy while improving the CPU cost. However, for higher damage level, a fine temporal mesh must be used which will give adequate accuracy in the expense of calculation time. Therefore, an optimal choice would be to use an adaptive temporal discretisation (e.g. table 3), where the length of the time elements will be dictated by the damage level.

5.2. Influence of the “training stage”

However, before using any adaptive temporal scheme, the influence of the number of cycles in the “training stage” is also checked. Using $U_0 = 0.063$ mm for $1 \times 10^5$ cycles in the structure given in fig. 4, four virtual tests are conducted with 50, 20, 10, and 5 cycles for the “training stage”. Thereafter, a uniform temporal mesh of 100 cycles per time element is considered for all the four cases. The solution obtained using 50 cycles in the “training stage” is considered as the reference, and solutions of the other three cases are compared using the error criterion of eq. (71). This error criterion increases as the number of cycles in the “training stage” decreases (see fig. 9). The inaccuracy with respect to the reference solution, is however extremely low (see fig. 9), and using higher number of cycles in the “training stage” is not necessary and can also be detrimental as far as the numerical cost is concerned.

<figure><img src="image_10.png" /><figcaption>Figure 9: Variation of accuracy with respect to the number of cycles of the “training stage”</figcaption></figure>

The distribution of damage obtained after the end of loading using 50 cycles in the “training stage” is shown in fig. 10, with maximum damage at the central part and minimum at the periphery. The corresponding evolutions of damage for these two regions are given in fig. 11.
<figure><img src="image_11.png" /><figcaption>Figure 10: Distribution of damage at the end of loading, i.e. after $1 \times 10^5$ cycles</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Damage evolution in the plate for a total of $1 \times 10^5$ cycles</figcaption></figure>
## 5.3. Simulation of large number of cycles

With the effect of the “training stage” being sorted, the aim now is to observe the damage behaviour involving large number of cycles. An adaptive temporal meshing scheme given in table 3 is used for the numerical simulation of the simple structure given in fig. 4. The choice of this type of simple structure is to avoid any localisation and stress raisers that may lead to regional concentration of damage.

<table><caption>Table 3: Adaptive temporal discretisation scheme with respect to damage level</caption><thead><tr><th>Damage range</th><th>Number of cycles per time element</th></tr></thead><tbody><tr><td>$D_m^{max} < 0.001$</td><td>500</td></tr><tr><td>$0.001 < D_m^{max} \le 0.01$</td><td>200</td></tr><tr><td>$0.01 < D_m^{max} \le 0.1$</td><td>100</td></tr><tr><td>$D_m^{max} > 0.1$</td><td>50</td></tr></tbody></table>

The size of the time element [$\Theta_m, \Theta_{m+1}$] i.e. the value of $\Theta_{m+1}$ is estimated using the maximum damage value at the end of “nodal cycle” $m$ denoted by $D_m^{max}$, according to the scheme presented in table 3.

Using a “training stage” comprising of only 10 cycles, a total of $5 \times 10^5$ cycles are simulated with load amplitude $U_0 = 0.060$ mm. The adaptive scheme calculates only 2185 “nodal cycles” where the ROB obtained from the “training stage” containing three PGD modes, is only updated.

<figure><img src="image_13.png" /><figcaption>Figure 12: Distribution of damage at the end of loading, i.e. after $5 \times 10^5$ cycles</figcaption></figure>

The distribution of damage at the end of loading is observed to be much higher in the central part where the material is weaker, compared to the periphery (see fig. 12). The corresponding damage evolution is depicted in fig. 13, where the central part shows non-linear increase in damage, and the peripheral part shows damage saturation after certain load cycles.

## 5.4. Pre-damaged structure

The adaptive temporal scheme of the two-time scale discretisation is now used to examine the damage behaviours of virgin and pre-damaged structures. The two-dimensional geometry consists of a rectangular plate of length $L = 20$ mm, width $W = 10$ mm, with two semicircular notches of equal radii $\phi = 4$ mm as illustrated in fig. 14. Uniformly distributed sinusoidal displacements of amplitude $U_0 = 0.0057$ mm and time period $\Delta T = 1$ s are applied on the longitudinal ends for 20000 cycles. In fig. 14 is also shown the quarter of the plate generated out of symmetry and the region of interest where damage localisation is anticipated. Unlike the previous cases, there is localisation of the quantities of interest and the ROB is allowed to enrich by a maximum of one PGD mode per “nodal cycle”.

For the first analysis, the simulation is conducted assuming the complete structure to be virgin. The “training stage”, consisting of 10 cycles, generates 7 PGD modes. A total of 171 “nodal cycles” are subsequently calculated generating a total of 11 PGD modes, i.e. 4 supplementary modes are generated after the “training stage”.
<figure><img src="image_14.png" /><figcaption>Figure 13: Damage evolution in the plate for a total of 5 × 10⁵ cycles</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: A plate with semicircular notches subjected to distributed loads</figcaption></figure>

The next analysis is to assign certain Gauss points (as shown in fig. 16) initial damage values of 0.005, and then recalculate the problem. The complete behaviour has been simulated using a total of 186 “nodal cycles” with a total of 13 PGD modes. The “training stage” similar to the previous case is of 10 cycles, generating 7 PGD modes.

<figure><img src="image_16.png" /><figcaption>Figure 15: FE discretisation showing virgin(●) and pre-damaged (■) Gauss points</figcaption></figure>

The spread of damage in the region of interest is shown in fig. 16 for both virgin and pre-damaged structures after certain specific cycles. For both cases the most affected Gauss point is identical, the respective damage evolution is depicted in fig. 17.

5.5. *Variable amplitude loading*

The focus now is to simulate, using the two-time scale adaptive discretisation scheme, the different behaviour of a structure when subjected to constant amplitude (CA) load and variable amplitude (VA) load. For the VA load, the effect of load sequence, i.e. high to low (H-L) load and low to high (L-H) load, is also examined.
<figure><img src="image_17.png" /><figcaption>Figure 16: Damage distribution at the end of loading for the notched plate</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Damage evolution at the weakest Gauss point for the notched plate</figcaption></figure>
A square plate of length $L = 40$ mm with an elliptical hole with semi-major axis $a = 10$ mm and semi-minor axis $b = 5$ mm is used for the analysis (see fig. 18). Sinusoidal displacements of the form $U_d = U_0 \sin \frac{2\pi t}{\Delta T}$ with time period $\Delta T = 1$ s are uniformly applied on the ends, as illustrated in fig. 18 for 10000 cycles. The quarter of the plate along with the symmetric boundary conditions and the anticipated region of interest is also represented in fig. 18.

<figure><img src="image_19.png" /><figcaption>Figure 18: A square plate with an elliptical hole subjected to distributed loads</figcaption></figure>

The three cases of loading, which are constant amplitude (CA), uniformly increasing (L-H) and uniformly decreasing (H-L) are depicted on fig. 19. The three cases are such that the mean values of the load amplitudes are the same.

<figure><img src="image_20.png" /><figcaption>Figure 19: Variation of amplitude with number of load cycles for three different load cases</figcaption></figure>

For these three loading cases, the spread of damage is shown in fig. 20, at certain cycles, in the region of interest. The L-H loading has the minimum damage and the H-L loading has the maximum, with the CA loading being in between. In fig. 21 the evolution of damage for the weakest GP is plotted for three load cases. It can be observed that the evolution of damage is much higher for H-L than for L-H loading, whereas both cases have the same mean load amplitude. The damage curve for the CA loading is enveloped by the two VA loading cases.

In terms of numerical behaviour of the algorithm, in the three cases, the “training stage” consists of 10 cycles. It yields 8 PGD pairs for CA and H-L case, but 7 PGD pairs for L-H case. After the “training stage”, the adaptive scheme is employed, it leads to the computation of 95 “nodal cycles” for CA application, 90 “nodal cycles” for L-H and 98 for H-L. As mentioned before, a maximum of one space-time pair may be added per “nodal cycle”. This leads to the extension of the basis to 19, 18, and 19 PGD modes for CA, L-H, and H-L, which means that approximately 10 supplementary modes have been computed for all the cases.
<figure><img src="image_21.png" /><figcaption>Figure 20: Spread of damage at certain load cycles in the region of interest for the plate with elliptical hole</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21: Evolution of damage with respect to number of load cycles for the weakest GP of plate with elliptical hole</figcaption></figure>
## 6. Conclusion

An innovative technique has been proposed in this article for the simulation of fatigue problems. The time domain has been discretised in a finite element like description using two time scales, a short time that describes the rapid evolution of the quantities of interest within one cycle and a long time that indicates their slow evolution along the cycles. The description of the time domain introduced certain pre-defined cycles called the “nodal cycles”, where the quantities of interest are calculated and from which an interpolation is used within each time element. The usage of LATIN-PGD technique to calculate the quantities of interests at the “nodal cycles” has been proved to be robust and cost effective as shown in the academic examples, with applications up to half a million cycles. A non-homogeneous long time discretisation has been investigated to obtain higher accuracy for lower CPU cost. This method has on a whole proved effective to deal with material non-linearities for large number of cycles and could be used in the future for random or even dynamic fatigue simulations.

## Acknowledgement

The authors are grateful to the German Research Foundation (DFG) for funding the research through International Research Training Group 1627. The funding from the German-French University (DFH-UFA) is also acknowledged.

## Appendix A. Finite element operators

The construction of the FE operators $E$ and $C$ used in section 4.3.3 is here briefly summarised. For the derivation of the finite element operator $E$, the spatial problem eq. (43) is considered, i.e.

$$ \int_{\Omega} C\bar{\varepsilon} : \varepsilon(\bar{u}^*) d\Omega = \int_{\Omega} C\bar{\varepsilon}^p : \varepsilon(\bar{u}^*) d\Omega, \quad \forall \bar{u}^* \in \mathcal{U}_0. \qquad (A.1) $$

Classical finite element scheme is used for the spatial discretisation. The left hand side of eq. (A.1) for one elemental domain $\Omega_{el}$ can be written as

$$ \left[ \int_{\Omega_{el}} \mathbf{BCB}^T d\Omega \right] \{\bar{u}_{el}\}, \qquad (A.2) $$

and the corresponding right hand side reads

$$ \left[ \int_{\Omega_{el}} \mathbf{BC} d\Omega \right] \{\bar{\varepsilon}_{el}^p\}, \qquad (A.3) $$

where $\{\bar{u}_{el}\}$ and $\{\bar{\varepsilon}_{el}^p\}$ denote the nodal unknowns and the unknowns at all Gauss points, for one element respectively. The matrix $\mathbf{B}$ comprises of the derivatives of the shape functions. Assembling the elemental equations leads to the following global problem

$$ \mathbf{K}_{tot} \{\bar{u}_{tot}\} = \mathbf{L}_{tot} \{\bar{\varepsilon}_{tot}^p\}, \qquad (A.4) $$

where $\{\bar{u}_{tot}\}$ and $\{\bar{\varepsilon}_{tot}^p\}$ are the total number of nodal unknowns and Gauss point quantities respectively, $\mathbf{K}_{tot}$ and $\mathbf{L}_{tot}$ are the assembled matrices. An operator

$$ \mathbf{O}_{tot}^{red} = (\mathbf{K}_{tot}^{red})^{-1} \mathbf{L}_{tot}^{red} \qquad (A.5) $$

thereafter is calculated. $\mathbf{K}_{tot}^{red}$ and $\mathbf{L}_{tot}^{red}$ denote the resulting assembled matrices after the removal of the rows and columns corresponding to the known nodal degrees of freedom. The corresponding rows of $\mathbf{O}_{tot}^{red}$ are thereby filled with zeros to obtain $\mathbf{O}_{tot}$. Thus, eq. (A.4) can be re-written as

$$ \{\bar{u}_{tot}\} = \mathbf{O}_{tot} \{\bar{\varepsilon}_{tot}^p\}. \qquad (A.6) $$
The derivatives of the shape functions contained in the total assembled matrix are used to express eq. (A.6) as

$$B_{tot}^{T} \{\bar{u}_{tot}\} = B_{tot}^{T} O_{tot} \{\bar{\varepsilon}_{tot}^{p}\} \Rightarrow \{\bar{\varepsilon}_{tot}\} = E \{\bar{\varepsilon}_{tot}^{p}\}, \qquad (A.7)$$

with E being the finite element operator that relates the space functions corresponding to plastic strain {$\bar{\varepsilon}_{tot}^p$} and total strain {$\bar{\varepsilon}_{tot}$}.

Using the strain partition relationship,

$$\{\bar{\varepsilon}_{tot}^{e}\} = \{\bar{\varepsilon}_{tot}\} - \{\bar{\varepsilon}_{tot}^{p}\} \qquad (A.8)$$

the space functions related to stress can be expressed using eq. (A.7) as

$$\{\bar{\sigma}_{tot}\} = C_{tot}(E - I) \{\bar{\varepsilon}_{tot}^{p}\} = C \{\bar{\varepsilon}_{tot}^{p}\}, \qquad (A.9)$$

with $C_{tot}$ being the complete elasticity matrix for all the Gauss points, I the identity matrix, and C the FE operator relating the space function of plastic strain {$\bar{\varepsilon}_{tot}^p$} and stress {$\bar{\sigma}_{tot}$}.

## References

Aguado, J. V., Huerta, A., Chinesta, F., Cueto, E., 2015. Real-time monitoring of thermal processes by reduced-order modeling. International Journal For Numerical Methods In Engineering 102 (5), 991–1017.

Allix, O., Ladevèze, P., Gilletta, D., Ohayon, R., 1989. A damage prediction method for composite structures. International Journal for Numerical Methods in Engineering 27 (2), 271–283.

Ammar, A., Zghal, A., Morel, F., Chinesta, F., 2015. On the space-time separated representation of integral linear viscoelastic models. Comptes Rendus Mécanique 343 (4), 247–263.

Aubry, D., Puel, G., 2010. Two-timescale homogenization method for the modeling of material fatigue. IOP Conference Series: Materials Science and Engineering 10 (1), 012113.

Bellenger, E., Bussy, P., 1998. Plastic and viscoplastic damage models with numerical treatment for metal forming processes. Journal of Materials Processing Technology 80–81, 591–596.

Bhamare, S., Eason, T., Spottswood, S., Mannava, S. R., Vasudevan, V. K., Qian, D., 2014. A multi-temporal scale approach to high cycle fatigue simulation. Computational Mechanics 53 (2), 387–400.

Bhattacharyya, M., Fau, A., Nackenhorst, U., Néron, D., Ladevèze, P., 2017. A latin-based model reduction approach for the simulation of cycling damage. Computational Mechanics.

Bhattacharyya, M., Fau, A., Nackenhorst, U., Néron, D., Ladevèze, P., 2018. A model reduction technique in space and time for fatigue simulation. In: Sorić, J., Wriggers, P., Allix, O. (Eds.), Multiscale Modeling of Heterogeneous Structures. Springer International Publishing, Cham, pp. 183–203.

Chaboche, J. L., Lesne, P. M., 1988. A non-linear continuous fatigue damage model. Fatigue and Fracture of Engineering Materials and Structures 11 (1), 1–17.

Chinesta, F., Ammar, A., Cueto, E., 2010a. Proper generalized decomposition of multiscale models. International Journal for Numerical Methods in Engineering 83 (8-9), 1114–1132.

Chinesta, F., Ammar, A., Cueto, E., 2010b. Recent advances and new challenges in the use of proper generalized decomposition for solving multidimensional models. Archives of Computational Methods In Engineering 17 (4), 327–350.

Chinesta, F., Ladevèze, P., Cueto, E., 2011. A short review on model order reduction based on proper generalized decomposition. Archives of Computational Methods in Engineering (18), 395–404.

Cognard, J., Ladeveze, P., Talbot, P., 1999. A large time increment approach for thermo-mechanical problems. Advances in engineering software 30, 583–593.

Cognard, J.-Y., Ladevèze, P., 1993. A large time increment approach for cyclic viscoplasticity. International Journal of Plasticity 9 (2), 141–157.

Cojocaru, D., Karlsson, A., 2006. A simple numerical method of cycle jumps for cyclically loaded structures. Internation Journal of Fatigue 28, 1677–1689.

De Vuyst, F., 2016. Efficient solvers for time-dependent problems: a review of imex, latin, paraexp and parareal algorithms for heat-type problems with potential use of approximate exponential integrators and reduced-order models. Adv. Model. and Simul. 3 (8).

Desmorat, R., Cantournet, S., 2008. Modeling microdefects closure effect with isotropic/anisotropic damage. International Journal of Damage Mechanics 17 (1), 65–96.

Desmorat, R., Kane, A., Seyed, M., Sermage, J., 2007. Two scale damage model and related numerical issues for thermo-mechanical High Cycle Fatigue. European Journal of Mechanics - A/Solids 26 (6), 909–935.

Desmorat, R., Pauget, F., Sermage, J., 2006. Damage 2005: a post-processor for high cycle fatigue under complex thermomechanical loading. In: ASME Pressure Vessels and Piping Conference. Vancouver, Canada.

Fish, J., Shek, K., Pandheeradi, M., Shephard, M. S., 1997. Computational plasticity for composite structures based on mathematical homogenization: Theory and practice. Computer Methods in Applied Mechanics and Engineering 148 (1–2), 53–73.

Gander, M. J., 2015. 50 Years of Time Parallel Time Integration. Springer International Publishing, pp. 69–113.

González, D., Alfaro, I., Quesada, C., Cueto, E., Chinesta, F., 2015. Computational vademecums for the real-time simulation of haptic collision between nonlinear solids. Computer Methods in Applied Mechanics and Engineering 283 (1), 210–223.
Guennouni, T., 1988. Sur une méthode de calcul de structures soumises à des chargements cycliques : l'homogénéisation en temps. Modélisation mathématique et analyse numérique 22 (3), 417-455.

Guennouni, T., Aubry, D., 1986. Réponse homogénéisée en temps de structures sous chargements cycliques. Comptes rendus de l'Académie des sciences. Série 2, Mécanique, Physique, Chimie, Sciences de l'Univers, Sciences de la Terre 303 (20), 1765-1768.

Ladevèze, P., 1989. The large time increment method for the analyse of structures with nonlinear constitutive relation described by internal variables. Comptes Rendus Académie des Sciences Paris 309 (II), 1095-1099.

Ladevèze, P., 1999. Nonlinear computational structural mechanics. Mechanical Engineering Series. Springer New York.

Ladevèze, P., 2014. Separated Representations and PGD-Based Model Reduction: Fundamentals and Applications. Springer Vienna, Ch. PGD in linear and nonlinear Computational Solid Mechanics, pp. 91-152.

Ladevèze, P., 2016. On reduced models in nonlinear solid mechanics. European Journal of Mechanics - A/Solids 60, 227-237.

Ladevèze, P., Lemaitre, J., 1984. Damage effective stress in quasi-unilateral conditions. In: The 16th international congress of theoretical and applied mechanics.

Ladevèze, P., Lorong, P., 1991. A large time increment approach with domain decomposition for mechanical non linear problem. In: Proceedings of the 10th International Conference on Computing Methods in Applied Sciences and Engineering on Computing Methods in Applied Sciences and Engineering. Nova Science Publishers, Inc., Commack, NY, USA, pp. 569-578.

Lemaitre, J., 1996. A course on damage mechanics. Springer, Berlin.

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

Lemaitre, J., Sermage, J., Desmorat, R., Apr 1999. A two scale damage concept applied to fatigue. International Journal of Fracture 97 (1), 67.

Lesne, P. M., Savalle, S., 1989. An efficient cycles jump technique for viscoplastic structures involving large number of cycles. In: Second International Conference on Computational Plasticity. Barcelona, pp. 591-602.

Néron, D., Ladevèze, P., 2010. Proper generalized decomposition for multiscale and multiphysics problems. Archives of Computational Methods In Engineering 17, 351-372.

Nesnas, K., Saanouni, K., 2000. A cycle jumping scheme for numerical integration of coupled damage and viscoplastic models for cyclic loading path. Revue Européenne des Éléments Finis 9, 865-891.

Nouy, A., 2010. A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations. Computer Methods in Applied Mechanics and Engineering 199 (23), 1603-1626.

Oskay, C., Fish, J., 2004. Fatigue life prediction using two-scale temporal asymptotic homogenization. International Journal for Numerical Methods in Engineering 61, 329-359.

Puel, G., Aubry, D., 2014. Efficient fatigue simulation using periodic homogenization with multiple time scales. Journal for Multiscale Computational Engineering 12 (4), 291-318.

Relun, N., Néron, D., Boucard, P.-A., 2011. Multiscale elastic-viscoplastic computational analysis. European Journal of Computational Mechanics 20 (7-8), 379-409.

Ruhe, A., 1983. Numerical aspects of gram-schmidt orthogonalization of vectors. Linear Algebra and its Applications 52, 591-601.

Sobczyk, K., Spencer, B. F., 1992. Random Fatigue: From Data to Theory. Academic Press.

Suresh, S., 2001. Fatigue of materials, 2nd Edition. Cambridge University Press, Cambridge.

Vandoren, B., Proft, K. D., Simone, A., Sluys, L., 2013. A novel constrained large time increment method for modelling quasi-brittle failure. Computer Methods in Applied Mechanics and Engineering 265, 148-162.