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

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.

A multi-temporal scale model reduction approach for the computation of fatigue damage 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(Desmorat et al., , 2006Desmorat 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 4 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.
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, highfrequency 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 i before investigating the subsequent time t i+1 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 problems 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.

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 Ω having a boundary ∂Ω. This structure may be subjected to volumetric forces f d , prescribed displacements U d over a part ∂ 1 Ω of the boundary ∂Ω and to prescribed traction forces F d over the complementary part ∂ 2 Ω of the boundary ∂Ω. The field of internal stresses thus generated is permissible with the applied external forces if the static admissibility condition defined as ∀u * ∈ U * is satisfied. Here σ denotes the statically admissible stress tensor, and U * is the vector space associated with the space U of kinematically admissible displacement fields 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 ∀σ * ∈ F * where F * is the vector space associated with the space F of statically admissible fields σ. Under the assumption of infinitesimal strain theory, the kinematically admissible total strain tensor field ε can be additively decomposed into an elastic part ε e and a plastic part ε 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) where the Macauley brackets • represents the non-negative part of •. Equation (3a) is the elastic state law where E is the modulus of elasticity, ν 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 α i j and its corresponding thermodynamic force β i j 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), From the damage function f d and the damage multiplierλ 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λ 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 where with σ D i j being the deviatoric part of the stress tensor, σ y being the yield stress, and a is a material parameter for kinematic hardening. Contrary to plasticity,λ p has a closed form solution for visco-plasticity, which is defined by the Norton's law aṡ 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 aṡ 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 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.

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 = ε p , ε e ,α,Ḋ, σ, β, 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 Γ 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 ∈ A d and the approximation at the local stageŝ i+1/2 ∈ Γ continues till convergence is reached. The convergence, theoretically, is obtained when the exact solution s ex ∈ A d ∩ Γ can be calculated. However, in a numerical point of view, the convergence criterion is based on the distance between s i+1 andŝ i+1/2 , which should be lower than a pre-defined tolerance to achieve convergence.

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 ∈ A d , the estimation ofŝ i+1/2 ∈ Γ satisfies the search direction equation where the direction of ascent B + is assumed to be vertical , i.e. B + −1 = 0 (see Ladevèze, 1999). The approximation 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ŝ i+1/2 is calculated, it is used thereafter to solve the global stage to obtain s i+1 .

Global stage including PGD
The solution set s i+1 ∈ 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 The operator H − lies in the tangent space to the solution set s i+1/2 in the manifold Γ, C denotes the Hooke tensor for the undamaged material. It should be noted that H − is of 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 The search direction equation for hardening variables along with the state equation for kinematic hardening can be expressed as where Λ 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 To begin with the developments, consider the search direction equation relating elastic strain and stress written in terms of corrections, with the residual stress term ∆R i+1 is given by Introducing a residual strain term ∆ε R i+1 given by eq. (12) can be represented as, Using the additive strain decomposition relationship (infinitesimal strain theory), i.e.
eq. (15) is transformed into 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 where both the quantities have to be kinematically admissible to zero. Plugging this into eq. (17) such that with the complete description of stress is given by, The quantities ∆σ ′ i+1 , ∆ε ′ i+1 , ∆ε p i+1 will satisfy the static admissibility through PGD representation, and the quantities ∆σ i+1 , ∆ε i+1 , will satisfy the static admissibility bereft of PGD as ∆ε R i+1 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 where∆ i+1 provides information from the local stage and is given as 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 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 H σ and H β are given in Bhattacharyya et al. (2017) as, with J 2 = 3 2 ̺ · ̺ and I being the identity matrix.
PGD is then used to estimate the correction terms corresponding to the plastic behaviour ∆σ ′ i+1 and ∆ε p i+1 in a separable form.

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, 1999Ladevèze, , 2014. Considering any separable quantity of interest defined over the complete space-time domain represented by χ 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 ι is the number of terms chosen to approximate the quantity of interest, χ can be written as 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, ∆σ ′ and ∆ε 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.

Updating stage
Considering µ 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 . If this updating stage does not improve the accuracy of the approximation significantly, the greedy algorithm allows to enlarge the space and time bases.

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 whereλ µ+1 ,ε p µ+1 , andσ µ+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.

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.

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θ is defined on the domain [0, T ] representing the beginning of every cycle, and a fine temporal discretisation τ defined within a cycle. Any instant of interest t can be discretised by the pair θ , τ as where N is the number of cycles included in the time domain.

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" [Θ m , Θ m+1 ] containingp cycles and demarcated by the "nodal cycles" m and m + 1, a local variable θ k , ∀k ∈ 0,p − 1 is introduced which represents the coarse discretisation for the given "time element". This representation provides the time t to be defined as For the given "time element" [Θ m , Θ m+1 ] temporal shape functions ν m and ν m+1 are defined as nodal cycle m + 1 nodal cycle m If any temporal quantity of interest χ is known at the "nodal cycles" m and m + 1, the approximation on the whole "time element" [Θ m , Θ m+1 ], containingp cycles is given by with τ m = τ m+1 = τ k belonging to [0, ∆T ] as the time period is constant and k = {0, 1, · · · ,p − 1}. The discretisation is defined such that θ k = Θ m + k∆T , θ 0 = Θ m and θp −1 = Θ 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". 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 [Θ 0 , Θ 0 + ∆T ]. Then, from "nodal cycle" 0, the goal is to calculate "nodal cycle" 1 defined over [Θ 1 , Θ 1 + ∆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.

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.

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 λ j µ j=1 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 ϑ (Θ m + τ m ) defined ∀τ m ∈ [0, ∆T ] and at Θ m . This quantity should be transformed intoθ (Θ m + τ m ), such that A linear transformation rule can be assumed which is of the form where a ϑ and b ϑ are parameters obtained from the periodicity condition of eq. (33). The initialisation ϑ 0 (Θ m+1 + τ m+1 ) for the "nodal cycle" m + 1 is thus based on the transformed quantity computed for the "nodal cycle" m as 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 Θ m + ∆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.

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 (Θ 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 allp − 1 cycles are known (see fig. 3), D (Θ m = θ 0 ) can be used to calculate where ∆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. ∆D ini and ∆D f in , 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 where ∆D ini and ∆D f in 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.

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 ∀ u * ∈ U 0 as with 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 The test field in eq. (38a) when represented in a separated form, i.e. u * = λ u * ū + λ uū * withū * ∈ U 0 and λ 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 λ u so that ∀λ u * , The simplistic assumption from eq. (41) is λ u = λ p = λ (see Relun et al., 2011, for instance), which results in the corrective terms of stress ∆σ ′ i+1 be written as The spatial problem thereafter attempts to findū ∈ U 0 such that This allows to define operators E and C (see Appendix A for details) in the finite element space, that linksε p with the other space functions asε The quantities of interest depending on the plastic strain thereby become The corrective terms ∆σ i+1 and ∆ε i+1 are obtained from the weak form of eq. (38b) along with eq. (39b), i.e.
Knowing ∆ε R i+1 from the local stage, the corrective terms of strain and the corresponding stress tensor can be estimated using the operators E and C as 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 µ 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 The updates of the time functions λ j µ j=1 are evaluated from the minimisation of the constitutive relation error e 2 CE defined for the given "nodal cycle" as the norm of the search direction equation given by eq. (22) as which leads to 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 spacetime pair will be added. The quality of the ROM is evaluated through an error saturation indicator, which is given by, The approximation for this LATIN iteration is considered acceptable if the indicator is larger than a particular predefined value ζ tol . If the indicator is lower than ζ 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 where λ µ+1 (t) ,ε p µ+1 x is the new space-time pair. The search direction equation is now written as 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ε p µ+1 , the temporal function λ µ+1 is calculated through the minimisation of the constitutive relation error e 2 CE λ µ+1 = arg min 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 where ∆ε p is obtained from eq. (53) and ∆σ ′ = λσ withσ = λCε p . The kinematic admissibility condition in rate form is given by which leads to with · = [Θ m+1 ,Θ m+1 +∆T ] · dt and σ * = λσ. Introducingε such that with W −1 = H σ λ 2 + λ λ C −1 andδ = ∆ λ leads to Ωε :σ * dΩ = 0, ∀σ * ∈ S 0 .
The static admissibility to zero ofσ is written as which may be reformulated as using eq. (60) and introducing displacementū ∈ U 0 such that ε(ū) =ε. The pseudo-displacementū along with the associated strainε is obtained by solving the finite element problem as usual. Thereby the space functionε p µ+1 is estimated as, 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 and It must be noted that the calculation of H β and H σ 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 with 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.

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 ζ < 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.

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 2πt ∆T , with time period ∆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.   Lemaitre and Desmorat, 2005) 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.   whereε p ms is the mono-scale reference solution andε p ts 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 cal ts /t cal ms , where t cal ms is the time needed for the mono-scale computation and t cal ts 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. 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σy) 2 2E . Parameters k D = 0.128 MPa −n d 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 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. 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.  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 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.

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×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. 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.

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.  The size of the time element [Θ m , Θ m+1 ] i.e. the value of Θ m+1 is estimated using the maximum damage value at the end of "nodal cycle" m denoted by D max m , according to the scheme presented in table 3. Using a "training stage" comprising of only 10 cycles, a total of 5 × 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. 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.

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 φ = 4 mm as illustrated in fig. 14. Uniformly distributed sinusoidal displacements of amplitude U 0 = 0.0057 mm and time period ∆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". 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. 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.

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. 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 2πt ∆T with time period ∆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. 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.

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.

AppendixA. 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. Ω Cε : ε ū * dΩ = Ω Cε p : ε ū * dΩ, ∀ū * ∈ U 0 . (A.1) Classical finite element scheme is used for the spatial discretisation. The left hand side of eq. (A.1) for one elemental domain Ω el can be written as (A.2) and the corresponding right hand side reads where {ū el } and ε p el denote the nodal unknowns and the unknowns at all Gauss points, for one element respectively. The matrix B comprises of the derivatives of the shape functions. Assembling the elemental equations leads to the following global problem The derivatives of the shape functions contained in the total assembled matrix are used to express eq. (A.6) as with E being the finite element operator that relates the space functions corresponding to plastic strain ε p tot and total strain {ε tot }.
Using the strain partition relationship, ε e tot = {ε tot } − ε p tot (A.8) the space functions related to stress can be expressed using eq. (A.7) as 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 ε p tot and stress {σ tot }.