A Model Reduction Technique in Space and Time for Fatigue Simulation

The simulation of mechanical responses of structures subjected to cyclic loadings for a large number of cycles remains a challenge. The goal herein is to develop an innovative computational scheme for fatigue computations involving nonlinear mechanical behaviour of materials, described by internal variables. The focus is on the Large Time Increment (LATIN) method coupled with a model reduction technique, the Proper Generalized Decomposition (PGD). Moreover, a multi-time scale approach is proposed for the simulation of fatigue involving large number of cycles. The quantities of interest are calculated only at particular pre-defined cycles called the “nodal cycles” and a suitable interpolation is used to estimate their evolution at the intermediate cycles. The proposed framework is exemplified for a structure subjected to cyclic loading, where damage is considered to be isotropic and micro-defect closure effects are taken into account. The combination of these techniques reduce the numerical cost drastically and allows to create virtual S-N curves for large number of cycles. M. Bhattacharyya (B) · A. Fau · U. Nackenhorst IBNM, Leibniz Universität Hannover, Appelstraße 9a, 30167 Hannover, Germany e-mail: mainak.bhattacharyya@ibnm.uni-hannover.de A. Fau e-mail: amelie.fau@ibnm.uni-hannover.de U. Nackenhorst e-mail: nackenhorst@ibnm.uni-hannover.de D. Néron · P. Ladevèze LMT, ENS Paris-Saclay, CNRS, Université Paris-Saclay, 61 avenue du Président Wilson, 94235 Cachan, France e-mail: neron@lmt.ens-cachan.fr P. Ladevèze e-mail: ladeveze@lmt.ens-cachan.fr


Introduction
The phenomenon of fatigue has been of great importance in the design of mechanical components and civil structures. Fatigue in the most generic case can be defined as the change in properties of a structure subjected to repetitive loading. The focus here is on mechanical fatigue involving fluctuations in externally applied load with respect to time. The nature of the loading however can be completely periodic, nonperiodic, or random. Whatever the nature of fatigue maybe, the failure of the structure is composed of three stages [1]: 1. nucleation and growth of micro-voids into a macrocrack, 2. stable propagation of the macro-crack, 3. unstable crack propagation that leads to complete failure.
There are three main engineering approaches to ensure the safety of a structure [1]. The "safe-life" approach demands that the structure remains safe under a prescribed load for a certain number of cycles. The "fail-safe" approach requires that a structure is able to withstand damage without a catastrophic failure of the whole structure. The "damage tolerance" approach investigates the ability of a structure to survive with damage before reparation. Whatever the engineering goal, researchers have to represent the complex phenomena of cyclic fatigue or random fatigue with pertinent and flexible models.

Different Modelling Approaches for Fatigue Analysis
Historical approaches for fatigue analysis are based on the exploitation of some empirical curves introduced by Wöhler [2]. The stress range is plotted against the experimentally observed number of surviving cycles giving rise to S-N curves. This may be a straightforward tool for "safe-life" analysis if a perfect periodic loading is considered with the same mean stress as the one of the S-N curve. From the asymptotic behaviour of the curve, an endurance limit, below which the material will never fail, may be defined. For loading inducing appreciable plastic deformation e.g. low cycle fatigue, it is better to introduce a strain-life approach as proposed by Coffin [3] and Manson [4] for metallic materials. Similar to S-N curves that are used for stress-life approaches, ε a -N or ε-N curves are generally used for strain-life approaches [5], where ε a and ε are the total strain amplitude and total strain range respectively. The empirical method may be extended to different sophisticated cases e.g.: • taking into account the mean stress effect on fatigue as proposed by Gerber [6], Goodman [7] and Soderberg [8]; • characterising the fatigue life of a structure when subjected to blocks of cyclic stresses of different amplitudes by using a cumulative damage rule coupled with S-N curve. The accumulation damage rule may be linear as the Palmgren-Miner damage rule [9,10], or non-linear accumulation rule as proposed by Marco and Starkey [11]; • for notched members using the Neuber rule [12], which takes into account the stress concentration factor and the fatigue notch factor; • for multiaxial stresses, an effective stress based for example on von Mises criterion, Crossland criterion [13], Sines criterion [14], or Dang Van criterion [15] is used together with the Basquin relation. An effective strain or effective plastic strain may either be estimated according to the maximum shear strain theory or according to the distortional energy theory; • for random loading, by cycle counting approach the random load history is representing as a discrete number of equivalent cycles, for example using rainflow counting.
The aforementioned methods, are phenomenological in nature and deal with the usage of empirical relations that are based on experimental findings. For example, they do not offer the flexibility to investigate the order in which different levels of load are applied to the structure. To overcome this limitation, the changes in a given structure can be described by internal variables using continuum mechanics approaches. Damage, defined as an internal variable, is used to quantify the initiation phase of a macro-crack that reduces the load carrying capacity of the material [5]. Continuum damage mechanics to predict fatigue life was introduced by Lemaitre and Lesne [16] where a non-linear continuous fatigue damage model was used to describe the different phases of the deterioration process. Various modifications and developments have been done over the years for the modelling of fatigue damage to incorporate the physical phenomena as far as possible. For instance, when the load is large, the structure undergoes appreciable plastic deformation, leading to a fatigue life less than 10 5 cycles, which is referred as low-cycle or oligocyclic fatigue. On the contrary, under high-cycle fatigue (HCF), the loading is much less than the yield stress. Therefore, no macro-plastic deformation is involved. The structure may handle a very large number of cycles. The two-scale damage model [17] is an important development for modelling HCF, and represents the macro-elastic behavior while the damage is evaluated only at the micro-scale. Different micromechanical models may be proposed to represent the micro-void growth [18].
To investigate macro-crack propagation, the premise of fracture mechanics is used knowing the pre-existence of a crack within the material. The most traditional law to describe the crack growth is the Paris-Erdogan law [19]. More information on classical techniques and recent developments on fatigue modelling and simulations can be found in the review articles [20,21].
Here, a continuum damage approach is used for modelling the fatigue behaviour. This allows to consider the chronology of the different cycles or inertial effects due to high-frequency cycles. However, this approach may lead to very expensive computational cost. Developments have been done in the light of new and robust numerical techniques that reduce CPU cost. For example, this problem may be overcome by using model order reduction techniques.

Model Reduction Technique for Fatigue Computation
Model order reduction is a family of numerical strategies which has shown efficiency for many large size mechanical problems such as parametric studies or real-time computations [22]. The solution is approximated by solving the Galerkin problem within a reduced-order basis, whose dimension is much smaller than the size of the original high-dimensional model. Proper Orthogonal Decomposition (POD) as proposed e.g. in [23] is based on a training stage. From the solution of the full-order problem at some particular time instants and/or parameter values arbitrarily chosen, a reduced-order basis is built as the truncation of a singular value decomposition [24]. In Reduced-Basis approach a greedy algorithm is employed to define the most relevant calculations within the parametric space to optimise the enrichment of the reduced-order basis [25]. Using Proper Generalised Decomposition (PGD), the problem is also solved in a relevant reduced-order basis, but the basis is defined on-the-fly by a greedy algorithm [26][27][28].
As computations for continuum damage problems may face strain localisation, numerical response is highly sensible to any modification of the model. Therefore damage computational strategies using model order reduction are challenging and may be hazardous. Controlling the accuracy of POD computations with circumspection is recommended and adaptive schemes such as A Priori Hyper Reduction Method [29] or POD coupled with Newton-Krylov algorithms [30] are preferred. Delamination has been modelled by cohesive zone and PGD by Metoui et al. [31], PGD-based multi-scale computations for rate-dependent damage model have been recently proposed by El Halabi et al. [32].
The LATIN method [26] in its classical sense is employed such that an approximation of the solution on the whole time-space domain is defined at every iteration, and the gobal equilibrium is tackled as a linearised problem and the non-linear material behaviours are considered separately. Therefore, this approach offers a convenient framework to include model-order reduction techniques even for non-linear computations. It has been developed for solving plasticity and visco-plasticity problems and shown a drastic decrease of the computational cost compared to a classical approach [33][34][35]. Here an extended version is used for (visco-)plastic problems with damage that incorporates micro-defects closure effects [36].
For fatigue computation applied to a large number of cycles, advanced numerical strategies may be required to work around the computational cost due to the large time domain.

Efficient Time Schemes for Fatigue Computations
For high or very high cycle fatigue or for combined cycle fatigue in which the load is a combination of large amplitude low frequency and small amplitude high frequency loads [1], the computational cost of the time integration scheme may be extremely high. The computational cost may be decreased by benefiting from the cyclic loading to not compute explicitly the full number of cycles.
The jump cycle procedure is a very robust technique [5] which works around full blocks of cycles. After computing in details a set of cycles, they are used to establish a trend line and extrapolate the quantities of interest for further evolution. Then, the extrapolated state is used as initial condition for a further computation such that the whole time domain is spanned. A control function has been proposed such that the lengths of cycle jumps are monitored to ensure an accurate approximation [37]. This method is most suitable for quasi-linear systems. However the control function allows to also consider non-linear behaviour and cycles jumps are automatically shortened or even cancelled [38,39].
For combined cycle fatigue, time homogenisation techniques are based on a hypothesis of time scale separation between a large time scale associated with the low-frequency load and a small time scale associated with the high-frequency load. The ratio between both of them is assumed to be small enough such that the two scales can be considered as independent, and the behaviour due to the high frequency is homogenised for a time discretisation scheme consistently defined for the low-frequency load. Developed initially for plasticity and quasi-static computations [40], this method has been extended to damage and dynamic effects [41].
In the framework of LATIN method, as the whole time domain is available at each iteration, a two-time scale in a "finite element like" discretisation of the time intervals has been proposed for visco-plasticity problems [26,42]. Some cycles of interest, called "time nodes" or "nodal cycles", are computed in the classical time scheme. The evolution of the quantities of interest in the time slot between two time nodes is interpolated using shape functions as for space finite-element method.
As empirical methods are based on expensive and cumbersome experiment, the goal of this project is to provide virtual experiments based on continuum mechanics approach. To overcome the computational cost due to these models, a sophisticated model order reduction scheme is proposed which reduces drastically the cost and allows to compute virtual test even for a very large number of cycles. The proposed numerical framework is flexible concerning the damage model, the global framework and equations considered are summarised in Sect. 2. Then, the model order reduction based on the separation of variables is introduced in Sect. 3. Finally, the numerical treatment of the time problem, which allows to reduce the numerical cost, is detailed in Sect. 4.

Continuum Damage Mechanics Approach
The focus herein is on the methods that deal in quantifying the change in material properties due to fatigue loading using internal damage variable in a continuum mechanics framework.
A continuous reference structure for any quasi-static analysis can be considered in a spatial domain . For simplistic case, the evolution of the state of the structure can be considered to be isothermal within time domain [0, T ]. Such a reference structure in general is subjected to prescribed body forces f d , to traction forces F d over a part ∂ 2 of the boundary ∂ , and to prescribed displacements u d over the complementary part ∂ 1 . Thereby the compatibility of the applied forces with the internal stress is given by the equilibrium equation which in a weak form becomes the static admissibility condition, which is defined such that ∀u * ∈ U * (1) with the stress tensor σ being statically admissible, and U * being the homogeneous vector space associated to the space U of kinematically admissible field u. Also, the compatibility of the prescribed displacement with the generated strain within the structure is given by the strain-displacement relationship which in weak form becomes the kinematic admissibility condition, which is defined such that ∀σ * ∈ F * with ε being the total strain tensor which is kinematically admissible and can be split into an elastic part ε e and a plastic part ε p additively. F * is the homogeneous vector space associated to the space F of statically admissible field σ . The mechanical properties of the materials are described by a set of constitutive relations. The equations of state for elasto-(visco)plastic materials subjected to unilateral damage are obtained from a free energy function, and are given as [5] ε e i j = Here, Eq. (3a) represents the elastic state law, which because of damage is nonlinear in nature, as it is not described by a linear operator, with E and ν being the modulus of elasticity and Poisson ratio respectively. D is the isotropic damage variable and h is the closure parameter representing micro-defects closure effect, which indicates that the effect of damage is more predominant in tension than in compression. During compression, some of the micro-defects are closed, thereby the effective area is increased, the material regains some of the stiffness. This can be represented by the effective modulus of elasticity during tensionẼ + , which is given , and by the effective elastic modulus during compressionẼ − given byẼ − = E (1 − h D). The closure parameter h has values between 0 (complete stiffness recovery) and 1 (no stiffness recovery). Equations (3b) and (3c) give the relationships between the internal variables α i j and r describing kinematic and isotropic hardening of the material and their corresponding thermodynamic forces β i j and R, by means of material parameter C for kinematic hardening and through a function g for isotropic hardening. Equation (3d) defines the strain energy release rate Y , which is the thermodynamic force corresponding to damage and is non-linear with respect to the stress tensor and the damage variable. The evolution equations obtained from potential functions are written aṡ This is the normality rule for standard materials. Equation (4a) gives the evolution of the plastic strain with respect to time, withλ p being the plastic multiplier which is measured from the consistency condition, and F is a potential which can be obtained from experimental findings. This function has to be convex, non-negative and should pass through the origin. ∂ F ∂σ i j is the flow vector indicating plastic flow is normal to the potential function F. This potential F, for associative plasticity is considered to be the yield function. Equations (4b) and (4c) similarly, give the evolution of the internal variables for kinematic and isotropic hardening respectively. Equation (4d) describes the evolution of damage with respect to time, whereλ D is the damage multiplier and F D is a potential which is also identified from experimental findings. All the conditions required to formulate F must also be taken into account for F D . These evolution equations, essentially, take into account the history dependency of the material.
The initiation of the macro-crack is indicated by the critical damage D c and the material is assumed to fail when D c is reached.

LATIN-Based Model Order Reduction Approach for Damage Computation
The LATIN approach tackles the set of equations on the whole time-space domain at every iteration. The equations are considered iteratively, namely the global equilibrium of the structure on one hand, the non-linear elastic law and the non-linear evolution equations on another hand. Model reduction is used in the sense of separation of time and space variables. LATIN-PGD approach allows for visco-plasticity case to define both quantities of interest which are stress and plastic strain using a unique time basis. Including damage, a new quantity of interest is added, which is the elastic strain, the total set of the solution field can be represented as s = ε p , ε e ,Ẋ,Ḋ, σ , Z, Y , where X represents the set of hardening variables (kinematic and isotropic) and Z is the conjugate variable of X. The non-linear state law could be included in the global stage and then the separation of variables tackled by using different time functions for stress and strain. Otherwise, as presented here, the non-linearity due to the state law can be tackled in the local stage, and the stress can be decoupled in two parts, one defined from the local stage and the other one which can be written as a time-space decomposition form using the same time functions as the plastic strain one. The algorithm is only briefly overviewed in this contribution. Detailed explanation of every step and details about the numerical implementation may be found in [36]. The algorithm is initialised by solving the problem considering the boundary conditions of the exact problem but assuming that the material behaviour is perfectly elastic for the whole loading conditions. Then, plastic and damage corrections are added to the elastic solution at each subsequent iteration. The set of equations is divided in two sub-groups, one comprising the global and linear equations whereas the other one comprises the local and non-linear ones. One LATIN iteration consists of two parts: • the global and linear problem is solved in the space A d which belongs to the manifold of the admissibility conditions Eqs. (1) and (2), the linear state laws Eqs. (3b) and (3c), and the non-linear state law for damage Eq. (3d); • the local and non-linear problem is solved in the space Γ which belongs to the manifold of the evolution equations Eq. (4) and the elastic state law Eq. (3a) which was not linearisable due to damage.
It can be noted that Eq. (3d) although being non-linear is tackled with the group of linear equations as a post-processing step from the knowledge of the stress tensor and the damage variable at the end of each iteration. The exact solution s ex of the problem is defined as the intersection of the two manifolds by The approximation of the solution is looked alternatively in the two manifolds until reaching convergence. From the knowledge of one step, the approximation in the following manifold is looked for by using certain linear operators called search direction operators.

Local Stage
In the local stage, evolution equations for internal variables, which are local in space and non-linear, are solved. The elastic state law, being non-linear, is also tackled in this stage. From the solution set s i ∈ A d at LATIN iteration i, the approximation s i+1/2 ∈ Γ is estimated such that the local search directions are satisfied Here, B + is the direction of ascent. Following [26], the search direction is considered to be vertical such that B + −1 = 0.
The solution of the search direction equation (Eq. 6) along with the evolution equations (Eq. 4) and the non-linear elastic law (Eq. 3a) constituteŝ i+1/2 . From the approximation at the local stageŝ i+1/2 , the solution set s i+1 is estimated in the global stage.

Global Stage Including Model Order Reduction
In the global stage, the solution set s i+1 ∈ A d satisfies the state laws, the admissibility conditions and the descent search directions ⎡ ⎢ ⎣ε The operator H − belongs to the tangent space associated with the solution setŝ i+1/2 in the manifold Γ , and C is the undamaged Hooke tensor. Considering the damage variable is not updated in the linear stage, the search direction operator b − is defined as zero.
The first step being the calculation of the hardening variables, the state equations are combined in the form where is a linear operator containing the sate law parameters. The search direction equation for hardening variables Eq. (8a) combined with the state equation (Eq. 9) can be written as with H Z being the decoupled part of H − that relates the internal variables with the corresponding associated variables. The hardening variables are then obtained by solving Eq. (8a) in time at each Gauss point. The difficulty to calculate the stresses and strains, compared to former works with the LATIN method, is that the elastic state law Eq. (3a) is non-linear due to the presence of damage, leading to solve a non-linear problem at the global stage. This point is particularly tricky as it prevents the introduction of a model reduction strategy at this stage. The idea proposed herein is to transform this non-linear problem into separate linear equations by decomposing stress and total strain into two parts depending on plastic deformation and damage respectively.
The quantities of interest at this point σ i+1 , ε e i+1 andε p i+1 are represented in a corrective form at iteration i + 1 as The stress and total strain corrections in the global stage at iteration i + 1 are separated into parts depending on plastic deformation ( σ i+1 , ε i+1 ) and on damage ( σ i+1 , ε i+1 ), From these separations and the search direction equation along with the additive strain decomposition relation, it can be established that where ε R i+1 can be interpreted as a residual strain obtained from non-linear state law at iteration i + 1. σ i+1 and ε i+1 are thereby obtained from the equilibrium equation, directly.
On the other hand, if only the plastic part is considered, the search direction can be re-written as ε with¯ i+1 is a plastic corrective term from the local stage and H σ represents the decoupled part of H − that relates stress to plastic strain rate.
The correction terms linked with the plastic behaviour σ i+1 and ε i+1 are then written in a separable form using Proper Generalised Decomposition.

Separation of Variables
The Proper Generalised Decomposition (PGD) is a flexible model order reduction technique, which is not based on a training stage. As at every LATIN iteration, the quantities of interest are approximated on the whole space-time domain by a linear form of the mechanical equilibrium, the usage of PGD coupled with LATIN is convenient. As any function dependent on several independent variables can be approximated as an infinite sum of products of one-variable functions [26,43], PGD looks for an approximation of any quantity of interest as a finite sum of products of low-dimensional functions defined by a greedy algorithm. Therefore, this approximation includes an error due to the truncation of the series of separable forms. Here the plastic strain and stress part due to plastic deformation dependent on space and time variables are approximated aṡ where μ is the number of pairs involved in the decomposition, and C is a linear operator which relates the space functions of stress and plastic strain.

Updating Stage
The greedy algorithm is such that after defining on-the-fly a first pair of space and time functions, at every iteration a first decomposition is looked using the previously defined space functions and updating the time functions. This step is equivalent to a Proper Orthogonal Decomposition on the current space basis. Considering μ spacetime modes have been generated to approximate the stress and plastic strain rate at iteration i, the corrections of stress and the plastic strain rate at iteration i + 1 are given as The updates of the time functions are calculated by minimising a mechanical residual which is defined as the norm of the search direction operator, i.e.
Then, if the improvement of the approximation is not efficient enough, a new pair is added to the decomposition.

Enrichment of Space-Time Bases
The objective of the enrichment phase is to add a new space-time pair. The corrections of stress and plastic strain for this case are written as with the intention of calculating the separable quantitiesλ μ+1 andε p μ+1 . A hybrid strategy is used, the space functionε p μ+1 is calculated from a Galerkin formulation, by using the kinematic admissibility condition (Eq. 2) such that ∀σ * which is statically admissible to zero [0,T ] × ε : σ * d dt = 0 (19) and the static admissibility condition (Eq. 1) such that ∀u * which is kinematically admissible to zero Subsequently the time function λ μ+1 is solved similarly as in the update stage by minimising a mechanical residual This fixed point iteration between space and time problems converges quickly. Once the stress tensor is known at iteration i + 1, the strain energy release rate for damage is calculated from Eq. (3d).

Numerical Example of a Plate Under Cyclic Loading
The proposed usage of LATIN-PGD for damage problem is exemplified for a twodimensional problem as depicted in Fig. 1. Material considered is a Cr-Mo steel at 580 • C [5]. The variation in material property is represented by the yield stress, with σ y = 80 MPa at x = 0, ∀y ∈ [0, W ] and σ y = 85 MPa at x = L , ∀y ∈ [0, W ] and with a linear variation along the length of the structure.
The distribution of the damage variable D at t = T is depicted in Fig. 2. It shows localisation, with maximum near (x, y) = (L , W ) and minimum near (x, y) = (L , 0).
For fatigue computation, calculation of the time dependent quantities is expensive. Therefore a peculiar effort is done to reduce the cost of the time-computation.

A Two-Time Scale Approach
In the case of a large number of cycles, the previous strategy can be enhanced by adding a multi-time scale feature. To describe the time dependent quantities defined on the whole time domain [0, T ], two time scales are introduced, • a long time discretisation θ defined on the interval [0, T ] that represents the slow evolution along the cycles, • a short time discretisation τ describing the rapid evolution within a cycle.
The idea is to introduce a finite element like description of the temporal quantities which are calculated only at certain chosen cycles called the "nodal cycles" (Fig. 3). For any time element θ m , θ m+1 once the nodal cycles m and m + 1 are known a linear one-dimensional interpolation formula can be used [42]. If χ represents any temporal quantity over time element θ m , θ m+1 , then where χ m (τ ) and χ m+1 (τ ) are defined ∀τ ∈ [θ m , θ m + T ] and ∀τ ∈ θ m+1 , θ m+1 + T respectively. The first few cycles are computed classically. After that the nodal cycles are calculated progressively. The last classically computed cycle defined over [θ 0 , θ 0 + T ] becomes the nodal cycle 0 and thereby the idea is to calculate nodal cycle 1, defined over [θ 1 , θ 1 + T ]. Thereafter knowing nodal cycle 1, nodal cycle 2 is calculated. This computation is continued till the last nodal cycle is calculated. The initialisation of the quantities of interest at nodal cycle m, knowing them at nodal cycle m − 1 depends on the quantities that are being initialised. The quantities that are cyclic, namely the stress, elastic strain as well as the kinematic variables, are duplicated from the nodal cycle m − 1 and a transformation is considered such that they are periodic, preserving the continuity from cycle m − 1. The time functions representing the plastic strain λ j μ j=1 are also duplicated in a similar manner. The quantities that are non-cyclic and non-decreasing, namely damage and isotropic variables, are initialised as constant over nodal cycle m, with the magnitude being that obtained at θ m−1 + T . The strain energy release rate for damage is calculated from the damage and the stress tensor. Thereafter the solution field over the nodal cycle m is calculated iteratively using the two-step algorithm till a convergence is obtained.

Local Stage
In the local stage, all the quantities of interest except damage do not need any time integration and can be calculated directly. The only problem in this stage is while integratingḊ to obtain the damage variable. To integrate Eq. (4d) over the nodal cycle m the initial condition at θ m needs to be known. Considering a general first order ODE [26] dχ dt defined over the complete time domain, with κ and υ being time-dependent known quantities. The idea is to calculate χ (θ m ) from χ (θ m−1 ). The time element θ m−1 , θ m is discretised into certain instances k such that k = θ m−1 + k T , with k = 0, 1, 2, · · · , p − 1, where p is the number of cycles in the time element θ m−1 , θ m . This provides 0 = θ m−1 and p−1 = θ m . Knowing χ θ k , Eq. (24) can be solved to obtain χ θ k+1 as whereχ represents the solution of the ODE with zero initial condition, and represents the "resolvent" operator [26]. The challenge henceforth is to calculateχ and with minimum numerical cost. The easiest way is to calculate the quantities only at the nodal cycles m − 1 and m and then use linear interpolations to obtainχ k+1 , k and k+1 , k , i.e.
with ν being the linear shape functions defined as These interpolated values can be used to rewrite Eq. (25) as

Global Stage
The spatial modes, calculated for the initial cycles, are reused to compute the time functions of the PGD modes over the nodal cycle m. The initialisation of the time functions {λ j } μ j=1 is such that continuity is maintained with respect to the nodal cycle m − 1. Thereafter corrections of the time functions { λ j } μ j=1 are computed by solving Eq. (17) using zero initial conditions. The next concern are the kinematic variables, which also being cyclic can be treated in a similar way. Equation (10) has to be solved for the kinematic variables over the nodal cycle m. However, an exact measurement of the initial conditions using the "resolvent" technique, being numerically expensive, is not necessary. The initialisation of the kinematic variables has been done to maintain continuity with respect to the nodal cycle m − 1. Thereby the quantities are calculated in terms of corrections by solving the ODE with zero initial condition.
The isotropic variable, however, being non-cyclic needs an accurate measurement of the initial condition for time integration of Eq. (10) over the nodal cycle m. This first order ODE is solved using the "resolvent" technique previously described.
This approach provides a drastic reduction in cost compared to classical LATIN technique.

A Numerical Example of a Bar Under Fatigue Loading to Build Virtual ε a -N Curves
The one-dimensional numerical example considered here is a bar under traction to build virtual ε a -N curves. The material considered is Cr-Mo steel at 25 • C with  Experimental ε a -N are generally obtained when a particular specimen is loaded under a given ε a and the number of cycles needed for the specimen to rupture is measured. This experiment is repeated for several values of ε a to obtain different values of N . For the numerical tests described here, a critical damage value of 0.2 is considered as a failure point. Similar to the physical experiments, several numerical tests are conducted by varying ε a to obtain different values of N needed by the structure to reach the critical damage level. Some virtual ε a -N curves are depicted in Figs. 4 and 5. The time domain for each numerical test is discretised uniformly, with a lesser number of cycles per time element for larger ε a . The range of the number of cycles per element is between 10 and 200. The computation of each curve requires approximately 1 hr.
The influence of yield stress σ y is depicted in Fig. 4. It is observed that, with the increase in σ y , the structure becomes less susceptible to damage for a given strain amplitude. As the damage threshold considered in the model is directly proportional to the yield stress, more number of cycles are needed to reach the critical damage. The influence of elastic modulus E is also shown in Fig. 4. It is witnessed that, with the increase in E, the stress in the structure increases for a given strain amplitude, resulting in increased susceptibility to damage. Also, the damage threshold, considered in the model is inversely proportional to the modulus of elasticity, resulting in a higher damage for a higher value of E at a given ε a . The aforementioned reasons culminate in the decrease in the number of cycles to failure with increase in E for a given strain amplitude. It is also noticed that this decrease is profound at lower strain amplitudes.
The next sets of tests investigate the influence of initial damage and mean strain on the ε a -N curves, as shown in Fig. 5. The first set of tests investigates the influence of the presence of initial damage. It can be observed that, with the increase in magnitude of the initial damage, less number of cycles are needed to reach the critical damage for a given strain rate. The subsequent set of tests consists of different total mean strains ε m . For a positive mean strain, the stress is more in the tensile part than in the compressive part, thereby a higher damage is obtained resulting into a lower N , for a given ε a , than compared to the zero mean strain case. For a negative mean strain, the stress is more in the compressive region, thereby damage evolution is less compared to the zero mean stress case, resulting into a higher N for a given ε a .
To evaluate the accuracy and the efficiency of the two-time scale approach a case with U 0 = 1.40 mm for 2000 cycles is investigated. A mono-scale LATIN-PGD computation obtained in a CPU time of 19 hrs is considered as a reference. The performance of the two-time scale approach is analysed using various uniform time discretisations.
The evolution of damage for different sizes of time elements is plotted in Fig. 6. For decreasing size of time elements the evaluation converges to the reference solution. The relative error of the two-time scale approach with respect to the reference solution is defined as where D ms and D ts are the damage variables computed using the mono-scale computation and the two-scale methods respectively. The evolution of this error and the computational time with respect to the size of time elements is depicted in Fig. 6 too. The multi-scale computational times are represented as percentages of the monoscale computational time in Fig. 6. Using only one time element containing 2000 cycles, the computational cost has been decreased to 0.06% of the reference cost, but the error is 16% of the reference. By increasing the number of time elements, the computational cost increases but the error decreases. For 40 elements of 50 cycles, the error is only 0.36% for a computational cost 0.3% of the reference solution. The calculation time is drastically less compared to a mono-scale LATIN-PGD approach for acceptable accuracy.

Conclusion
An innovative numerical approach has been here presented for the computation of fatigue damage. A non-incremental technique has been used as numerical framework and the computational cost has been reduced by the usage of PGD that converts the actual high dimensional problem into much lower dimension. A multi-scale approach in time domain has been suggested to extend the LATIN-PGD method for the simulation of fatigue damage behaviour involving large number of cycles. By that, the increased numerical efficiency has been demonstrated which paves a way for physically based high cycle fatigue simulations.