Keywords: LATIN method, damage mechanics, (visco)plasticity, metal fatigue, model order reduction, PGD, multi-temporal scale méthode LATIN, mécanique de l'endommagement, (visco)plasticité, fatigue des métaux, réduction

The motivation of the research project is to predict the life time of mechanical components that are subjected to cyclic fatigue phenomena. The idea herein is to develop an innovative numerical scheme to predict failure of structures under such loading. The model is based on classical continuum damage mechanics introducing internal variables which describe the damage evolution. The challenge lies in the treatment of large number of load cycles for the life time prediction, particularly the residual life time for existing structures.

Traditional approaches for fatigue analysis are based on phenomenological methods and deal with the usage of empirical relations. Such methods consider simplistic approximations and are unable to take into account complex geometries, and complicated loadings which occur in real-life engineering problems. A thermodynamically consistent continuum-based approach is therefore used for modelling the fatigue behaviour. This allows to consider complicated geometries and loads quite efficiently and the deterioration of the material properties due to fatigue can be quantified using internal variables. However, this approach can be computationally expensive and hence sophisticated numerical frameworks should be used.

The numerical strategy used in this project is different when compared to regular time incremental schemes used for solving elasto-(visco)plastic-damage problems in continuum framework. This numerical strategy is called Large Time Increment (LATIN) method, which is a non-incremental method and builds the solution iteratively for the complete space-time domain. An important feature of the LATIN method is to incorporate an on-the-fly model reduction strategy to reduce drastically the numerical cost. Proper generalised decomposition (PGD), being a priori a model reduction strategy, separates the quantities of interest with respect to space and time, and computes iteratively the spatial and temporal approximations. LATIN-PGD framework has been effectively used over the years to solve elasto-(visco)plastic problems. Herein, the first effort is to solve continuum damage problems using LATIN-PGD techniques.

Although, usage of PGD reduces the numerical cost, the benefit is not enough to solve problems involving large number of load cycles and computational time can be severely high, making simulations of fatigue problems infeasible. This can be overcome by using a multi-time scale approach, that takes into account the rapid evolution of the quantities of interest within a load cycle and their slow evolution along the load cycles. A finite element like description with respect to time is proposed, where the whole time domain is discretised into time elements, and only the nodal cycles, which form the boundary of the time elements, are calculated using LATIN-PGD technique. Thereby, classical shape functions are used to interpolate within the time element.

This two-scale LATIN-PGD strategy enables the reduction of the computational cost remarkably, and can be used to simulate damage evolution in a structure under fatigue loading for a very large number of cycles.

Zusammenfassung

Das Ziel des Forschungsprojektes ist es, die Lebensdauer mechanischer Tragwerke unter zyklischen Ermüdungsphänomenen vorherzusagen. Die zugrunde liegende Idee ist, ein innovatives numerisches Schema zu entwickeln, um das Versagen der Strukturen unter bestimmten Belastungen zu berechnen. Für das Kontinuumsmodell werden innere Variablen eingeführt, um die Entwicklung der Materialermüdung zu beschreiben. Die Herausforderung liegt dabei in der numerischen Behandlung der Belastungszyklen für die Voraussage der Lebensdauer, insbesondere für die Voraussage der Restlebensdauer der Strukturen.

Traditionelle Ansätze zur Ermüdungsanalyse basieren auf phänomenologischen Methoden welche aus emperischen Beobachtungen abgeleitet werden. Solche Verfahren berücksichtigen einfache Annäherungen und können komplexe Geometrien und komplizierte Belastungen, die in realen technischen Problemen auftreten, nicht abbilden. Zur Modellierung des Ermüdungsverhaltens wird daher ein thermodynamisch konsistenter kontinuumsbasierter Ansatz verwendet. Dieser ermöglicht es, komplizierte Geometrien und Lasten recht effizient zu betrachten, und die Verschlechterung der Materialeigenschaften aufgrund von Ermüdung kann unter Verwendung interner Variablen quantifiziert werden. Dieser Ansatz kann jedoch numerisch teuer sein und daher sollten effiziente numerische Techniken entwickelt werden.

In diesem Projekt wird eine innovative numerische Strategie basierend auf einer nicht-zeitinkrementellen Lösung von elasto-(visko)plastischen Schädigungsproblemen in der Kontinuumsmechanik vorgestellt. Diese numerische Strategie basiert auf der sogenannten LATIN Methode (Large Time Increment method). Bei dieser nicht-inkrementellen Methode wird iterativ im vollständigen Zeit-Raum Bereich eine Lösung gesucht. Eine wichtige Eigenschaft der LATIN Methode ist eine ad hoc Modellsreduktionsstrategie um die Rechnenkosten zu reduzieren. Die Proper Generalised Decomposition (PGD) ist eine a priori Modellreduktionsstrategie, die in Raum und Zeit die Variablen in zwei unabhängige Teile auftrennt, einen zeitabhängigen und einen raumabhängigen. Iterativ wird eine Approximation dieser beiden Anteile gesucht. Die Effektivität des LATIN-PGD Verfahrens ist seit Jahren für die Lösung von elasto-(visko)plastischen Problemen nachgewiesen. Der erste Teil dieses Projektes zielt auf die Erweiterung dieser Methode für Schädigungsmodelle ab.

Obwohl die Anwendung des LATIN-PGD Verfahrens die numerischen Kosten reduziert, ist die Ersparnis nicht groß genug, um Probleme mit vielen Belastungszyklen lösen zu können. Dieses Problem soll mit einem Mehrskalenansatz in der Zeit gelöst werden. Damit können schnelle Entwicklung der Variablen während eines Belastungszyklus erfasst und ihre langsame Evolution während vielen Zyklen abgebildet werden. Der Zeitbereich wird dabei mit diskontinuierlichen Finiten Elementen diskretisiert, während die Knotenzyklen, die die Grenzen der Zeitelemente bilden, mit der LATIN-PGD Technik beschrieben werden. Dafür werden klassische Ansatzsfunktionen benutzt, um die Variablen innerhalb der Zeitelemente zu interpolieren.

Mit dieser innovativen Strategie werden die Rechnenkosten signifikant reduziert und die Simulation von Ermüdungsbelastung unter großen Lastzyklenzahlen ermöglicht.

Introduction

Failure of structures have been a major consideration among engineers, over the decades. Failure in a general sense can be defined as rupture or excessive deformation. Among many different loading phenomena that induce failure of structures, such as static loading or monotonic loading, one important phenomenon is cyclic loading or fatigue loading, where the load is lower than the static strength but causes failure due to repetition. Fatigue forms one of the most important failure mechanisms for civil, mechanical, aviation and automotive components. It is believed that 50-90 % of all the mechanical failures especially for metallic structures occur due to fatigue (see [START_REF] Sobczyk | Random Fatigue: From Data to Theory[END_REF]. The failure mechanism for structures subjected to cyclic fatigue is complicated in a sense that failure does not happen instantaneously, but material is degraded over a period of time and this degradation is often not apparent and in many cases structures may fail without any prior warning.

Although, the experimental studies of fatigue had been under progress for a long time, these experimental investigations, in general are very time consuming and restricted to specimens suitable for laboratory conditions. Thereby, a lot of attentions were focused on numerical studies of fatigue later on. The initial numerical studies of fatigue made in the 20th century, were empirical in nature, and could only simulate very simple cases. In the later part of the century, thermodynamically consistent models were utilised to describe the phenomenological aspects of fatigue. This idea gave a better description of the physical behaviour of the problem, which also paved the way for efficient time integration schemes to simulate large number of cycles. A major boost in the numerical modelling of fatigue process was in the late 20th century, when continuum damage mechanics (CDM), a branch of continuum mechanics, was used to model fatigue behaviour. This branch of continuum mechanics, essentially quantifies the degradation of material properties, in terms of a separate internal variable, called the damage variable. A lot of developments have been proposed from the first usage of continuum damage mechanics in [START_REF] Chaboche | A non-linear continuous fatigue damage model[END_REF] to predict fatigue, both in the aspects of modelling different types of fatigue processes, and also on different improved numerical schemes to handle large number of cycles.

The goal of this thesis is to provide a novel simulation approach to compute fatigue damage using continuum damage mechanics. The challenge in this case is obviously to handle the temporal part of the problem, which involves a large number of cycles (e.g. 10 5 cycles). Therefore, instead of using a traditional time incremental scheme, the large time increment (LATIN) method is used. The advantage of LATIN method is that it separates the difficulties into a part where only the linear global admissibilities are solved and another where the local non-linear constitutive relations can be solved. The solution can be built over the whole space-time domain, and is improved successively, with each iteration between the linear and local stages. This kind of description allows the user to introduce a model reduction strategy in space and time which not only allows reduction in numerical expenditure, but also incorporates separate treatment of the temporal part, which is necessary to tackle large number of cycles.

The invention and popularity of model order reduction (MOR) or model reduction (MR) techniques stem from the fact that simulation of full physical models consisting of thousands of degrees of freedom with structural evolutions spanning for a large time interval, can be extremely expensive. In simple words, the idea of MOR is to reduce a complete high fidelity model into several low-fidelity models which are relatively less expensive to solve. One of the most widely used model reduction methods is to find through a training or learning phase, the full solution of the problem at hand and then perform singular value 1 decomposition (SVD) of the full solution to extract the relevant modes that can be used to approximate it as much as possible. The learning phase can sometimes be computed for a certain short interval of time and the spatial basis functions that are extracted are reused to approximate the solution over a larger time interval by only recomputing the temporal basis. Otherwise, if the learning phase is for the whole time domain, the space-time bases can be reused to solve a similar problem but with certain changes, for instance small changes in boundary conditions or material parameters. This kind of model reduction method is a posteriori technique that relies on the solution of the full model and is popularly referred to as proper orthogonal decomposition (POD), Karhunen Loève expansion (KLE) or principal component analysis (PCA). However it is desirable to build an a priori method which does not rely on the complete solution of the full model. The proper generalised decomposition (PGD) is one of such techniques, which does not require the knowledge of the full model, but build, "on the fly" separable approximation of the problem. The governing equations are separated into equations depending on the number of independent variables to be considered (in most classical cases space and time are the two independent variables), and then these separated equations are solved independently using fixed point method. If the approximation is not satisfactory, a new basis vector can be added spontaneously. This method is very helpful and flexible as no training phase is required.

PGD was first introduced in the LATIN framework by Ladevèze (see [START_REF] Ladevèze | The large time increment method for the analyse of structures with nonlinear constitutive relation described by internal variables[END_REF]Ladevèze, , 1985a,b) ,b) as the "radial loading approximation" to solve the global equilibrium. PGD basically separates the quantities of interest into space and time which leads to sovereign time and space problems that can then be solved using a staggered algorithm. There have been many versions of LATIN method over the years, to solve (visco)plasticity, viscoelasticity, contact problems, however all of them have the two key features: separation of difficulties, and separation of variables using PGD.

Therefore the innovation of this research is divided into two parts:

-The first part is to incorporate CDM in LATIN framework. Although there have been certain strategies proposed in the past to solve CDM using LATIN method, a unique way to include CDM in LATIN-PGD method is suggested in this thesis.

-The next part is to introduce an innovative multi-temporal scheme so that a large number of cycles can be computed. The model reduction strategy that is used in the LATIN framework enables the temporal quantities to be treated using this time scheme. Similar strategy has been used in the past for cyclic visco-plasticity. This strategy incorporates a massive reduction in numerical expenditure.

The structure of the thesis is described in the following paragraphs.

Chapter 1 introduces the definition and classification of fatigue. A brief history of fatigue over the past two centuries is given and different methods of fatigue calculations from traditional till current developments are also detailed.

In the current work, the focus is on continuum damage mechanics (CDM), developed out of classical continuum mechanics, where the deterioration of material properties is represented in terms of internal variables in a thermodynamically consistent framework. Chapter 2 introduces the fundamental ideas of continuum mechanics, (visco)plasticity, and damage. It also highlights a few methods that have been developed over the years to incorporate CDM in fatigue problems.

The most challenging part in the solution of thermodynamic based mechanical problem is the numerical expense that arises while solving a full model having a large spatial degrees of freedom for the complete temporal domain. For fatigue simulation this issue is accentuated due to the enhancement of the temporal domain which involves large number of cycles. To solve this issue, model reduction techniques may be used which can reduce the numerical cost drastically by separating the spatial and temporal part and tackling them independently. Chapter 3 gives a brief overview of different model reduction techniques that are available and their numerical importance. Chapter 3 also gives a brief introduction of large time increment (LATIN) method, a non-incremental technique which is used as a solution framework in this research.

The usage of LATIN has been extensive over the years for the solution of non-linear problems such as plasticity or visco-plasticity, however a novel approach has been adopted in this research to solve CDM problems using LATIN method. The numerical modification in the LATIN method to include CDM is portrayed in chapter 4 along with certain academic examples, depicting the numerical robustness of the approach.

The simulation of large number of cycles for fatigue problems is an obvious numerical challenge and has been addressed within the framework of LATIN technique in chapter 5. This approach uses the idea of two separate time discretisation schemes, one within a cycle and the other along the cycles. Several numerical examples are presented to depict the efficiency of the two-time scale formulations.

Chapter 6 finally summarises the important discoveries in the research and provides a foresight on the future scope and development.

Chapter 1

Fatigue: an overview Fatigue simulation being the goal of this research, a classification of fatigue, and established approaches are presented here. It has to be noted that no novel approaches are presented in this chapter. The only purpose here is to provide the readers with a brief but comprehensive overview of fatigue. The chapter begins with a concise historical perspective of fatigue, leading on to different forms of fatigue and their engineering importance. Thereafter, the mechanism of fatigue and the classical numerical approaches towards it are reviewed.

The term "fatigue" originating from the Latin word fatigare, meaning "to tire", has been used extensively in engineering to describe the degradation of materials subjected to fluctuating loads (see [START_REF] Suresh | Fatigue of materials[END_REF]. Fatigue in a generic mechanical sense can be defined as the "change in properties" of materials under the influence of time-varying loading. This "change in properties" leads to the detoriation of the materials which ultimately culminates into the inability of the structure to sustain the intended loading. The most noteworthy feature of fatigue is that failure at macro-scale does not occur immediately, but happens after a certain number of load fluctuations has taken place (see [START_REF] Cui | A state-of-the-art review on fatigue life prediction methods for metal structures[END_REF].

History of fatigue

The first study of metal fatigue was done by a German mining engineer named Albert around 1829 who studied fatigue effect of conveyor chains (see [START_REF] Albert | Über Treibseile am Harz[END_REF]. The most noticeable activity in record after Albert was from a British railway engineer named Rankine. Although he became more famous in the field of thermodynamics, around 1842 he discussed the importance of stress concentration in the fatigue strength of railway axles (see [START_REF] Rankine | On the causes of the unexpected breakage of the journals of railway axles, and on the means of preventing such accidents by observing the law of continuity in their construction[END_REF]. The first idea of "safe life concept" in perspective of railway coaches was given by a French engineer named Morin (see [START_REF] Morin | Lec ¸ons de mécanique pratique-résistance des matériaux[END_REF]. The credit of coining the term "fatigue" to describe cracking of metals under repeated loading goes to a British researcher named Braithwaite (see [START_REF] Braithwaite | On the fatigue and consequent fracture of metals[END_REF]. The most important name in that era however was of a German engineer named Wöhler, who conducted systematic experiments on fatigue failure of railway axles and found that the fatigue strength is much below the static strength (see [START_REF] Wöhler | Versuche zur Ermittlung der auf die Eisenbahnwagenachsen einwirkenden Kräfte und die Widerstandsfähigkeit der Wagen-Achsen[END_REF][START_REF] Wöhler | Über die Festigkeitsversuche mit Eisen und Stahl[END_REF]. His self designed machines were able to perform rotating-bending as well as combined bending and torsion fatigue. His experimental results were plotted first by his successor, another German engineer named Spangenberg (see [START_REF] Spangenberg | Über die Festigkeitsversuche mit Eisen und Stahl[END_REF], which came to be known as S-N curves or "Wöhler curves" since 1936. Gerber and Goodman studied the effect of mean stress in the late 19th century (see [START_REF] Gerber | Bestimmung der zulässigen Spannungen in Eisen-constructionen[END_REF][START_REF] Goodman | Mechanics Applied to Engineering[END_REF]. Another important contribution in the late 19th century was by another German researcher named Bauschinger (see [START_REF] Bauschinger | Über die Veränderungen der Elastizitätsgrenze und der Festigkeit des Eisens und Stahls durch Strecken, Quetschen, Erwärmen Abkühlen und durch oftmals wiederholte Belastung[END_REF], who discovered the difference in the elastic limits of metals between reverse and monotonic loadings. This was one of the earliest investigations on cyclic hardening and softening phenomenon.

In the early 20th century, first important research contribution was from a British researcher named Erwing, who along with his colleagues investigated the crack nucleation and propagation mechanism in polycrystalline materials (see Ewing and[START_REF] Ewing | The fracture of metals under rapid alterations of stress[END_REF][START_REF] Ewing | Experiments in micro-metallurgy-effects of strain. preliminary notice[END_REF]. The American engineer Basquin also became famous in the early 20th century, after his proposal of empirical laws to characterise S-N curves in 1910. He discovered that in a logarithmic scale the stress has a linear relationship with number of fatigue cycles over a large stress range (see [START_REF] Basquin | The exponential law of endurance tests[END_REF]. The early part of the 20th century was dominated by the British and American engineers, with notable contributions from [START_REF] Smith | A fatigue testing machine[END_REF], [START_REF] Haigh | A new machine for alternating load tests[END_REF], [START_REF] Gough | The Fatigue of Metals[END_REF], [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF], [START_REF] Inglis | Stresses in a plate due to the presence of cracks and sharp corners[END_REF], [START_REF] Kommers | Repeated stress tests of steel[END_REF], [START_REF] Moore | The fatigue of metals[END_REF], and many more. The most important of them was Griffith, who with his investigations on glass, laid the foundation of fracture mechanics. By the early quarter of the century, fatigue became one of the most important fields of study. Another noteworthy contribution of the period was from a Swedish scientist named Palmgren who investigated and proposed a theory on accumulation of damage in material during fatigue, especially for variable amplitudes. His contribution (see [START_REF] Palmgren | Die Lebensdauer von Kugellagern[END_REF] was extended and formulated by the American scientist Miner (see [START_REF] Miner | Cumulative damage in fatigue[END_REF], who developed what is now called the famous "Palmgren-Miner linear damage rule". Another landmark in this period was by the Swedish engineer Weibull, who extended his theory of static strength to fatigue and conducted several tests to obtain numerical data of number cycles to failure, to quantify his proposed distribution, which later on became famous as "Weibull distribution" (see [START_REF] Weibull | A Statistical Theory of the Strength of Materials[END_REF]. One of the most important German researchers of this era was Thum, who pioneered the idea of "Gestaltfestigkeit", which according to Thum was "a strength value which depends on the magnitude and type of loading as well as on the material, and especially on the component's shape" (see [START_REF] Schütz | A history of fatigue[END_REF]. He had a major contribution in the area of fatigue and his studies in the areas of corrosion fatigue, hardenings, residual stresses, fretting and such others were revolutionary (see [START_REF] Thum | Festigkeitsprüfung bei schwingender Beanspruchung[END_REF]. The concept of fatigue stress concentration factor was developed by him. Another important personality during this time was the German engineer Gassner who was the pioneer in variable amplitude fatigue tests (see [START_REF] Gassner | Auswirkung betriebsähnlicher Belastungsfolgen auf die Festigkeit von Flugzeugbauteilen[END_REF]. The idea of "Betriebsfestigkeit", which translates to operational fatigue strength was first formulated by him. Another important scientist of this time was a German engineer named Neuber, who published the first comprehensive book (see [START_REF] Neuber | Kerbspannungslehre, Grundlage für eine genaue Spannungsrechnung[END_REF] about the theoretical calculations of stress concentration factors. He was one of the first to study notch effects on fatigue. The next landmark came in the later half of the century, when two American engineers Coffin and Manson, independently investigated the dependency of plasticity on cyclic damage (see [START_REF] Coffin | A study of the effects of cyclic thermal stresses on a ductile metal[END_REF][START_REF] Manson | Behaviour of materials under conditions of thermal stress[END_REF]. Their works were the first efforts in the investigations of low cycle fatigue and led to the strain-based measurements, compared to the stress-based measurements that had been used traditionally. The ideas of Griffith used for brittle materials were extended by Irwin (see [START_REF] Irwin | Analysis of stresses ans strains near the end of a crack traversing a plate[END_REF] for ductile materials, giving the concept of stress singularity using stress intensity factor. Paris, with his co-workers in the 1960s found that the fatigue crack growth rate can be related to the stress intensity factor range (see Paris andErdogan, 1963, Paris et al., 1961). Japanese researchers Matsuishi and Endo developed a method in 1968 to represent a non-cyclic loading in equivalent number of cycles (see [START_REF] Matsuishi | Fatigue of metals subjected to varying stress[END_REF]. This method became famous as rainflow counting. During the similar time frame, a German scientist named Elber, experimentally found out that under cyclic tensile loading, crack closes before the load reaches zero (see [START_REF] Elber | Fatigue crack propagation: some effects of crack closure on the mechanisms of fatigue crack propagation under cyclic tensile loading[END_REF]. His proposed "crack closure" ideas has paved the way for crack-growth hypotheses since his discovery. A notable person during this period was the Dutch scientist Schijve, who greatly contributed to the fatigue studies of aircraft, even for variable loading (see Schijve and[START_REF] Schijve | Crack propagation tests based on a gust spectrum with variable amplitude loading[END_REF][START_REF] Schijve | Fatigue crack growth in lugs and the stress intensity factor[END_REF]. A contemporary British scientist Pearson was the first to identify the "short crack problem", which later on became an important are in fatigue crack growth (see [START_REF] Pearson | Initiation of fatigue cracks in commercial aluminium allys and the subsequent propagation of very short cracks[END_REF]. During the late 70s and early 80s, an Indian scientist named Suresh along with his co-workers investigated various types of crack growth and crack closure mechanisms (see [START_REF] Ritchie | Near-threshold fatigue crack growth in 2 1/4 Cr-1 Mo pressure vessel steel in air and hydrogen[END_REF][START_REF] Suresh | A geometrical model for fatigue crack closure induced by fracture surface morphology[END_REF]. The numerical simulation of fatigue got a boost when researchers began using continuum damage mechanics to model and simulate fatigue phenomena. The onset of continuum damage mechanics although happened in 1958 through a Russian scientist named Kachanov (see [START_REF] Kachanov | In time to rupture in creep conditions (in Russian)[END_REF], the usage of CDM for fatigue modelling and simulations took place a lot later. The initial usage of CDM for fatigue computation started during late 80s and early 90s. Two French scientists, Lemaitre (see [START_REF] Lemaitre | Coupled elasto-plasticity and damage constitutive equations[END_REF][START_REF] Lemaitre | Damage 90: a post processor for crack initiation[END_REF] and Chaboche (see [START_REF] Chaboche | A non-linear continuous fatigue damage model[END_REF] were pioneers in the usage of damage mechanics in the prediction of fatigue life. Later on a Japanese researcher named Murakami, who also played a part in the development of CDM in the 70s (see [START_REF] Murakami | A constitutive equation of creep damage in polycrystalline metals[END_REF], used damage mechanics for fatigue computation (see [START_REF] Tanaka | Study on Evolution of Internal Damage in CFRP in Fatigue Process[END_REF]. Thereafter a new idea of cycle jumping methods has been proposed to compute fatigue. Initially in Lesne and Savalle (1989) and later on in [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF], jump-in-cycles method has been used extensively in the 90s and 2000s. Fish and his colleagues since the early 2000s have used an alternative method to simulate fatigue with CDM (see Fish and[START_REF] Fish | A nonlocal multiscale fatigue model[END_REF]Fish, 2004a,b) using the temporal homogenisation technique developed by a French scientist called Aubry in the late 80s (see [START_REF] Guennouni | Réponse homogénéisée en temps de structures sous chargements cycliques[END_REF].

Although, the complete history of fatigue was not covered, only certain landmarks are highlighted, which according to the author are most significant. For more comprehensive historical perspective on fatigue, [START_REF] Schütz | A history of fatigue[END_REF], [START_REF] Suresh | Fatigue of materials[END_REF] are suggested.

Different forms of fatigue in the engineering world

The driving factor of fatigue is cyclic or fluctuating stresses or strains, which can arise due to many factors in real-life engineering structures. One of the most classical form is the mechanical fatigue, in which the internal repetitive stresses occurs under the fluctuations of externally applied mechanical loading (forces or displacements). Most of the civil structures like bridges, dams, and such others, or mechanical components like machine tools, pumps, vibration dampers etc. are susceptible to mechanical fatigue. Thermal fatigue occurs under the influence of fluctuating temperature. This induces internal fluctuating stresses and strains which result into fatigue failure. This can happen in high temperature super heaters, nuclear reactors, and such others. A fluctuation of both mechanical and thermal quantities, for instances in internal combustion engines, steam or gas turbines, results into thermo-mechanical fatigue. Electrical fatigue occurs generally for electro-active materials such as piezoelectrics under the influence of fluctuating electric field. This kind of materials produces fluctuating internal mechanical stresses or strains due to electromechanical coupling resulting into fatigue of the material. Thermo-electrical fatigue and electro-mechanical fatigue occurs for electro-active materials when electrical fatigue is combined with thermal and mechanical fatigue respectively. Creep fatigue takes place when a structure is subjected to mechanical fatigue at high temperature. Corrosion fatigue occurs at chemically aggressive environment. There can still be further different varieties depending on the different forms of loading.

Phases of fatigue life

Irrespective of the type of loading, fatigue phenomenon is considered to be originating from local yielding of the material (see [START_REF] Cui | A state-of-the-art review on fatigue life prediction methods for metal structures[END_REF][START_REF] Schijve | Fatigue of Structures and Materials[END_REF][START_REF] Sobczyk | Random Fatigue: From Data to Theory[END_REF]. Due to local stress concentration at grain level of a material, dislocations or invisible micro-cracks occur in slip bands. These dislocations, under the influence of cyclic stress, migrate and coalesce to form a micro-crack. The microcracks grow to form macro-crack. This macro-crack then propagates through the structure, generally perpendicular to the maximum principal stress direction, which ultimately leads to the failure of the structure (see [START_REF] Cui | A state-of-the-art review on fatigue life prediction methods for metal structures[END_REF][START_REF] Ottosen | Continuum approach to high-cycle fatigue modeling[END_REF][START_REF] Schijve | Fatigue of Structures and Materials[END_REF][START_REF] Sobczyk | Random Fatigue: From Data to Theory[END_REF]. Thereby the initiation and propagation of macro-crack are the two paramount stages of fatigue process.

Various steps of fatigue life within the two principal stage as mentioned in [START_REF] Oller | A continuum mechanics model for mechanical fatigue analysis[END_REF], Schijve (2008), [START_REF] Sobczyk | Random Fatigue: From Data to Theory[END_REF] are concisely depicted in fig. 1.1. The two principal stages and the physical phenomena behind them are herein detailed in brief (see Schijve, 2008, Sobczyk and[START_REF] Sobczyk | Random Fatigue: From Data to Theory[END_REF], for more details).

Crack initiation phase of fatigue is a surface phenomenon. Plasticity, being a local phenomenon, for initial cycles, is restricted to a small number of grains and can occur more effectively if there are no surrounding grains to hinder the plastic deformation, i.e. at the surface. This micro-plasticity is a consequence of cyclic slip originating from cyclic shear stress, which is not homogeneous but varies from grain to grain. Cyclic loading creates a non-reversible slip band on parallel slip planes. The first initiation of a micro-crack is experienced to develop along a slip band in the form of intrusion or extrusion (see [START_REF] Schijve | Fatigue of Structures and Materials[END_REF][START_REF] Suresh | Fatigue of materials[END_REF], for more details). These micro-defects will be initiated depending on

Crack initiation period

Nucleation

Micro-crack growth

Crack propagation period

Macro-crack growth Failure

Total fatigue life

Figure 1.1: Different stages of fatigue process surface imperfections, inherent material defects, geometric discontinuities, high stress concentrations, and such other conditions that facilitate cyclic slip. These micro-cracks grow into adjacent grains with the growth directions deviating from initial slip band orientations with the grain boundaries providing resistance against growth. When the number of grains along the micro-crack front becomes sufficiently large, the growth happens continuously with the growth rate depending on the crack growth resistance of the material. The completion of the crack initiation phase is assured when the growth of the micro-crack becomes independent of surface conditions. Crack propagation phase starts when the growth of crack is controlled by the crack growth resistance which is a bulk property and not a surface phenomenon. The crack generally grows perpendicular to the direction of the maximum principal stress through the structure until it fails completely. For brittle materials the propagation can be rapid, however for ductile material it is slower and accompanied by local plastic zones around the crack tip which provides crack growth resistance.

The percentage of total life that is contributed by these two phases depends on the type of material along with loading itself. This total life can be used as quantity that classifies different forms of fatigue.

Different domains of fatigue

The classical way to split the domains of fatigue is with respect to the number of cycles to failure (see [START_REF] Kazymyrovych | Very high cycle fatigue of engineering materials : A literature review[END_REF]. Structures that are subjected to low amplitude loading such that the induced stress range is below the yield stress, are macroscopically elastic. However, plasticity and damage can occur at the micro-scale. In such cases, the number of load cycles to failure is generally more than 10 6 (see [START_REF] Kazymyrovych | Very high cycle fatigue of engineering materials : A literature review[END_REF], and this is referred to as high cycle fatigue (HCF). For high amplitude loading, however, the induced stresses exceed the elastic limit and the dominating mechanism is plasticity at the macro-scale. In these cases, the number of cycles is relatively low and generally below 10 6 and they are called low cycle fatigue (LCF). These demarcations are also relative to the authors, as for instance, [START_REF] Pook | Metal Fatigue: What It Is, Why It Matters[END_REF] considers LCF if number of cycles to failure N f < 10 4 and above that it is HCF, [START_REF] Davoli | Principles of current methodologies in high-cycle fatigue design of metallic structures[END_REF] however introduces 10 5 as the demarcation point between LCF and HCF. [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF] on the other hand considers LCF if N f < 10 4 and HCF if N f > 10 5 and in between either of them can be considered depending on the particular case.

Apart from these distinctions, if the number of cycles to failure is extremely high, it is referred to as very high cycle fatigue (VHCF). The demarcation point in this case is also debatable, [START_REF] Kazymyrovych | Very high cycle fatigue of engineering materials : A literature review[END_REF] considers it to be more than 10 8 , but [START_REF] Pyttel | Very high cycle fatigue -Is there a fatigue limit[END_REF] considers VHCF if N f > 10 7 . In this context there also exist what is called very low cycle fatigue (VLCF), also referred to as extremely or ultra low cycle fatigue. This is applicable when N f is extremely low. [START_REF] Skrzypek | Modeling of material damage and failure of structures: theory and applications[END_REF] defines N f to be less than 10, however according to [START_REF] Xue | A unified expression for low cycle fatigue and extremely low cycle fatigue and its implication for monotonic loading[END_REF] this is less than 100. Irrespective of these variable demarcations, it is universally agreed that LCF and VLCF are governed by macro-plasticity and HCF and VHCF are governed by macro-elasticity (especially for metals). It has to be noted that VLCF and VHCF are subsets of LCF and HCF respectively, and the differentiation is based more on the number of cycles instead of any difference in the governing physical phenomena.

Although fatigue can be classified with respect to number of cycles, however for many practical cases "the cycles" prescribed by the load history are not well defined. Hence it is necessary to introduce different forms of load fluctuations.

Different types of load fluctuations

It is necessary at this point to introduce the different types of fluctuations or the load forms that are generally considered for fatigue studies. The most primitive form is constant amplitude (CA) loads where the loads are completely repetitive after certain time instances called the time period. All the cycles in such type of loading are identical (see [START_REF] Pook | Metal Fatigue: What It Is, Why It Matters[END_REF]. The loading in such cases is defined in general by certain terminologies which are classical in CA fatigue. Considering for instance a stress-controlled case where the loading is between maximum stress σmax and minimum stress σ min , the stress amplitude is then σa = σmax-σ min 2 and the mean stress is σm = σmax+σ min 2 . The stress range S or ∆σ is defined as ∆σ = 2σa = σmax -σ min . The loading may also be expressed in terms of the stress ratio R = σ min σmax , and the amplitude ratio A = σa σm . Certain types of CA loadings are shown in fig. 1.2. Although, these terminologies are classically used for stress-based loading, similar expressions exist in terms of strains as well (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF].

Variable amplitude (VA) loading however does not have a fixed definition, but encompasses everything which does not fall under CA loading (see [START_REF] Schijve | Fatigue under variable-amplitude loading[END_REF]. [START_REF] Schijve | Fatigue under variable-amplitude loading[END_REF] states that the most simple form of VA loading is when several blocks of CA loadings are applied. Similarly, loadings that have constant amplitude but variable mean stress or loadings which have constant mean values but variable amplitudes also fall under VA loading. [START_REF] Richard | Fatigue Crack Growth: Detect -Assess -Avoid[END_REF]

Existing fatigue approaches

Various authors over the years have used different ways of classifying different methods that have been developed for fatigue. In [START_REF] Suresh | Fatigue of materials[END_REF], the methods have been grouped under "total-life approaches" and "defect-tolerance approaches". However, in [START_REF] Cui | A state-of-the-art review on fatigue life prediction methods for metal structures[END_REF], the methods have been grouped under "cumulative fatigue damage theories" and "fatigue crack propagation theories". All the approaches where crack behaviour is not taken into consideration, or to be more specific, all the approaches which do not consider a pre-existing crack can be categorised under "cumulative fatigue damage theories", because essentially these methods quantify the deterioration of material properties without involving fracture mechanics. "Fatigue crack propagation theories" or "defect-tolerance approaches" both essentially indicate those methods which assume the presence of macroscopic flaws and compute the number of cycles needed to propagate them. The "total-life approaches" described in [START_REF] Suresh | Fatigue of materials[END_REF] encompass all the phenomenological continuum-based approaches where the total life is expressed as a function of measurable quantities like stress or strain range. The total life in this case is the time to initiate a dominant crack and the time to propagate the rupture to failure.

It is perhaps important at this point to introduce what is known as "safe-life" and "fail-safe" design approaches for fatigue. These concepts came into existence in the 1950s for the design of aircraft components (see [START_REF] Schütz | A history of fatigue[END_REF]. In the case of safe-life approach, the component in question is subjected to a load spectrum that mimics the actual service conditions and a fatigue life is obtained which is invariably modified using a factor of safety. Several tests can be considered using various influential factors to get an idea of the useful fatigue life. The component is thereby rejected from service when it reaches its "safe life", even though there is no failure. The "fail-safe" concept however stems from the idea that even if an individual component fails, the other components should have enough structural integrity so that the whole structure can operate safely. Periodic inspections of individual components are mandatory in this design philosophy, such that damaged members can be repaired or replaced (see [START_REF] Suresh | Fatigue of materials[END_REF].

The following classification uses the designation attributed by [START_REF] Cui | A state-of-the-art review on fatigue life prediction methods for metal structures[END_REF].

Cumulative fatigue damage theories

The basic idea behind these theories is the fact that the inherent microscopic parameters governing the damage of the material is related to macroscopic quantities like stress or strain. Different damage theories have been developed over the years based on different macroscopic quantities.

Stress-based approach or stress life approach

This was first introduced by Wöhler in 1860. The basic idea of this method is to generate what is called the S-N curves. The specimen is loaded at specified cyclic stresses till the specimen ruptures completely, and the number of cycles to failure (N f ) is observed. The stress range (∆σ) or the stress amplitude (σa) is plotted with respect to the number of cycles to failure to obtain what is called the S-N curve (fig. 1.5). The concept of endurance limit or fatigue limit (σe) is an outcome of the asymptotic behaviour of the S-N curves, which signifies that if the specimen is loaded with a stress amplitude below the endurance limit (σe) then it will have infinite life.

σ a log (N f ) σ e Figure 1.5: A typical S-N curve
The numerical formulation of the S-N curves is given by the Basquin's equation as

σa = σ ′ f 2N f b , (1.1)
where σ ′ f is the fatigue strength coefficient and b is the fatigue strength exponent. However, the Basquin's equation can not predict the complete S-N curve especially for the low cycle regime. A relatively new reformulation of the Basquin's was made in [START_REF] Kohout | New functions for a description of fatigue curves and their advantages[END_REF] to describe the complete S-N curve. The S-N curve is influenced by a lot of factors such as frequency of loading, ambient temperature, notches in the specimens, and such others. The most important influencing factor is the "mean stress effect", where with the increase in the mean stress, the number of cycles to failure decreases for a given stress amplitude. The classical way of representing the "mean stress effect" is by using constant life diagrams (fig. 1.6). The most important mean stress models used are developed in Gerber (1874), [START_REF] Goodman | Mechanics Applied to Engineering[END_REF][START_REF] Soderberg | Factor of safety and working stress[END_REF]. Although these diagrams do not provide any information on the number of cycles to failure, they can predict any potential to fatigue failure. These diagrams basically depict the maximum mean stress that can be achieved for a given stress amplitude without fatigue failure. The Basquin's equation has been modified in [START_REF] Morrow | Fatigue properties of metals[END_REF] to include mean stress effect as,

σa = σ ′ f -σm 2N f b . (1.2)
The next important consideration in stress-based approaches is to deal with multiaxial cyclic stresses. The most easy way to analyse multiaxial fatigue is to consider effective stresses. For instance, if σ i,a are the amplitudes of the principal stresses σ i , with i = 1, 2, 3, the effective stress amplitude can be computed using the von Mises criterion as

σe,a = 1 √ 2 (σ 1,a -σ 2,a ) 2 + (σ 2,a -σ 3,a ) 2 + (σ 3,a -σ 1,a ) 2 . (1.3)
This effective stress can be used in the Basquin's equation to obtain the number of cycles to failure. In the same way the effective mean stress can be calculated which can then either be used in constant life approach or in Morrow's equation. One of the major drawbacks of this type of effective stress approach is the fact that the effect of tension and compression might not be adequately captured. Thereby the critical plane approach is generally used. Critical plane for multiaxial fatigue can be loosely defined as the plane on which damage and subsequently cracks will occur. Although, the orientation of the critical plane for general multiaxial and non-proportional stress state is not straightforward, several simplified theories have been proposed over the years for multiaxial loading. The first one was from [START_REF] Crossland | Effect of large hydrostatic pressure on the torsional fatigue strength of an alloy steel[END_REF] who proposed a safety criterion based on effective stress,

J 2,a + ᾱσ H,max ≤ β. (1.4)
J 2,a is the amplitude of the second invariant of the deviatoric stress and σ H,max = σ H,m + σ H,a , with σ H,m and σ H,a being the mean and amplitude of the hydrostatic stress. ᾱ and β are material parameters. This criterion was modified by [START_REF] Sines | Behavior of metals under complex static and alternating stresses[END_REF] to use the mean hydrostatic stress instead of the maximum hydrostatic stress, J 2,a + ᾱσ H,m ≤ β.

(1.5)

Several other theories exist especially in Dang [START_REF] Van | Sur la résistance à la fatigue des métaux[END_REF] with the maximum shear stress and instantaneous hydrostatic stress, and in [START_REF] Papadopoulos | A comparative study of multiaxial high-cycle fatigue criteria for metals[END_REF] where volume averaged stress quantities are used.

Strain-based approach or strain life approach

The idea of stress-based analysis is mainly restricted to elastic deformation. However, for practical cases, structures generally undergo localised plastic deformation, therefore it is more appropriate under these circumstances, to consider a strain-based analysis. The first formulated idea of strain-based analysis is given by the Coffin-Manson relationship,

ε p a = ε ′ f 2N f c , (1.6)
where, ε ′ f is the fatigue ductility coefficient, ε p a is the plastic strain amplitude, and c is the fatigue ductility exponent. Considering additive strain decomposition relation

εa = ε e a + ε p a , (1.7)
with εa and ε e a being the total and elastic strains respectively, and the linear elastic law

ε e a = σa E , (1.8)
with E being the modulus of elasticity, the Basquin's equation and the Coffin-Manson rule can be combined to give

εa = σ ′ f E 2N f b + ε ′ f 2N f c . (1.9)
This equation can be plotted as ε-N curves as shown in fig. 1.7. This figure defines a transition life N f,t which is the fatigue life when elastic and plastic strain amplitudes are equal.

ε a ε ′ f σ ′ f E ε p a ε e a b 1 c 1 log (2N f ) log (2N f,t ) Figure 1.7: A typical ε-N curve
Many modifications have been proposed for the Coffin-Manson relation to replicate accurately experimental ε-N curves, initially in [START_REF] Manson | Fatigue: a complex subject-some simple approximations[END_REF] and later on in [START_REF] Roessle | Strain-controlled fatigue properties of steels and some simple approximations[END_REF].

In cases of multiaxial loading, similar to the stress-based analysis, the simplest way is to calculate an effective strain amplitude as

εe,a = 1 √ 2 (1 + ν) (ε 1,a -ε 2,a ) 2 + (ε 2,a -ε 3,a ) 2 + (ε 3,a -ε 1,a ) 2 , (1.10)
with ν being the Poisson's ratio, ε i,a the principal strain, and i = 1, 2, 3. This can then be used in eq. (1.9). Similar criterion can be derived using Tresca theory. However these effective strain measures do not capture the mean stress effect, hence an energy-based criterion was proposed in [START_REF] Smith | A stress-strain function for the fatigue of metals[END_REF],

σmaxεa = σ ′ f 2 E 2N f 2b + ε ′ f σ ′ f 2N f b+c .
(1.11)

Cumulative damage rules

The first damage rule is the linear summation of fatigue damage that is being accumulated with load level. If n i is the number of cycles for which a load level of σ a,i is applied in a sequence of m load levels, the Palmgren-Miner linear damage rule (LDR) states that,

D = m i=1 r i = m i=1 n i N f,i , (1.12)
where D is the damage variable, N f,i is the number of cycles to failure for load level σ a,i which can be obtained from S-N curves, and r i is the cycle ratio corresponding to load level σ a,i . The assumption of LDR is that at each load level the damage is accumulated in the form of cycle ratio and failure occurs if

m i=1 r i = 1. (1.13)
This theory does not take into account the order of the load levels σ a,i , and does not comply with experimental findings.

The first important modification of the LDR was in [START_REF] Marco | A concept of fatigue damage[END_REF], where non-linear damage rule (NLDR) is represented by a power relationship,

D = m i=1 r x i i = m i=1 n i N f,i x i , (1.14)
where x i is an exponent depending on load level σ a,i . This power law mimics the experimental behaviour of m i=1 r i < 1 for high-to-low (H-L) loading and m i=1 r i > 1 for low-to-high (L-H) loading. The comparison of LDR and NLDR for H-L and L-H loadings is shown in fig. 1.8. Another important cumulative damage theory is the two stage linear damage theory developed by the combined work of Langer and Grover (see [START_REF] Grover | An observation concerning the cycle ratio in cumulative damage[END_REF][START_REF] Langer | Fatigue failure from stress cycles of varying amplitude[END_REF]. The idea is basically to separate the damage process into crack initiation, N I = αN f ; and crack propagation, N P = (1 -α) N f ; where α is a life fraction factor. The LDR can then be applied to either stage. [START_REF] Manson | Fatigue: a complex subject-some simple approximations[END_REF] provided the double linear damage rule (DLDR) through which the equations can be separated into

0 0 1 D 1/3 2/3 1 A B C D E F G H ni Nf,i σa,1 LDR σa,2 σa,3 σa,1 > σa,2 > σa,3 σa,1 -→ σa,2 -→ σa,3 3 i=1 ni Nf,i = AC + DF + GH < 1 σa,3 -→ σa,2 -→ σa,1 3 i=1 ni Nf,i = AC + BF + EH > 1
N I = N f -P N 0.6 f ,
(1.15a)

N P = P N 0.6 f , (1.15b)
where P is a fatigue life coefficient. The LDR can be applied thereby one for the initiation of crack The connection between the hysteretic area and fatigue behaviour, first discovered by Inglis (see [START_REF] Inglis | Hysteresis and fatigue of Wöhler rotating cantilever specimen[END_REF], has led to a lot of research in describing cumulative damage variables in terms of energy-based parameters. Energy-based damage models can predict the damage accumulation in notched specimen, based on Glinka's rule; which basically states that, if the plastic zone around a notch is surrounded by elastic zone, the elasto-plastic strain energy is approximately equal to the strain energy if the material is assumed to be completely elastic. Energy-based models are able to take mean stress effects and multiaxial behaviour into account. One of the earliest energy-based model was in [START_REF] Kujawski | A cumulative damage theory of fatigue crack initiation and propagation[END_REF], where plastic strain energy density as a parameter was used. Later on it was found that mean stress effect can not be directly incorporated using plastic strain energy. To rectify this problem, a total strain energy-based approach was introduced in [START_REF] Golos | A total strain energy density theory for cumulative fatigue damage[END_REF]. Irrespective of the type of model used, the main idea is to formulate a power law similar to the Basquin's equation to obtain energy-life relation (see [START_REF] Fatemi | Cumulative fatigue damage and life prediction theories: a survey of the state of the art for homogeneous materials[END_REF]. Another important development was in [START_REF] Leis | A nonlinear history-dependent damage model for low cycle fatigue[END_REF], where the damage parameter similar to Smith (eq. (1.11)) was linked as

D = 4σ ′ f E 2N f 2b 1 + 4ε ′ f σ ′ f 2N f b 1 +c 1 , (1.18)
with b 1 and c 1 being obtained from instantaneous strain hardening exponent. It is found in [START_REF] Niu | Hardening law and fatigue damage of a cyclic hardening metal[END_REF] that the strain hardening coefficient tends to change while the strain hardening exponent remains the same during cyclic process. Thereby, a damage function was proposed as

D = r 1/(n ′ +α) , (1.19)
where α = (∆σ∆εp/4) 2b √ a with a and b being constants, and n ′ is the cyclic strain hardening exponent. Another important method was postulated in [START_REF] Kliman | Fatigue life prediction for a material under programmable loading using the cyclic stress-strain properties[END_REF], who proposed a damage model for repeated blocks of harmonic load cycles. Considering W b is the hysteretic energy for a block and W f,b is the energy at fracture for the block, the damage fraction for the block is given by

D = W b W f,b .
(1.20)

1.6.1.5 Continuum damage mechanics (CDM) theories

The mechanics of damage is the study, through mechanical variables, of the mechanisms involved in the deterioration when the materials are subjected to loading. The idea of degradation of material being quantified by internal variable in a thermodynamically consistent continuum framework was first developed by [START_REF] Kachanov | In time to rupture in creep conditions (in Russian)[END_REF][START_REF] Kachanov | Introduction to continuum damage mechanics[END_REF]. This is a relatively new branch, which is an extension of the classical continuum mechanics. The first usage of CDM for fatigue computation was in [START_REF] Chaboche | A non-linear continuous fatigue damage model[END_REF]. The proposed damage law is of the form

D = 1 -1 -r 1/(1-α) 1/(1+β) , (1.21)
where, β is a material parameter, α is a function of the stress state, and r is the cycle ratio. Many notable works have followed since, and well-established continuum damage mechanics theories can be found in [START_REF] Lemaitre | A course on damage mechanics[END_REF], [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF], [START_REF] Murakami | Continuum Damage Mechanics: A Continuum Mechanics Approach to the Analysis of Damage and Fracture[END_REF]. CDM will be dealt in details in Chap. 3.

Fatigue crack propagation theories

Although many authors (see [START_REF] Fatemi | Cumulative fatigue damage and life prediction theories: a survey of the state of the art for homogeneous materials[END_REF] do not distinguish between fatigue damage theories and crack propagation theories, a distinction between damage mechanics and fracture mechanics concepts following the demarcation from [START_REF] Cui | A state-of-the-art review on fatigue life prediction methods for metal structures[END_REF], [START_REF] Suresh | Fatigue of materials[END_REF] has been made here. However the counter argument in [START_REF] Fatemi | Cumulative fatigue damage and life prediction theories: a survey of the state of the art for homogeneous materials[END_REF] is also valid as crack is initiated at the very beginning of the loading and the entire life of the structure is spent in propagating the crack. The initiation and propagation phase can be distinguished from the length of the crack (which is again not fixed). However based on classical fracture mechanics theories, a distinction is made on the fact that crack propagation theories are based on the presence of a pre-crack in the structure.

Based on the theories explained in [START_REF] Cui | A state-of-the-art review on fatigue life prediction methods for metal structures[END_REF], [START_REF] Fatemi | Cumulative fatigue damage and life prediction theories: a survey of the state of the art for homogeneous materials[END_REF], [START_REF] Suresh | Fatigue of materials[END_REF], it can be divided into three parts based on the crack length.

Long crack growth

The first law about fatigue crack growth is based on the findings of Paris and Erdogan (see [START_REF] Paris | A critical analysis of crack propagation laws[END_REF] [START_REF] Donahue | Crack opening displacement and the rate of fatigue crack growth[END_REF],

da dN = C (∆K -∆K th ) m , (1.24)
where ∆K th is the threshold factor below which there will be no crack growth, and C and m are material parameters. For the region of unstable crack growth, the growth rate is extremely high leading to either brittle fracture or ductile tearing. A relation expressing unstable crack growth rate is proposed in [START_REF] Forman | Numerical analysis of crack propagation in cyclic-loaded structures[END_REF] as,

da dN = C (∆K) m (1 -R) Kc -∆K , (1.25)
where Kc is the fracture toughness and R is the load ratio. These three regions are schematically represented in fig. 1.10. Other works include the contribution of Elber (see [START_REF] Elber | Fatigue crack closure under cyclic tension[END_REF] who introduced crack closure concepts which was later on improved in [START_REF] Donald | An evaluation of ∆K ef f estimation procedure on 6061-T6 and 2024-T3 aluminum alloys[END_REF], [START_REF] Hertzberg | Crack closure: Correlation and confusion[END_REF], [START_REF] Vasudevan | A review of crack closure, fatigue crack threshold and related phenomena[END_REF] and especially in [START_REF] Kujawski | Enhanced model of partial crack closure for correlation of R-ratio effects in aluminum alloys[END_REF], where the concept of partial crack closure was introduced. The main ideas and developments about long crack growths are based on linear elastic fracture mechanics (LEFM).

Physically small crack growth

The growth of physically small cracks is governed by elastic plastic fracture mechanics (EPFM). The first idea of EPFM was proposed in [START_REF] Tomkins | Fatigue crack propagation-an analysis[END_REF], who equated the crack growth rate as a function of the plastic strain, for high strain fatigue, i.e.

da dN = B ∆εp √ πa m -Ξ, (1.26)
where B and m are material parameters and Ξ is the threshold condition. Many theories have been proposed especially by modifying the long crack growth laws to include plastic deformation. They were able to predict different EPFM processes quite effectively. One such is the model proposed in [START_REF] Mcevily | A modified constitutive relation for fatigue crack growth[END_REF] which matched quite well with experimental findings.

Microstructurally small crack growth

This method was developed to measure crack growth at the level of microstructure. The first idea of micro-structural fracture mechanics (MFM) was first introduced in [START_REF] Hobson | Two phases of short crack growth in a medium carbon steel[END_REF] and later on in [START_REF] Navarro | A microstructurally short fatigue crack growth equation[END_REF]. The concept of MFM is also based on EPFM and the growth of the crack can be expressed as

da dN = B (∆γ) α (d -a) , (1.27)
where B and α are material parameters, ∆γ is range of shear strain, and d is a microstructural dimension. As both microstrually small crack (MSC) and physically small crack (PSC) are governed by EPFM, [START_REF] Miller | Metal fatigue-a new perspective[END_REF]Hobson (see Hobson, 1986, Miller, 1992) proposed combined laws which are of the form

da dN = B (∆γ) α (d -a) , for MSC, a 0 ≤ a ≤ at, (1.28a) da dN = A (∆γ) β a -C, for PSC, at ≤ a ≤ a f , (1.28b)
where a 0 is the initial defect size, at is the transition crack length, a f is the final crack length, A and α are material parameters, C is the crack growth rate at threshold condition.

Although several other methods exist, only the most popular techniques of traditional or contemporary approaches have been covered. For a more extensive review, [START_REF] Cui | A state-of-the-art review on fatigue life prediction methods for metal structures[END_REF], [START_REF] Fatemi | Cumulative fatigue damage and life prediction theories: a survey of the state of the art for homogeneous materials[END_REF] are suggested.

Cyclic elasto-(visco)plasticity

Structures when subjected to cyclic elasto-(visco)plasticity, three typical phenomena may occur (see [START_REF] Maugin | The Thermomechanics of Plasticity and Fracture[END_REF][START_REF] Weiß | Simulation of ratcheting and low cycle fatigue[END_REF]. One of the behaviour corresponds to the case where the stress-strain response curve does not close and the strain increases with each cycle to induce a ratcheting effect (see fig. 1.11a). Otherwise, the response curve may end up in a stable hysteresis loop leading to plastic shakedown or accumulation (see fig. 1.11b). It is also possible in some instances where this stable loop can completely flatten such that the structure behaves elastically, generating elastic shakedown or adaptation (see fig. 1.11c). In many cases the idea is to find the stabilised cycle, and thereafter the Coffin-Manson rule can be utilised to find the number of cycles to failure. The search for the stable loop however can be expensive, and hence some sophisticated numerical methods are used to accelerate the simulation. In general these acceleration techniques are divided into three groups (see [START_REF] Nasri | Dimensional reduction for the simulation of metal fatigue[END_REF] -Methods where the evolutions of the quantities of interest are estimated at each point of the structure, e.g. direct cyclic method, cycle skip method, Zarka method, LATIN method;

-Methods where the evolutions of the quantities of interest are estimated only at the critical points of the structure, e.g. Neuber type methods;

-Methods based on model reduction approaches.

The direct cyclic method as introduced in Pommier ( 2003) is based on calculation of a single cycle with imposed periodic condition. The global admissibilities and local behaviour laws are iteratively calculated until convergence is achieved. Cycle skip methods, for plasticity, essentially deal with the fact that full blocks of cycles can be skipped by successive extrapolations [START_REF] Saï | Modèles à grand nombre de variables internes et méthodes numériques associées[END_REF]. The Zarka method basically considers a structural transform parameter to find the stabilised cycle from the history of its elastic state, both for elastic and plastic shakedown (see Zarka and Caiser, 1979, for details). The variant of the LATIN method as introduced in Cognard and [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF] used temporal model reduction technique to simulated large number of cycles for elasto-viscoplastic materials. The Neuber method is similar to the energy-based approach where most of the structure is considered elastic and Neuber correction is applied only on the local plastic zones (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. Nasri (2017) used several PGD-based model reduction methods to obtain the stabilised cycle.

Concluding remarks

From all the methods that have been discussed, the S-N curve, or ε-N curve approaches might seem straightforward, but they are extremely time consuming. Experimental S-N or ε-N curves are generally built using uniaxial tests, for simple specimens. In practical cases involving mixed loading, with localised plasticity with hardening or softening behaviour, therefore resulting in the stress or strain amplitude being depended on space and time, making the S-N or ε-N curves irrelevant. The cumulative damage theories, like the LDR, NLDR, or DLDR although seem to be efficient, rely on the information of S-N or ε-N curves, which for practical structures are not available. Although, NLDR and DLDR give better approximations for block loadings than LDR, the damage variable for all these cases is a function of the cycle ratio and does not represent the physical mechanism of fatigue damage. Energybased approaches although are thermodynamically consistent models, the damage variable is calculated as a post process of the hystertic area, and not in the form of internal variable. CDM however considers damage to be an intrinsic variable in the fundamental lemmas of thermodynamics, and thereby provides a thermodynamically consistent quantification of the degradation of the material due to fatigue. Fracture mechanics, although provides a consistent framework to model and simulate the crack propagation phase, especially with the advent of the J integral and EPFM, which takes into account crack tip plasticity, is beyond the scope of this research. The focus of this research is to use a novel approach for CDM simulation such that the crack initiation phase can be predicted. A brief review of CDM especially in the context of fatigue loading is presented in the next chapter.

Chapter 2

Continuum damage mechanics

This chapter inherently discusses about the basic theory and modelling of quasi-static evolution of structures satisfying the momentum balance, the energy balance, and the entropy balance. No new formulations are presented here but classical theories of modelling elasto-plastic processes involving damage as an internal variable are discussed, details of which can be found in [START_REF] Lemaitre | A course on damage mechanics[END_REF], [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. Lastly certain innovations based on damage mechanics to tackle fatigue processes are highlighted.

Continuum mechanics: an overview

To obtain the global admissibility conditions, it is necessary that a reference problem is introduced in a spatial domain Ω, as shown in fig. 2.1. The region Ω has a boundary ∂Ω, which is subdivided into ∂Ω 1 and ∂Ω 2 , such that ∂Ω = ∂Ω 1 ∪ ∂Ω 2 . The body is subjected to specified surface force per unit surface area F d on ∂Ω 2 and to specified displacement U d on ∂Ω 1 . The outward normal vector n is defined at any material point M on the boundary ∂Ω. The body is also subjected to body force per unit volume f d . The evolution of the structure is considered to be quasi-static within time t ∈ [0, T ], with T being the total time.

Ω F d U d f d ∂Ω 2 ∂Ω ∂Ω 1 Figure 2.1: Reference problem in domain Ω

Admissibility conditions

In this thesis, infinitesimal deformation is considered, with no distinction between Lagrangian and Eulerian description. Considering u ( x, t) is the displacement field of the body Ω, with x representing the position vector, the velocity field can be written as v ( x, t) = ˙ u ( x, t), and the corresponding strain tensor is ε. If the body is subjected to a virtual displacement field δ u, the total internal virtual power can be written as (see [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF])

W i = - Ω×[0,T ] σ : δ εdV dt, (2.1)
where δ ε is the virtual strain rate corresponding to the virtual velocity δ ˙ u such that

δ ε = 1 2 ∇ δ ˙ u + ∇δ ˙ u T , (2.2)
and σ is the Cauchy stress tensor. The external virtual power is given by

We = Ω×[0,T ] f d • δ ˙ u dV dt + ∂Ω 2 ×[0,T ] F d • δ ˙ u dSdt. (2.3)
From the virtual work principle, it can be established that

W i + We = 0, =⇒ - Ω×[0,T ] σ : δ εdV dt + Ω×[0,T ] f d • δ ˙ u dV dt + ∂Ω 2 ×[0,T ] F d • δ ˙ u dSdt = 0.
(2.4)

It has to be noted that δ u is kinematically admissible, i.e.

ε = 1 2 ∇ u + (∇ u) T , in Ω, (2.5a) u = U d , i.e. δ u = 0 on ∂Ω 1 . (2.5b)
Using these conditions, along with the divergence theorem, eq. (2.4) can be re-written as

Ω×[0,T ] ∇ • σ + f d • δ ˙ u dV dt + ∂Ω 2 ×[0,T ] F d -σ • n • δ ˙ u dSdt = 0, ∀δ ˙ u.
(2.6) Equation (2.6) or eq. (2.4) is known as the static admissibility condition. From eq. (2.6), the static equilibrium equation and the static boundary condition are obtained as

∇ • σ + f d = 0, in Ω, (2.7a
)

F d = σ • n, on ∂Ω 2 .
(2.7b)

In the previous case, the variation principle was applied to the rates of displacements. If similar formulations are performed on stress, complementary virtual work can be obtained (see [START_REF] Wunderlich | Mechanics of Structures: Variational and Computational Methods[END_REF]. The principle of complementary virtual work is dual to the principle of virtual work. If Ω is subjected to a virtual stress field δσ, the internal complementary virtual power can be written as

W c i = - Ω×[0,T ]
δσ : ε dV dt.

(2.8)

Considering a virtual force δ Fu is acting on the boundary ∂Ω 1 such that δ Fu = δσ • n, on ∂Ω 1 , (2.9) the external complementary virtual power can be written as,

W c e = ∂Ω 1 ×[0,T ] δ Fu • ˙ u dSdt = ∂Ω 1 ×[0,T ] δσ • n • ˙ U d dSdt.
(2.10) Similar to the previous case, the principal of complementary virtual work states,

W c i + W c e = 0, =⇒ - Ω×[0,T ] δσ : ε dV dt + ∂Ω 1 ×[0,T ] δσ • n • ˙ U d dSdt = 0. (2.11)
It has to be noted at this point, that δσ has to be statically admissible, i.e.

δ (∇ • σ) = 0, in Ω, (2.12a) σ. n = F d i.e. δσ • n = 0 on ∂Ω 2 .
(2.12b) Equation (2.11) is called the kinematic admissibility condition (see [START_REF] Wunderlich | Mechanics of Structures: Variational and Computational Methods[END_REF]. Similar to the static admissibility condition, it is possible to obtain the corresponding strong form, i.e. the strain-displacement relationship and the kinematic boundary condition from eq. (2.11). These can be written as,

ε = 1 2 ∇ u + (∇ u) T , in Ω, (2.13a) u = U d , on ∂Ω 1 . (2.13b)

Constitutive relations

Apart from the admissibility conditions described before, the structure also needs certain local laws which describe the material properties and are called constitutive relations. The starting point of any constitutive relation is the first law of thermodynamics, or the conservation of energy. Without going into the details of thermodynamics, it is possible to write the first law in the form,

ρ ė = σ : ε + w -∇ • q, (2.14)
where ρ is the mass density of the material, e is the specific internal energy, w is the volumetric density of the internal heat production, and q is the heat flux vector.

The next is the second law of thermodynamics, which states that the total entropy production rate is always greater than or equal to the rate of heating divided by the temperature T . This can be written as

ρ ds dt - w T + ∇ • q T ≥ 0, (2.15)
where s is the specific entropy. Equation (2.14) can be used to rewrite eq. (2.15) as

ρ T ds dt - de dt + σ : ε -q • ∇T T ≥ 0.
(2.16)

The next step is to define Helmholtz specific free energy Ψ = e -T s, to rewrite eq. (2.16) as

σ : ε -ρ Ψ + s Ṫ -q • ∇T T ≥ 0.
(2.17)

This equation is called the Clausius-Duhem inequality. It has to be considered now, that the thermodynamic state of a material medium at a given point and instant is completely defined by the knowledge of certain variables at that particular point and instant. These variables are called the state variables, and the physical phenomena they describe will be thermodynamically admissible if the Clausius-Duhem inequality is satisfied. Certain state variables can be measured directly and are called observable variables, e.g. the total strain ε and the temperature T . Also, for dissipative phenomena, the current state depends on the past history, which is represented by the values at each instant of certain other variables called internal variables. Restricting to infinitesimal strain theories, the strain can be decomposed additively, into two internal variables,

ε = ε e + ε p . (2.18)
where ε p is the plastic or permanent strain, and ε e is the elastic strain. To describe dissipative phenomena other than plasticity, e.g. damage or hardening, other internal variables have to be introduced. The choice of the type or number of internal variables is not restricted. However, they are dictated by the physics of the phenomena that these variables are intended to describe. For the time being, let Υ i , i ∈ Z, be the set of all internal variables used for the following formulation. The concepts of plasticity along with other internal variables that dictates the dissipative phenomenon will be dealt later on. The general choice of Helmholtz free energy should be such that the Clausius-Duhem inequality can be rewritten as

Ψ = Ψ (ε -ε p , T , Υ i ) = Ψ (ε e , T , Υ i ) , ( 2
σ -ρ ∂Ψ ∂ε e : εe + σ : εp -ρ s + ∂Ψ ∂T Ṫ -ρ ∂Ψ ∂Υ i Υi -q • ∇T T ≥ 0. (2.22)
Considering uniform and constant temperature conditions (no thermal effects), the Clausius-Duhem inequality reduces to

σ -ρ ∂Ψ ∂ε e : εe + σ : εp -ρ ∂Ψ ∂Υ i Υi ≥ 0. (2.23)
As eq. (2.23) is valid for any εe , the elastic state law can be written as

σ = ρ ∂Ψ ∂ε e .
(2.24) Equation (2.20) along with eq. ( 2.24) can be combined to give

σ = ρ ∂Ψ ∂ε e = ρ ∂Ψ ∂ε = -ρ ∂Ψ ∂ε p . (2.25)
This defines the stress tensor σ as the associated or conjugate variable of ε e , ε, and -ε p . In a similar manner, thermodynamic forces associated with the set of the internal variables Υ i can be defined as

Λ i = ρ ∂Ψ ∂Υ i . (2.26)
Using the state laws (eq. ( 2.24) and eq. (2.26)), the Clausius-Duhem inequality (eq. (2.23)) can be rewritten as

σ : εp -Λ i Υi ≥ 0, (2.27)
which basically states that the internal intrinsic dissipation is positive. In order to define the laws related to the dissipation process, a dissipation potential is defined of the form

Φ = Φ εp , Υi . (2.28)
This function Φ is a positive convex function with zero at the origin of the space of the flux variables εp , Υi (see de Souza Neto et al., 2011, Lemaitre and[START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]. The complementary laws are thereby defined as

σ = ∂Φ ∂ εp , (2.29a) Λ i = - ∂Φ ∂ Υi . (2.29b)
However, it is more convenient to express the complementary laws as functions of the dual variables. To achieve that, a dissipation pseudo-potential Φ * (σ, Λ i ) is obtained through Legendre-Fenchel transformation of the dissipation potential Φ εp , Υi . This pseudo-potential Φ * , being the dual of Φ, helps to write the complementary laws in the form

εp = ∂Φ * ∂σ , (2.30a) Υi = - ∂Φ * ∂Λ i . (2.30b)
These equations are called the evolution equations of the flux variables as functions of the dual variables.

It can be seen in eq. (2.30) that the rates of change of the flux variables are normal to the surface Φ * . This normality rule is obeyed by all generalised standard materials. The only difficulty hereby is to define the pseudo-potential Φ * . For behaviour independent of velocity, Φ * becomes non-differentiable, and the normality rule can then be written in terms of another potential function F as

εp = λ ∂F ∂σ , (2.31a) Υi = - λ ∂F ∂Λ i , (2.31b)
where λ is the Lagrange multiplier which is obtained from the consistency condition eq. (2.33). For associative flow rules, the potential F is taken to be the plastic yield function f p , and in such cases the loading-unloading criterion is given by the Karush-Kuhn-Tucker condition λ ≥ 0, f p ≤ 0, λf p = 0.

(2.32)

The consistency condition is given by λ ḟ p = 0.

(2.33)

General constitutive behaviour

The most primitive form of constitutive behaviour is the linear elastic state law. To define a linear elastic state law, the elastic free energy function Ψ e is written as,

ρΨ e = 1 2 ε e : Cε e , (2.34)
where C is called Hooke tensor or elasticity tensor. This helps to rewrite eq. ( 2.24) as σ = Cε e .

(2.35)

However, such phenomenon is completely reversible and is valid till there are no permanent deformations of the structure. Certain materials, especially brittle materials, can be completely quantified using elastic stress-strain relationships, which fail without any permanent deformations. For ductile materials, e.g. metals, however, there are appreciable permanent deformations before rupture and the material behaviour are non-linear and dissipative after the elastic limit is reached. Such behaviour is essentially deemed to be plastic.

Concept of plasticity

The most idealistic plastic behaviour is perfect plasticity, which rarely happens for metals in practical scenario. Perfect plasticity or ideal plasticity results in infinite deformations after the elastic limit is reached, or in other words the load carrying capacity is completely lost after the elastic limit is crossed. Strain hardening, however in general, results in a non-linear increase in stress with respect to the strain after yielding. Due to hardening, the material does not loose the load carrying capacity at the onset of yielding. The stress-strain relationship is complicated and can in the most simplified case be defined in terms of hardening modulus which is a function of total or plastic strain, or can be described as power law, e. g. the classical case of Ramberg-Osgood equation (see [START_REF] Ramberg | Description of stress-strain curves by three parameters[END_REF]. Perfect plasticity and strain hardening are illustrated in fig. 2.2.

However, Ramberg-Osgood type equations do not capture the stress-strain behaviour accurately, especially for cyclic plasticity. For an accurate description of hardening, it is more convenient to use separate internal variables. Classically, hardening can be classified into kinematic and isotropic. For isotropic hardening, the radius of the yield surface (the surface of elastic limit) in the principal stress space increases. In this case the yielding during tension and compression is considered to be the same. This kind of hardening fails to take into account Bauschinger effect, which indicates the loss of isotropy and lower strength during compression, prevalent mostly in metals. Thereby, kinematic hardening is introduced, where the yield surface only translates in the principal stress space without any change in radius. These phenomena are shown in fig. 2.3. For metals, in general, either pure kinematic hardening or mixed hardening model is used. It has to be noted that reverse phenomena called softening also exist, which have been detailed in [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF].

Rate-independent plasticity

As the name suggests rate-independent plasticity or often referred as just "plasticity", is confined to those cases where the rate of loading does not have any influence. Before beginning any description of constitutive relations it is necessary to define the accumulated or cumulative plastic strain p as a nondecreasing quantity, i.e.

ṗ = 2 3 εp : εp 1/2 , =⇒ p = t 0 2 3 εp : εp 1/2 dτ.
(2.36)

Now, the free energy function mentioned in eq. (2.19), can be decoupled into parts, one with respect to elasticity and the other one with respect to hardening,

Ψ = Ψ e (ε e ) + Ψ p (X) , (2.37)
where Ψ e is the contribution from elasticity and Ψ p gives the contribution from the set of hardening variables X = {α, r}, with α being the flux variable for kinematic hardening and r being the flux variable for isotropic hardening. It is convenient at this point to use the Gibbs specific free enthalpy Ψ * , by a partial Legendre transformation of Ψ on the strain (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. This gives

Ψ * = Ψ * ,e (σ) + 1 ρ σ : ε p -Ψ p (X) .
(2.38)

The elastic state law can be obtained then by differentiating the Gibbs free enthalpy with respect to the associated variable (stress tensor), i.e.

ε = ρ ∂Ψ * ∂σ = ρ ∂Ψ * ,e ∂σ + ε p , =⇒ ε e = ρ ∂Ψ * ,e ∂σ .
(2.39)

The other state laws for the hardening variables are obtained as

Z = -ρ ∂Ψ * ∂X = ρ ∂Ψ p ∂X , (2.40)
with Z being the set of thermodynamic forces, Z = {β, R}, which are conjugate to X. In the most classical case, Ψ * ,e is defined as

ρΨ * ,e = 1 2 σ : C -1 σ. (2.41)
This gives the elastic state equation to be

ε e = C -1 σ.
(2.42)

As far as the state laws for the hardening variables are concerned, mostly, the kinematic relation is taken to be linear and the isotropic relation is considered to be exponential. This is achieved by using Ψ p as

ρΨ p = 1 2 α : Qα + g (r) (2.43)
with Q being tensor of material parameters similar to the Hooke tensor and g being a function depending on r such that

g ′ (r) = R∞ (1 -exp (-γr)) , (2.44)
where R∞ and γ are the material parameters. The tensor Q in Voigt matrix notation will be a diagonal matrix with equal diagonal terms Q. This would give the equations of state as

β = Qα, (2.45a) R = R∞ (1 -exp (-γr)) .
(2.45b) Equation (2.45a) shows a linear kinematic state law and eq. (2.45b) shows a non-linear (exponential) isotropic state law (see [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF]. These are commonly used formulations, however either non-linear kinematic law or linear isotropic law can also be used (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. The choice of the type of law just has to be consistent with the real physical behaviour of the material it mimics. At this point it is necessary to introduce what is called the plastic yield function f p . This is basically a surface that demarcates the elastic state and the plastic flow. The most classical form of yield function that is used for metals is of the von Mises type,

f p = 3J 2 (σ -β) -R -σy, (2.46)
with σy being the initial yield stress and J 2 being the second invariant of the deviatoric part of σβ.

Considering, β to be completely deviatoric, eq. ( 2.46) can be re-written as

f p = 3 2 [(σ D -β) : (σ D -β)] -R -σy, (2.47)
where σ D is the deviatoric stress tensor. The idea of plasticity is that this function cannot be greater than zero. If f p < 0, it is considered to be in the elastic regime; if f p = 0 and ḟ p = 0, plastic flow is said to occur; and if ḟ p < 0 with f p = 0, it is considered to be elastic unloading. The evolution equations (eq. (2.31)) can then be written as

εp = 3 2 σ D -β 3 2 [(σ D -β) : (σ D -β)] λ, (2.48a) α = 3 2 σ D -β 3 2 [(σ D -β) : (σ D -β)] λ, (2.48b) ṙ = λ. (2.48c)
It is evident from eqs. (2.48a) to (2.48b) that for von Mises associative plasticity, the plastic strain and the back strain are the same, and the plastic multiplier and the isotropic internal variable are the same. Equation (2.48a) can be used to rewrite eq. ( 2.36) as ṗ = λ, (2.49) which basically gives ṙ = ṗ = λ.

(2.50)

Remaining challenge is to calculate the plastic multiplier λ, such that λ = 0 when f p < 0 or ḟ p < 0; and λ = 0 when f p = 0 or ḟ p = 0, for which λ can be calculated from the consistency condition ḟ p = 0.

Although, yield surface of von Mises type has been considered here, other forms such as Tresca yield surface, Mohr-Coulomb yield surface, Drucker-Prager yield surface are also possible (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. However for metals, von Mises yield function is the standard. Many variations of the von Mises type visco(plasticity) models exist, one of the most important is the Marquis-Chaboche visco(plasticity) model (see [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF] which includes an extra softening term in the yield function, i.e.

f p = 3 2 [(σ D -β) : (σ D -β)] + a 2 β : Q -1 β -R -σy.
(2.51)

Viscoplasticity

This is a more general case, where the evolution equations depend on the rate of loading. This is more profound in metals at high temperature (e.g. creep). In the case of viscoplasticity the plastic domain is given by f p > 0, which is of course not possible for plasticity; and the elastic domain is given by f p ≤ 0. Plasticity can be interpreted as a limiting case of viscoplasticity (see [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF]. For the formulation of viscoplasticity, the equations of state remain the same as given in eq. (2.42) and eq. (2.45). The only change that needs to be incorporated is in the evolution equations. In the case of visco-plasticity, the pseudo-dissipation potential Φ * introduced in eq. ( 2.30) has a closed form and is differentiable. The most common is to introduce a pseudo potential of the form

Φ * = k n + 1 f p n+1 + , (2.52)
where k and n are viscosity parameters. Thereby, eq. (2.48a) can be re-written as

εp = k f p n +   3 2 σ D -β 3 2 [(σ D -β) : (σ D -β)]   .
(2.53)

In a similar manner, the evolution equation for the hardening variable can also be written as

α = k f p n +   3 2 σ D -β 3 2 [(σ D -β) : (σ D -β)]   , (2.54a) ṙ = k f p n + .
(2.54b)

Comparing these equations with eq. ( 2.48), the term k f p n + can be interpreted as the viscoplastic multiplier and according to the definition of + this viscoplastic multiplier λ = k f p n + > 0 when f p > 0 and zero otherwise. Equations (2.53) to (2.54) are known as Norton's law, and although there exist a lot of other viscoplastic laws (see [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF], it remains one of the most popular.

Concept of damage

The concept of damage as a continuum variable was first introduced by Kachanov (see [START_REF] Kachanov | In time to rupture in creep conditions (in Russian)[END_REF], considering a representative volume element (RVE) on which all properties were represented by homogenised variables. Damage, especially for metals, is always related to plastic strain or plastic dissipation, either at mesoscale (scale of RVE) or at microscale (scale of discontinuity) (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. As damage can be interpreted as the growth of micro-cavities or micro-cracks, the intrinsic variable was defined by Kachanov as

D = δS D δS , (2.55)
where δS is the total area of cross section of an RVE, and δS D is the total area of micro-voids in the cutting plane of the RVE which basically represents the surface density of microvoids (see fig. 2.4). In this thesis the damage is considered as isotropic phenomenon described by a scalar variable. For information regarding tensorial damage variables, [START_REF] Lemaitre | A course on damage mechanics[END_REF], [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF], [START_REF] Murakami | Continuum Damage Mechanics: A Continuum Mechanics Approach to the Analysis of Damage and Fracture[END_REF] are suggested. The internal variable D introduces a new quantity of interest, called the effective stress σ, which is defined as

σ = σ 1 -D .
(2.56) Equation (2.56) shows the increase in effective stress because of damage. The effect of damage can also be interpreted as the decrease in the modulus of elasticity. The effective modulus of elasticity Ẽ is defined as

Ẽ = E (1 -D) .
(2.57)

It is obvious that the scalar variable can take values between 0 and 1, with 0 being the virgin state and 1 indicates complete failure.

Damage with elasto-(visco)plasticity

The concept of damage being already formulated, the only thing remaining is to couple the scalar damage variable D with elasto-(visco)plasticity.

The first is to reformulate eq. ( 2.41), as

ρΨ * ,e = 1 + ν 2E σ ij σ ij 1 -D - ν 2E σ 2 kk 1 -D . (2.58)
This gives the elastic state law as

ε ij = 1 + ν E σ ij 1 -D - ν E σ kk 1 -D δ ij . (2.59)
This also introduces a new variable Y which is the thermodynamic force conjugate to D. This variable can be obtained as

Y = ρ ∂Ψ * ,e ∂D = 1 + ν 2E σ ij σ ij (1 -D) 2 - ν 2E σ 2 kk (1 -D) 2 . (2.60)
For a loading that is cyclic, the damage law has to be modified. This is due to what is called "quasiunilateral conditions of microdefects closure" (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF], resulting in different material behaviour in tension and compression (see [START_REF] Ladevèze | Damage effective stress in quasi unilateral conditions[END_REF]. The difference in the elastic state equation during tension and compression is due to the fact that during compression the material regains some stiffness, as some of the micro-defects are closed (fig. 2.5), thereby increasing the effective area. This can be represented by the effective modulus of elasticity during tension Ẽ+ , which σ σ σ σ

Ẽ+ = E (1 -D)
, where E is the true modulus of elasticity, and by the effective elastic modulus during compression Ẽgiven by Ẽ-= E (1 -hD). The closure parameter h has values between 0 (complete stiffness recovery) and 1 (no stiffness recovery). This phenomenon is described by modifying eq. ( 2.58) such that

ρΨ * ,e = 1 + ν 2E σ + ij σ + ij 1 -D + σ - ij σ - ij 1 -hD - ν 2E σ kk 2 1 -D + -σ kk 2 1 -hD .
(2.61) Equations (2.59) to (2.60), therefore become

ε e ij = 1 + ν E σ + ij 1 -D + σ - ij 1 -hD - ν E σ kk 1 -D + -σ kk 1 -hD δ ij , (2.62)
and

Y = 1 + ν 2E σ + ij σ + ij (1 -D) 2 + h σ - ij σ - ij (1 -hD) 2 - ν 2E σ kk 2 (1 -D) 2 + h -σ kk 2 (1 -hD) 2 , (2.63)
respectively. The term δ ij in eq. ( 2.62) is the Kronecker delta, i.e. δ ij = 1 for i = j and 0 otherwise. It should be mentioned that the decomposition of the stress tensor σ ij into a positive part σ + ij and a negative part σ - ij is based on the information of the principal stresses and principal directions. The positive and negative principal stresses along with their corresponding normalised principal directions are used to build σ + ij and σ - ij respectively. Thereby, σ kk and -σ kk can be interpreted as traces of σ + ij and σ - ij respectively. Detailed information can be found in [START_REF] Lemaitre | A course on damage mechanics[END_REF], [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF].

The next is to update the evolution equation. Damage has a significant influence on the yield function f p previously defined by eq. (2.47), which is henceforth modified as

f p = 3 2 σ D 1 -D -β : σ D 1 -D -β -R -σy. (2.64)
To incorporate damage in the evolution equations, non-associative flow rules have to be considered. The potential F mentioned in eq. (2.31) is now written as

F = f p + F D , (2.65)
where F D is a non-linear function of Y . The most classical choice is to use (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF])

F D = S (1 + s) (1 -D) Y S s+1 ,
(2.66)

with s and S being material parameters. Equation (2.48) can then be re-written as

εp = 3 2 σ D 1 -D -β 3 2 σ D 1 -D -β : σ D 1 -D -β λ 1 -D , (2.67a) α = 3 2 σ D 1 -D -β 3 2 σ D 1 -D -β : σ D 1 -D -β λ, (2.67b) ṙ = λ, (2.67c)
where λ is the plastic or viscoplastic multiplier. Equation (2.50) becomes

ṙ = λ = ṗ (1 -D) .
(2.68)

The rate of growth of damage is then given by

Ḋ = ∂F D ∂Y = Y S s ṗ.
(2.69)

The threshold of damage for damage law of the form of eq. ( 2.69) is generally given by a certain level of cumulative plastic strain p D (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. Thereby, Ḋ is given by eq. ( 2.69) if p > p D and zero otherwise. The damage will continue to evolve till a critical damage level Dc is achieved, which marks the initiation of a macro-crack.

Apart from the particular damage model which has been illustrated here, several other damage models exist. The choice of the damage model depends not only on the material (e.g. ductile, brittle), but also on the type of loading (e.g. fatigue loading, creep loading). For different types of damage models, [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF] is suggested.

Another damage law formulation similar to Bellenger andBussy (1998, 2000) can also been proposed. A viscous pseudo damage potential is proposed which is of the form

Φ * ,d = k d n d + 1 f d n d +1 + , (2.70)
where k d and n d are material parameters, f d can be interpreted as damage function and in the most simplistic case is given by

f d = Y -Y th , (2.71)
where Y th describes the damage threshold. The evolution equation for damage is then given by

Ḋ = ∂φ * ,d ∂Y = k d f d n d + .
(2.72)

Although damage mechanics is able to model the initiation point of a macro-crack, the drawback of this is the localisation of the quantities of interest. The solution is in most cases mesh dependent, and thereby localisation limiters are necessary. The earliest theories on this are to use constant mesh size or even use mesh size as a parameter. However, later on, a material internal length was introduced in the constitutive relations. This basically boils down to the fact that localised quantities, like plasticity or damage, for a material point also depend on a small domain around it (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. Therefore, non-local theories are proposed. One of the earliest ideas is to replace each variable by the corresponding averaged homogenised variable considering a small domain surrounding each material point (see [START_REF] Bažant | Measurement of characteristic length of nonlocal continuum[END_REF]. Another, alternate way is to introduce gradient theories (see [START_REF] Aifantis | The physics of plastic deformation[END_REF], where plasticity and damage at any material point depend on the corresponding thermodynamic variables and also on the first n derivatives in the Taylor expansion series. However using non-local or gradient-enhanced damage models are beyond the scope of the thesis and only localised damage parameters are used.

Although classical damage formulations encompass the damage variable between 0 and 1, there have been other formulations, for instance, in Bellenger andBussy (1998, 2000), where a damage variable ranging between 0 and ∞ has been used. However, in this thesis, classical definition of D between 0 and 1 will be used.

Another type of damage model which has been popular off late, is the Gurson-Tvergaard-Needleman (GTN) model. This type of model gives a physically improved formulation of micro-void nucleation, growth, and coalescence process. The damage variable in this case, defined by the void volume fraction, is considered to be evolving separately during void nucleation and void growth (see [START_REF] Acharyya | A complete GTN model for prediction of ductile failure of pipe[END_REF]. Although, this gives an improved representation of the micro-mechanical process, it is not free from localisation problem. Hence, localisation limiters are used, e.g. [START_REF] Fish | A nonlocal multiscale fatigue model[END_REF] uses a nonlocal GTN model. Although GTN model is beyond the scope of this thesis, it is worthwhile to mention it.

CDM approaches in fatigue

Different methodologies have been been adopted over the years to model as well as to simulate the fatigue mechanisms using continuum damage mechanics. A few most important approaches are presented here.

One of the most novel approaches used is the formulation of the "two-scale damage model" by Lemaitre (see [START_REF] Lemaitre | A two scale damage concept applied to fatigue[END_REF]. This is mainly used for high cycle fatigue. The global behaviour remains elastic, which is obtained through a macroscopic elastic structural calculation. After that the Gauss points with maximum stress concentrations are identified. Eshelby-Kroner scale transition law is then applied to get the micro-strain at chosen Gauss points from the macro-strain. Thereafter, elasto-plastic analysis coupled with damage is performed at the micro-scale. The final outcome of this is the damage variable D, which is considered to be equal in both the scales.

An important development in the solution framework of fatigue is the popular "jump cycle" approach (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. The basic idea here is confined to loads that are periodic or at least periodic by blocks, where calculations of full blocks of load cycles are skipped. The quantities of interest are calculated for a certain number of initial load cycles. The damage and accumulated plastic strain increments over the last computed cycle are then estimated. From this information, an estimation of the number of cycles that can be skipped is made with an assumption that damage and accumulated plastic strain are linear with respect to the number of cycles. After that the quantities of interest are again calculated for one cycle and an estimation is made for the number of cycles to skip. This continues till Dc is achieved. The "two-scale damage model" is generally used together with the "jump cycle" technique, e.g. in [START_REF] Bhamare | A multi-temporal scale approach to high cycle fatigue simulation[END_REF], [START_REF] Lemaitre | A two scale damage concept applied to fatigue[END_REF].

Another important approach is the time homogenisation technique initially proposed in [START_REF] Guennouni | Réponse homogénéisée en temps de structures sous chargements cycliques[END_REF] for quasi-static elasto-viscoplastic problems. Later on Puel and Aubry extended this method for dynamics (see [START_REF] Puel | Material fatigue simulation using a periodic time-homogenisation method[END_REF] and also for CDM (see [START_REF] Devulder | Two-time scale fatigue modelling: application to damage[END_REF]. Fish and his co-workers, also proposed similar techniques for fatigue, using a viscoplasticity like damage model (see Oskay and Fish, 2004b) and also for GTN damage model (see [START_REF] Fish | A nonlocal multiscale fatigue model[END_REF]. The main idea is to divide the time into two separate scales: a fast scale (micro-time scale) and a slow scale (macro-time scale). Thereafter, the quantities at each slow time scale are represented as homogenised quantities of the fast time scale. This solves the problem of numerical simulation of all the quantities of interest at the fast time.

Concluding remark

A review of classical material non-linearities in a continuum framework has been presented in this chapter, with damage being the primary quantity of interest. Although, a few methods of computing large number of cycles have been established, a novel approach using LATIN method and model reduction technique will be used. An introductory review of classical LATIN approach and different model reduction techniques will be given in the next chapter.

Chapter 3

Reduced order modelling and large time increment method High degrees of fidelity in complex coupled multi-physics or multi-scale problems may result in large run-times and extreme memory consumption. Even with the availability of high performance computing platforms, such computations, in many industrial or academic cases can be rendered infeasible. Therefore, it is necessary to reduce these high-fidelity models, such that both memory-cost and run-time are reduced (see [START_REF] Cueto | Model order reduction based on proper orthogonal decomposition[END_REF]. The goal of model order reduction (MOR) methods or model reduction techniques is to extensively reduce computational cost and/or storage memory with a minimal loss of accuracy.

This chapter gives a brief overview of model reduction techniques for non-linear problems. The large time increment method in the light of material non-linearities is also introduced along with a succinct comparison with Newton-Raphson technique.

Classical incremental method

The classical or more popular approach to solve non-linear solid mechanics problems is to split the considered time domain [0, T ] into sufficiently small time intervals and then solve sequentially the boundary value problem at each time step. The equilibrium equation defined by eq. (2.6) for an incremental case can be written as a boundary value problem, which, after standard finite element discretisation, reduces to finding the nodal unknowns u n+1 such that the equilibrium equation represented by

r ( u n+1 ) ≡ F int ( u n+1 ) -F ext n+1 = 0, (3.1)
is satisfied, where F int denotes the nodal forces arising from internal stresses, F ext is obtained from the external applied body or/and surface forces, and r ( u n+1 ) is the residual vector (see [START_REF] De Souza Neto | Computational Methods for Plasticity: Theory and Applications[END_REF]. Newton-Raphson method is a particularly attractive solution for this kind of equation. Each iteration of the Newton-Raphson method consists in solving the linearised version of eq. (3.1). Considering u k-1 n+1 being known at the (k -1) th iteration of the Newton-Raphson method, the goal at iteration k is to solve the linear system given by

K T δ u k = -r k-1 , (3.2)
for δ u k , where K T is the global tangent stiffness matrix at (k -1) th iteration and the residual r k-1 is obtained as

r k-1 = F int u k-1 n+1 -F ext n+1 . (3.3)
The nodal displacement u k n+1 is then given as

u k n+1 = u k-1 n+1 + δ u k . (3.4)
Generally the Newton-Raphson procedure is continued till the residual r m at iteration m is low enough and satisfies a given convergence criterion. The schematic of the Newton-Raphson method is given in fig. 3.1. There are many varieties of Newton-Raphson methods which are the so-called modified Newton-Raphson methods. Some popular versions use the same tangent stiffness matrix for all the time increments, or use same stiffness matrix within each time increment and update only at the beginning of each increment, or update the stiffness matrix after certain iterations within each increment. For more details of classical incremental methods especially for material non-linearities, de Souza Neto et al. (2011), [START_REF] Reddy | An Introduction to Nonlinear Finite Element Analysis[END_REF], [START_REF] Wriggers | Nonlinear Finite Element Methods[END_REF] are some of the recommended books. High fidelity-problems with large number of degrees of freedom and encompassing large number of time steps, invariably lead to extremely high numerical expense and memory consumption hence it is a practical necessity to reduce the high-fidelity problem.

F int u F ext n+1 F ext n u 0 n+1 = u n u 1 n+1 u 2 n+1 u n+1 • • • δ u 2 δ u 1 r 1 r 2 K T K T

Model reduction techniques

Model reduction techniques can be classified into a posteriori model reduction techniques where the reduced solution is obtained after a similar high-fidelity solution (known as the training stage) is computed, and a priori methods which compute the approximation without knowing any high-fidelity solution (see Chinesta et al., 2014b). Whatever the type of model reduction technique may be, the purpose is to yield low-dimensional approximations for the full high-dimensional system (see [START_REF] Pinnau | Model reduction via proper orthogonal decomposition[END_REF].

For dynamic systems, especially for structural dynamics, the model reduction approaches are in general a posteriori as even though no high-fidelity pre-computations are needed, "basis vectors" are to be obtained before foregoing with any reduced computation. Some notable methods in this field are:

-Modal displacement method, where modal shape vectors are used as reduced basis. Eigenvalue analysis is performed for the structure in question and very few eigenvectors are extracted to build the reduced basis, on which the equation of motion is projected. Modal acceleration and modal truncation are varieties of modal superposition methods.

-Krylov subspace method, which reduces a system with higher degrees of freedom to a system with lower degrees of freedom but with similar input-output behaviour. A transfer function is calculated in the Laplacian domain from the time-invariant state-space system. Two bi-orthogonal projection matrices are calculated from the information of the moments of the transfer function, which are then used to reduce the original state-space system.

For a review of model reduction techniques in dynamics, [START_REF] Besselink | A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control[END_REF] is suggested. In this thesis, however, the focus is on systems exhibiting material non-linearities.

Proper orthogonal decomposition (POD): a posteriori model reduction technique

The basic idea of POD is to compute certain full-order problems in order to extract relevant information which can then be used to calculate similar problems more efficiently. The first ideas of POD were developed in [START_REF] Karhunen | Zur Spektraltheorie stochastischer Prozesse[END_REF], [START_REF] Kosambi | Statistics in function space[END_REF], [START_REF] Loève | Fonctions aléatoires de second ordre[END_REF], [START_REF] Pearson | On lines and planes of closest to points in space[END_REF]. The usage of POD in the field of mechanics deals with the creation of POD basis, from the snapshots obtained from the solution of the training phase, and then to use this basis to solve the intended problem in a reduced space (see [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF]. Also known as principal component analysis (PCA), Karhunen Loève expansion (KLE) and singular value decomposition (SVD) (see [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF], POD provides an optimally ordered set of basis functions in a least square sense for the full-order solution. These basis functions are called "proper orthogonal modes", "empirical eigenfunctions", or just "basis vectors". A reduced order model (ROM) or a surrogate model can then be generated by truncating the optimal basis (see [START_REF] Liang | Proper orthogonal decomposition and its applications -part I: Theory[END_REF][START_REF] Pinnau | Model reduction via proper orthogonal decomposition[END_REF].

A full solution given by u ( x, t) where x denotes the spatial coordinate and t being the temporal coordinate can be represented as

u ( x, t) = r i=1 T i (t) X i ( x) (3.5)
such that the approximation becomes exact as r → ∞. The space functions X i ( x) are chosen such that they are orthonormal and the determination of the time functions T i (t) should only be depended on the corresponding space functions X i ( x).

For a discretised case, if U are values at n space and m time points, the solution U will be an n × m matrix. One way to obtain the reduced basis is to perform a singular value decomposition (SVD) of the solution U such that U = UΣV T ,

(3.6)
where U is an n × n orthogonal matrix, V is an m × m orthogonal matrix and Σ is an n × m diagonal matrix. The diagonal Σ ii consists of r = min (n, m) non-negative numbers in decreasing order which are called the singular values of U. To establish the equivalence of eq. (3.5) and eq. (3.6), it is assumed that Q = ΣV T , where Q is an n × m matrix. Then eq. (3.6) can be written as

U = UQ = n i=1 u i q i , (3.7)
where u i and q i are ith column and ith row of the matrices U and Q respectively. Generally this summation needs not be for all n but can be for s columns and rows of U and Q respectively, with s ≪ n. The value of s can be determined from the magnitude of the singular values. Then U obtained from eq. (3.7) can be approximated as

U ≈ s i=1 u i q i . (3.8)
u i can be interpreted as function of space and q i can be interpreted as function of time. It has to be noted that singular value decomposition is a generalisation of eigenvalue decomposition for non-square matrices. If U is a square symmetric positive definite matrix, then the singular values (diagonal elements of Σ) are the same as eigenvalues, and the matrix of singular vectors are the same (i.e U = V) and equal to the matrix of the eigenvectors (both right eigenvectors and left eigenvectors) (see [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF]. From eq. (3.8), the first s columns of the U matrix (the {u i } s i=1 vectors) are used to construct a truncated matrix Up called the projection matrix, where U is of size n × s. (3.9)

It is also advisable to have each column of U be distributed around zero, i.e the mean of each column should be subtracted from every element of that column (see [START_REF] Kerschen | The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: An overview[END_REF][START_REF] Liang | Proper orthogonal decomposition and its applications -part I: Theory[END_REF]. It is obvious that C is an n × n matrix. Then the projection matrix Up is obtained from the eigenvectors i=1 are the eigenvectors that serve as basis vectors (see [START_REF] Cueto | Model order reduction based on proper orthogonal decomposition[END_REF]. KLE on the other hand can be considered to be an extension of PCA for infinite-dimensional spaces (e.g. spaces of continuous time and/or space functions). Although, KLE is not discussed here, [START_REF] Liang | Proper orthogonal decomposition and its applications -part I: Theory[END_REF] Once the projection matrix is built, it is used to solve the actual problem of interest in a reduced order space. If the non-linear problem defined by eq. (3.1), then at a particular time step n + 1, the vector of nodal unknowns u n+1 being of size n × 1 is projected in the reduced space by

u n+1 = Up g n+1 , (3.10)
where g n+1 is the vector of the reduced degrees of freedom and is of size s × 1. This definition is then used in eq. (3.2), and the result is pre-multiplied by the transpose of the projection matrix,

U T p K T Upδ g = -U T p r, =⇒ K T,r δ g = -rr, (3.11)
where the reduced tangent stiffness matrix K T,r is defined as

K T,r = U T p K T Up, (3.12)
and the reduced RHS rr is defined as rr = U T p r.

(3.13)

It might be noted that the iteration index k in eq. (3.2) has been dropped in eq. (3.11) just for convenience.

The reduced system of equation is computationally much less expensive compared to the full system. The problem however of classical POD-Galerkin formulation is that even if the global solution is computed in a reduced space, the computational effort to solve the non-linear constitutive relation at every local integration point remains high. To put that in perspective, the local computations are independent of the reduced model and the reduced approximation has no influence on the numerical expense with regard to the estimation of the internal variables that influence the tangent stiffness matrix. For linear problems this is not an issue as the stiffness matrix does not change and the time saving is directly proportional to the degree of reduction of the original finite element (FE) subspace to POD subspace.

To build an efficient ROM for non-linear problems, the hyper-reduction method has been proposed in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF], which tries to address this issue by introducing a reduced integration domain. The reduced integration domain is defined by choosing only a few elements of the complete mesh and considering the related balance condition. The internal variables are then not computed on the whole domain but only at the reduced integration domain, which are then extrapolated. This method gives robust computational time savings, and has been successfully used to simulate damage in [START_REF] Ryckelynck | A robust adaptive model reduction method for damage simulations[END_REF].

Another way of solving the issue of non-linear problems is to use discrete empirical interpolation methods as described in [START_REF] Radermacher | Pod-based model reduction with empirical interpolation applied to nonlinear elasticity[END_REF]. The discrete empirical interpolation method is a discrete version of the classical empirical interpolation method used to solve parametrised partial differential equations (see e.g. [START_REF] Lass | Pod galerkin schemes for nonlinear elliptic-parabolic systems[END_REF]. The idea of empirical interpolation was used first in [START_REF] Barrault | An "empirical interpolation" method: Application to efficient reduced-basis discretization of partial differential equation[END_REF] in the context of partial differential equations, which was later modified as the discrete empirical interpolation method in [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF]. The main idea described in [START_REF] Radermacher | Pod-based model reduction with empirical interpolation applied to nonlinear elasticity[END_REF] is to split the internal force vector of eq. (3.1) into a linear and a non-linear part, where the linear part can be expressed with the help of a constant stiffness matrix. Two separate projection matrices are used, one for the nodal unknowns and one for the non-linear part of the internal force. The dependency of the total number of DOFs on the non-linear term is thereafter minimised by means of an empirical interpolation. This ultimately leads to the reduced tangent stiffness matrix that consists of a linear constant part (calculated only at the initial time step) and another non-linear part whose calculation time is reduced by empirical interpolation of the internal force vector.

Apart from these techniques that have introduced modification in the classical Galerkin-POD method to solve non-linear solid mechanics problems, there have been other methods to improve the usage of POD, such as the generalised empirical interpolation method proposed in [START_REF] Maday | A generalized empirical interpolation method: Application of reduced basis techniques to data assimilation[END_REF] for data assimilation, or the missing point estimation procedure used in [START_REF] Astrid | Reduction of process simulation models: a proper orthogonal decomposition approach[END_REF] to reduce PDEs and applied in the process of heat transfer and computational fluid mechanics.

The usage of POD depends on the goal in hand. One usage may be to build the spatial basis by considering the training stage solution over a short time span and then use the reduced basis to calculate solution over a long time span. Another usage may be to calculate the reduced basis for a full-order problem for entire space-time domain in the training stage and then use the reduced basis to solve similar problems involving slight change in material properties or boundary conditions (see [START_REF] Cueto | Model order reduction based on proper orthogonal decomposition[END_REF].

Proper generalised decomposition (PGD): a priori model reduction technique

It is desired, ideally to have a technique where a reduced order approximation can be built without relying on any training stage. One would then be able to assess the accuracy of the reduced order approximation, and enrich the reduced order basis if necessary. This leads to the proper generalised decomposition (PGD), which was first introduced as "radial loading approximation" in the context of LATIN method (see [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF][START_REF] Ladevèze | The large time increment method for the analyse of structures with nonlinear constitutive relation described by internal variables[END_REF], 1999, and other works of Ladevèze) based on separated representation of the quantities of interest. For instance, the desired solution field V ( x, t) is obtained in the separated form as (see Chinesta et al., 2014a)

V ( x, t) = r i=1 T i (t) X i ( x) . (3.14)
Here the number of terms r needed for this finite sum decomposition is not known a priori. The functions T i (t) and X i ( x) are constructed by successive enrichments. For a particular enrichment stage v + 1, the functions {T i } v i=1 and {X i } v i=1 being known from the previous steps, T v+1 and X v+1 are sought. This is achieved by invoking the weak form of the problem. Considering a case where the problem is defined by

L (V ) = f (t, x) , (3.15) within the domain Ω × [0, T ].
The right hand side f is the external input, and the operator L (V ) can be defined as

L (V ) = Lt (V ) + Lx (V ) , (3.16)
where, Lt =

∂ ∂t

and Lx = ∇ 2 . The weak form corresponding to eq. (3.15) can be written as

Ω×[0,T ] δV • (L (V ) -f (t, x)) dV dt = 0.
(3.17)

The separated form represented in eq. (3.14) can be used to decompose eq. (3.17) into a space problem,

Ω δX v+1 • a x • Lx (X v+1 ( x)) + b x • X v+1 ( x) + w x v+1 ( x) dV = 0, (3.18)
and into a time problem,

[0,T ] δT v+1 • a t • Lt (T v+1 (t)) + b t • T v+1 (t) + w t v+1 (t) dt = 0. (3.19)
Here, a x , b x and a t , b t are constant terms. The terms w x v+1 ( x), and w t v+1 (t) are obtained by integrating

f (t, x) and L v i=1 T i (t) • X i ( x) over [0, T ]
and Ω respectively.

The spatial problem can be solved by any discretisation technique for elliptic boundary value problems, and the temporal problem can be solved by any discretisation technique for initial value problems.

These two problems are solved in a staggered manner, i.e. X v+1 is calculated knowing T v+1 and vice versa, till convergence is reached.

Apart from space-time separation, there are other forms of coordinate separation, e.g. fully separated decomposition for a three dimensional problems (fully separated PGD description has generally been used to address homogenisation problems in [START_REF] Chinesta | Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization[END_REF][START_REF] Prulière | On the solution of the heat equation in very thin tapes[END_REF] V (3.20) or for parametric problems, the solution can be separated using the parameter c as a separate coordinate (see [START_REF] Ammar | Solving parametric complex fluids models in rheometric flows[END_REF] 2013) which gives a detailed review on parametric PGD. PGD has been used extensively to solve many non-linear problems (see Chinesta and Cueto, 2014), for instance in fluid mechanics, the Fokker-Planck equation was solved using PGD in [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: Part II: Transient simulation using space-time separated representations[END_REF]. PGD was also used to solve stochastic equation within the Brownian Configuration Field framework in [START_REF] Chinesta | On the reduction of stochastic kinetic theory models of complex fluids[END_REF].

(x, y, z) = r i=1 X i (x) Y i (y) Z i (z) ,
, i.e. V ( x, t, c) = r i=1 X i ( x) T i (t) C i (c) .
In the case of solid mechanics, PGD has been used generally in the context of LATIN method especially in LMT-Cachan, which was initially developed for small-displacement (infinitesimal strain) problems, and will be discussed in details in the next section. Extramural to LATIN method, space-parameter PGD separation for multi-scale problems can be found in El [START_REF] Halabi | A PGD-based multiscale formulation for non-linear solid mechanics under small deformations[END_REF]. Apart from that, twodimensional space-space PGD separation has been used in [START_REF] Metoui | The proper generalized decomposition for the simulation of delamination using cohesive zone model[END_REF] to simulate delamination using cohesive zone model, and in [START_REF] Vidal | Composite beam finite element based on the proper generalized decomposition[END_REF] to solve beam problems. These are just a few examples of PGD-based work on solid mechanics that are extraneous to LATIN-based PGD frameworks.

Large time increment (LATIN) method

Contrary to classical incremental technique, Ladevèze first proposed this new non-incremental method in the context of non-linear mechanical problems (see Ladevèze, 1985a,b) more than 30 years ago. The LATIN method being a non-incremental iterative technique, in contrary to the classical incremental technique (previously described in section 3.1), is not built on the notion that the interval [0, T ] needs to be partitioned into small steps. The method generally starts with an approximation of the quantities of interest for the whole spatial domain Ω and the whole temporal domain [0, T ]. Thereafter, at each iteration, the quantities of interest are improved for the whole spatial and temporal domain till a convergence is reached. In general LATIN method is based on three basic principles (see [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF].

-The first principle is the separation of difficulties. The governing equations are separated into two groups. A group of equations which are local in space, possibly non-linear; and a group of equations which are linear but possibly global in space.

-The second principle is a two-step iterative algorithm, where at each iteration the solution is constructed alternatively for the first group of equations and then for the second group of equations.

-The last principle is to solve the global problem defined on the entire space-time domain Ω × [0, T ] using PGD based reduced order approximation such that the numerical cost is drastically lowered.

A pictorial representation of the LATIN scheme is given in fig. 3.5. Considering only the case of infinitesimal strain elasto-(visco)plasticity, the set of the quantities of interest is denoted by s = εp , Ẋ, σ, Z . For the time being, it is considered that the set of hardening variables follows a "normal formulation", i.e.

Z = ΠX, (3.22)
with Π being a linear operator of material properties relating the internal hardening variables X with the corresponding thermodynamic forces Z.

εp , Ẋ [σ, Z] s 0 s 1/2 s 1 s 3/2 s i s i+1/2 s ex H + H - A d Γ Figure 3.5: Schematic representation of LATIN method

First principle: separation of difficulties

To implement the first principle of LATIN method, two manifolds are introduced, based on a subdivision of the set of equations that must be solved.

-The first is the space A d which belongs to the manifold of the admissibility conditions (eqs. (2.6) and (2.11)) and the linear state laws (eqs. (2.35) and (3.22)). -The second is the space Γ which belongs to the manifold of the non-linear evolution equations (eq. (2.31), with Υ i = X and Λ i = Z in eq. (2.31b)). The exact solution is given by the intersection of these two manifolds, i.e. sex = A d ∩ Γ.

(3.23)

Second principle: two-step algorithm

The algorithm is initialised by a solution s 0 ∈ A d on the whole space-time domain, which is based on an elastic solution and εp , Ẋ, Z are considered to be zero. The two-step algorithm is then initiated which includes a local step, an intermediate step which involves calculating ŝi+1/2 ∈ Γ, knowing the solution set s i ∈ A d . This basically involves the solution of the evolution equations. The transfer of information from s i to ŝi+1/2 occurs through what is called the search direction operator. For the local step this equation is given by

εp i+1/2 -εp i -Xi+1/2 -Ẋi + H + σi+1/2 -σ i Ẑi+1/2 -Z i = 0, (3.24)
where the search direction operator H + represents the direction of ascent. The next step known as the global step is to calculate s i+1 ∈ A d , knowing the solution set ŝi+1/2 ∈ Γ by solving the linear state laws and the admissibility conditions. Another search direction equation is introduced which transfers the information from ŝi+1/2 to s i+1 ,

εp i+1 -εp i+1/2 -Ẋi+1 -Xi+1/2 -H -σ i+1 -σi+1/2 Z i+1 -Ẑi+1/2 = 0, (3.25)
where the search direction operator H -signifies the direction of descent. This iterative algorithm continues till a convergence is reached. LATIN method inherently defines a LATIN error indicator which basically measures the distance between ŝi+1/2 and s i+1 . If this indicator is below a pre-defined tolerance, convergence is said to have reached. It has also to be noted that search direction operators H + and H -only affect the rate of convergence of the algorithm (see [START_REF] Néron | Proper generalized decomposition for multiscale and multiphysics problems[END_REF] but not the converged solution. Theoretically these operators can be any symmetric, positive definite operators that may vary during the iterations. For elasto-(visco)plastic problems, the direction of ascent is chosen as

H + -1 = 0, (3.26)
i.e. the ascent direction is vertical (see [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF]Ladevèze, 1993, Relun et al., 2013). For the direction of descent, the fastest way to convergence is the usage of an operator which is associated to the tangent space of Γ, which is for viscoplastic problems the second derivative of the pseudo-potential Φ * (see [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF]Ladevèze, 1993, Relun et al., 2013). There can be other search directions that can be employed, readers are advised to consult [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF] for more information on search direction operators.

Third principle: model reduction method

Finally, the third principle essentially introduced a separation of variables, which was known back then as "radial loading approximation" (see [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF]. The idea is to solve the linear global stage (linear as all the non-linearities are dealt in the local stage) using certain space-time separated approximations which lead to the solution of the global admissibilities in a reduced space and generation of PGD-based ROM.

Galerkin formulation of the global stage

During the initial years of LATIN method (see [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF]Ladevèze, 1993, Ladevèze, 1999, as for examples) the solution process of the linear global stage was based on Galerkin formulation. The quantities of interest are represented in terms of corrections, i.e.

∆s i+1 = s i+1 -s i .
(3.27)

The global stage then essentially involves solving for ∆σ i+1 , ∆Z i+1 such that

Ω×[0,T ] C -1 ∆ σi+1 : δσ + Π -1 ∆ Żi+1 : δZ dV dt+ Ω×[0,T ] H -∆σ i+1 ∆Z i+1 : δσ δZ dV dt = Ω×[0,T ] ∆ε ∆X : δσ δZ dV dt, (3.28)
with the term on the right hand side being

∆ε ∆X = H - σi+1/2 -σ i+1 Ẑi+1/2 -Z i+1 - εp i+1/2 -εp i+1 -Xi+1/2 -Ẋi+1 (3.29)
∀ δσ which is kinematically admissible for homogeneous equation and ∀ δZ (see [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF], for the origin of eq. (3.28)). The quantities of interest are thereby written in a separated variable form, i.e. Detailed formulation can be found in [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF].

∆σ i+1 ( x, t) = m j=1 λ σ j (t) σj ( x) , (3.30a) ∆Z i+1 ( x, t) = m ′ j=1 λ Z j (t) Zj ( x) , ( 3

Minimisation formulation of the global stage

Apart from this Galerkin formulation, recently an alternative formulation has been developed to solve the global problem (see [START_REF] Ladevèze | PGD in linear and nonlinear computational solid mechanics[END_REF], Relun et al., 2011). This involves finding the solution set ∆s i+1 which minimises over A d the constitutive relation error e 2 CE . The constitutive relation error is associated with the search direction equation and is defined as

e 2 CE = ∆ εp i+1 -∆ Ẋi+1 -H -∆σ i+1 ∆Z i+1 + ∆ε ∆X 2 M , (3.31)
where the norm is defined as

2 M = Ω×[0,T ]
: M dV dt, (3.32)

and the operator M must be positive definite. The solution set ∆s i+1 is then calculated such that

∆s i+1 = arg min ∆s i+1 ∈A d ∆ εp i+1 -∆ Ẋi+1 -H -∆σ i+1 ∆Z i+1 + ∆ε ∆X 2 M . (3.33)
The quantities of interest are also written in a separated variable form, i.e. along with eq. (3.30),

∆ εp i+1 ( x, t) = m j=1 λ p j (t) εp j ( x) , (3.34a) ∆ Ẋi+1 ( x, t) = m ′ j=1 λ X j (t) Xj ( x) . (3.34b)
It has to be noted that along with the minimisation problem, the admissibility conditions should also be satisfied in a weak sense (as constraints to the minimisation weak form). Ultimately, separate time and space equations are obtained which are solved iteratively. For detailed solution of the minimisation problem, [START_REF] Néron | Proper generalized decomposition for multiscale and multiphysics problems[END_REF], [START_REF] Relun | Multiparametric strategy for robust design in fatigue[END_REF] are suggested. Generally, in the global stage of LATIN method, a combination of POD and PGD is used. At the beginning the reduced order basis is reused and only the time functions are updated, if the reduced order model hence formed is not satisfactory enough, the reduced basis is enriched by adding a space-time pair using the PGD technique (see [START_REF] Ladevèze | PGD in linear and nonlinear computational solid mechanics[END_REF]. At each LATIN iteration, generally the reduced basis is enriched by a maximum of one space-time mode (see [START_REF] Néron | Proper generalized decomposition for multiscale and multiphysics problems[END_REF].

A note on "normal formulation"

At the beginning of section 3.3, it was considered that the hardening variables follow a "normal formulation" given by eq. (3.22). Indeed, the efficient usage of LATIN algorithm requires that the state equations are defined by linear operators (see [START_REF] Bhattacharyya | A LATIN-based model reduction approach for the simulation of cycling damage[END_REF][START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF][START_REF] Ladevèze | The large time increment method for the analyse of structures with nonlinear constitutive relation described by internal variables[END_REF], 1999). The elastic state law (eq. (2.35)) and the kinematic state equation (eq. (2.45a)) introduced in chapter 2 already follow "normal formulations", however the isotropic state equation (eq. (2.45b)) does not. Henceforth, eq. (2.45b) is transformed into linear relations by a change of variable

R = R∞ r (3.35)
where, R and r are the new isotropic variables, such that

r = r 0 ∂R ∂r 1 R∞ 1/2 dr, =⇒ r = 2 γ 1/2 1 -exp - γr 2 .
(3.36)

This information of new isotropic internal variable along with eq. (2.45b) are used to establish the relation between R and R, i.e.

R = R∞ R R∞ γ 1/2 2 2 - R R∞ γ 1/2 2 .
(3.37)

The new isotropic variables can be used to interpret eq. (3.22) as

β R = Q 0 0 R∞ α r , (3.38) such that Z = β R , Π = Q 0 0 R∞ , X = α r .
(3.39)

A state law following such a "normal formulation" can easily be incorporated in the linear global step of the LATIN method.

LATIN method in a heuristic nutshell

To summarise, LATIN method can be viewed as a linearisation scheme, with inherent fully discrete formulation. Instead of classical solution schemes that use time step by time step strategy, LATIN method gives an approximation of the quantities of interest for all time. Separate linearisation techniques are not necessary because of the separation of the linear and non-linear part which directly provides a linear formulation of the global equations. The local stage (or the solution of evolution equations) is extremely cheap as no time integration is needed and the primal quantities of interest are obtained in rate form by solving algebraic equations and the integration is done only in the global stage. The cost of solving the global equilibrium is reduced due to the separated representation with respect to space and time. The POD phase of the global stage, i.e. where the spatial basis is reused, is extremely cheap as there are no iterations involved and only the temporal weak form is solved. The choice of enriching the basis is governed by inherent criteria, and this being one of the most costly part of the algorithm, is performed if it is absolutely essential. The calculation of the direction of descent is also expensive and thereby is not calculated at every LATIN iteration but only if necessary. The choice of this operator may depend on the degree of non-linearity. For instance, if the problem involves low (visco)plasticity then even the inverse of Hooke's tensor may provide appreciable convergence however for high degree of plasticity the search direction operator should be close to the exact tangent operator.

The comparison of Galerkin-based PGD formulation and formulation based on residual minimisation including the respective benefits and drawbacks has been well explained in [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]. Although his model problem was based on advection diffusion reaction equation, Nouy has drawn excellent insight on the two formulations in a generic sense. Formulation based on minimal residual criterion can be considered to be more robust than Galerkin-based PGD formulation as there is a strict monotonic decrease of the residual and hence a guaranteed convergence. However, the convergence rate is depended on the construction of a suitable norm of the residual, which in many cases not being straightforward, renders the convergence rate much lower than Galerkin-based PGD formulation. The computational effort is dramatically increased in the case of minimal residual formulation because the separation order is higher compared to Galerkin-based formulation. The temporal problem in case of minimal residual formulation also is a second order ordinary differential equation (ODE) with initial and final boundary conditions compared to the first order ODE with initial condition obtained for Galerkin formulation. Although it is recommended in [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF] to avoid residual minimisation-based PGD for cases where Galerkinbased PGD works, the guaranteed convergence of the minimal residual-based PGD cannot be undermined. Henceforth, in this research a hybrid formulation will be used that benefits from both types of PGD formulation.

Although LATIN method has been introduced in this chapter in the context of small displacement (infinitesimal strain) problems, LATIN-PGD methods have been used to solve large displacement (finite strain) problems (see [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF], in the case of buckling (see [START_REF] Boucard | A nonincremental approach for large displacement problems[END_REF] and also for elastomers (see [START_REF] Aubard | Modeling and simulation of damage in elastomer structures at high strains[END_REF]. LATIN-PGD has also been used in the field of parametric viscoplastic problems (see [START_REF] Relun | A model reduction technique based on the PGD for elasticviscoplastic computational analysis[END_REF] and also for multi-physics problems (see Néron and[START_REF] Néron | A computational strategy for poroelastic problems with a time interface between coupled physics[END_REF][START_REF] Néron | Proper generalized decomposition for multiscale and multiphysics problems[END_REF]. Computation of large frequency band in transient dynamics (see [START_REF] Barbarulo | Proper Generalized Decomposition applied to linear acoustic: A new tool for broad band calculation[END_REF] and solution inverse problems (see [START_REF] Allix | A new multi-solution approach suitable for structural identification problems[END_REF] has also been performed using LATIN-PGD algorithm. It has also been used for structural domain decomposition (see Ladevèze andNouy, 2003, Ladevèze et al., 2007) and to solve frictional contact problems (see [START_REF] Giacoma | Toward an optimal a priori reduced basis strategy for frictional contact problems with latin solver[END_REF]. A reference point method has also been introduced in [START_REF] Capaldo | The Reference Point Method, a hyperreduction technique: Application to PGD-based nonlinear model reduction[END_REF], where the idea of the hyperreduction method has been used in a LATIN-PGD framework. As far as damage is concerned, LATIN has been used in [START_REF] Allix | A damage prediction method for composite structures[END_REF] for delamination in composite structures, in Bellenger andBussy (1998, 2000) for metal forming processes. The main limitation of LATIN method, according to the author's point of view is its intrusive nature, especially when used in commercial finite element software or established finite element code. Most of the finite element software and academic finite element codes used in different academic and research institutions use classical time incremental schemes based on Newton-Raphson algorithms. Thereby, incorporation of LATIN method is not straightforward as the complete structure of the codes needs replacement and the complete notion of incremental finite element analysis needs to be re-structured. This according to the opinion of the author, is the reason behind the limited usage of the LATIN method in spite of its tremendous advantages. However, recently there have been works in LMT-Cachan to introduce LATIN in commercial finite element software.

Newton-Raphson technique in the light of LATIN method

With the background of LATIN method, the classical Newton-Raphson technique can be interpreted in the same perspective. In other words, Newton-Raphson methods can be viewed exactly as LATIN algorithms but not on whole space-time domain but at each time step (see [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF]Ladevèze, 1993, Ladevèze, 1999, for details). If a solution set -global stage: Seeking s k+1 n+1 from ŝ k n+1 such that the "search direction" equation In a certain sense both these algorithms can be viewed in the same light, but the two nested loops are swapped. In Newton-Raphson algorithm, the external loop sweeps the temporal domain from the first time point to the last, and at each time point the inner loop consists of iterations between the global admissibilities and local constitutive behaviour. On the contrary, in LATIN algorithm, the external loop iterates between the global admissibilities and local constitutive behaviour for all the time points, whereas the inner loops sweep the temporal domain.

s = [ε ( x) , σ ( x)] , ∀ x ∈ Ω (3.
σ k+1 n+1 -σ k n+1 = C T σ k+1 n+1 -σ k n+1 is

Concluding remark

To predict lifetime of structures under fatigue loading, model reduction strategies specific to LATIN method are important components in the solution philosophy, as the separated representation not only reduces the computation but also let the temporal part be treated separately. The strategy benefits from adaptive enrichment of the ROB and henceforth the ROM, dictated by the desired non-linear problem without any pre-computation.

In this research, the lifetime of a structure will be quantified by the evolution of damage variable, hence it is necessary to combine continuum damage mechanics and large time increment method to benefit from both of their advantages. Using damage as an extra internal variable in the context of LATIN algorithm is not straightforward as the elastic state law becomes non-linear. The detailed discussion regarding the inclusion of damage variable in LATIN framework will be carried out in the next chapter.

Chapter 4

LATIN-PGD technique for cyclic damage simulation

This chapter essentially deals with the strategy concerning the embodiment of damage as an internal variable in the LATIN framework. Unlike the usual elasto-(visco)plastic problems tackled by LATIN method as described in section 3.3, the elastic state law due to the presence of the damage variable is nonlinear and hence cannot be solved in the linear stage. The proposed idea is to solve the non-linear elastic state law in the local stage. This chapter is based on [START_REF] Bhattacharyya | A LATIN-based model reduction approach for the simulation of cycling damage[END_REF] with certain elaboration wherever necessary.

The chapter begins with a summary of all the equations that need to be satisfied, thereby moves on to the solution strategies in the local and global steps of the LATIN method. An elaborative formulation of the global step that includes a hybrid technique to solve the PGD problem and the novel idea to solve the non-linear elastic state law are also detailed. The chapter ends with a few illustrative numerical examples to display the functionality of the method.

The proposed problem

The proposed problem to be solved is a quasi-static evolution of a structure that is described in chapter 2. The solution set s defined in chapter 3 must be expanded to encompass all the field variables that define the state of the structure, i.e. s = εp , ε e , Ẋ, Ḋ, σ, Z, Y . The material chosen for this analysis is based on Marquis-Chaboche elasto-viscoplastic model with kinematic, isotropic hardening and coupled with damage. The set of all equations, with ̺ = σ D 1 -D β, that must be satisfied by the structure is summarised in table 4.1.

Table 4.1: The set of constitutive relations and admissibility conditions that need to be satisfied viscoplastic yield function

f p = 3 2 [̺ : ̺] + a 2 β : Q -1 β -R∞ R R∞ γ 1/2 2 2 - R R∞ γ 1/2 2 -σy damage function f d = Y -Y th evolution equations εp = k f p n +   3 2 ̺ 3 2 [̺ : ̺]   1 1 -D α = -k f p n +   - 3 2 ̺ 3 2 [̺ : ̺] + aQ -1 β   ṙ = -k f p n + R R∞ γ 2 -γ 1/2 Ḋ = k d f d n d + state laws ε e ij = 1 + ν E σ + ij 1 -D + σ - ij 1 -hD - ν E σ kk 1 -D + -σ kk 1 -hD δ ij β ij = Qα ij R = R∞ r Y = 1 + ν 2E σ + ij σ + ij (1 -D) 2 + h σ - ij σ - ij (1 -hD) 2 - ν 2E σ kk 2 (1 -D) 2 + h -σ kk 2 (1 -hD) 2 static admissibility Ω×[0,T ] σ : δ εdV dt + Ω×[0,T ] f d • δ ˙ u dV dt + ∂Ω 2 ×[0,T ] F d • δ ˙ u dSdt = 0 kinematic admissibility - Ω×[0,T ] δσ : ε dV dt + ∂Ω 1 ×[0,T ] δσ • n • ˙ U d dSdt = 0

Initialisation

The algorithm starts with the usual elastic initialisation. A solution set s 0 ∈ A d is calculated, considering the loading to be elastic. The elastic assumption ensures that all the internal variables X, Z, D and ε p are set equal to zero. The displacement and/or force boundary conditions are taken into account in this phase, and each subsequent LATIN iteration will involve correction to the assumed linear elastic solution.

Local stage

The local stage essentially involves solution of all the local non-linear equations, i.e. the evolution equations and the non-linear elastic state law (see table 4.1). The objective thereby, is with the knowledge s i ∈ A d , to find ŝi+1/2 ∈ Γ. The search direction equations that were introduced in eq. (3.24), are expanded to incorporate the ancillary variables of the solution set s, i.e.

  

εp i+1/2 -εp i -Xi+1/2 -Ẋi εe i+1/2 -ε e i    + B + σi+1/2 -σ i Ẑi+1/2 -Z i σi+1/2 -σ i = 0, (4.1a) Di+1/2 -Ḋi + b + Ŷi+1/2 -Y i = 0. (4.1b)
The new directions of ascent B + and b + are chosen based on eq. (3.26), i.e. to consider

B + -1 = 0, and b + -1 = 0. (4.2)
This choice of the search direction operators provides

σi+1/2 = σ i , Ẑi+1/2 = Z i , Ŷi+1/2 = Y i . (4.3)
The thermodynamic forces being obtained, the next step is to solve for the internal variables. The LATIN iteration index i+1/2 used for the local step is dropped henceforth for further formulations in this section to circumvent confusion and maintain clarity. The solution procedure begins with solving

D = k d f d n d + , ( 4.4)
where

f d = Ŷ -Y th , (4.5) to obtain D. If, the time domain [0, T ] is discretised into ñ+1 time instances i.e. [t 0 = 0, t 1 , t 2 , • • • , tn, • • • , t ñ-1 , t ñ = T ],
Dn at time instance tn is obtained using a backward Euler rule, i.e.

Dn = Dn-1 + k d Ŷ n -Y th n d + ∆tn, ∀n = {1, 2, • • • , ñ -1} , ( 4.6)
with ∆tn = tn -t n-1 . The quantity D0 is taken to be zero for virgin material, and contains the initial damage value for pre-damaged cases. The damage variable is calculated for all Gauss points. Thereafter, this information is used to calculate the elastic strain tensor, i.e.

εe ij = 1 + ν E σ + ij 1 - D + σ - ij 1 -h D - ν E σkk 1 - D + -σ kk 1 -h D δ ij . (4.7)
The algebraic formulation is solved to obtain εe for every Gauss point at all time instances. The evolution equations

εp = k f p n +   3 2 ̺ 3 2 [ ̺ : ̺]   1 1 - D (4.8a) α = -k f p n +   - 3 2 ̺ 3 2 [ ̺ : ̺] + aQ -1 β  (4.8b) r = -k f p n + R R∞ γ 2 -γ 1/2 (4.8c) with ̺ = σD 1 - D -β (4.9)
and the yield function

f p = 3 2 [ ̺ : ̺] + a 2 β : Q -1 β -R∞ R R∞ γ 1/2 2 2 - R R∞ γ 1/2 2 -σy (4.10)
are also solved algebraically at each Gauss point for every time instance. It has to be noted that integration of the evolution equations (eq. (4.8)) is not needed in the local stage as the plastic strain and the hardening variables are obtained in rate form. This culminates into the complete solution set ŝi+1/2 , the knowledge of which is quintessential in obtaining the global solution set s i+1 .

Search direction for global stage

The objective of the global stage is to construct a solution set s i+1 ∈ A d . The first set of equations that needs to be considered is the search direction equations. Equation (3.25) is thereby modified as

   εp i+1 -εp i+1/2 -Ẋi+1 -Xi+1/2 ε e i+1 -εe i+1/2    -B - σ i+1 -σi+1/2 Z i+1 -Ẑi+1/2 σ i+1 -σi+1/2 = 0, (4.11a) Ḋi+1 -Di+1/2 -b -Y i -Ŷi+1/2 = 0, (4.11b)
where B -= H -0 0 C -1 . The operator H -belongs to the tangent space associated with the solution set ŝi+1/2 in the manifold Γ. For mixed hardening, it is of the form

H -= Hσ H σβ H σR H βσ H β H βR H Rσ H Rβ H R .
(4.12)

To simplify the operator, only the diagonal terms are taken into account and the off-diagonal terms are considered to be zero, which results in

H -= Hσ 0 0 0 H β 0 0 0 H R = Hσ 0 0 H Z , (4.13)
with,

H Z = H β 0 0 H R . (4.14)
The diagonal terms of H -can be obtained from the pseudo-dissipation potential as

Hσ = ∂ ∂σ ∂Φ * ∂σ s=ŝ = kn f p n-1 1 1 -D 2 3/2 ̺ ς ⊗ 3/2 ̺ ς + k f p n 1 1 -D 2 3/2 ςI -3/2 ̺ ⊗ ̺ ς ς2 , (4.15a) H β = ∂ ∂β ∂Φ * ∂β s=ŝ = kn f p n-1 - 3/2 ̺ ς + a Q I β ⊗ - 3/2 ̺ ς + a Q I β + k f p n 3/2 ςI -3/2 ̺ ⊗ ̺ ς ς2 a Q I, (4.15b) H R = ∂ ∂ R ∂Φ * ∂ R s=ŝ = kn f p n-1 γ 1 - Rγ 1/2 2R∞ 2 + k f p n γ 2R∞ , (4.15c)
with ς = 3 2 ̺ : ̺ and I is identity matrix. It can be noted that the index i + 1/2 is dropped in eq. ( 4.15)

for convenience. The definitions of the search direction operators were established in [START_REF] Relun | A model reduction technique based on the PGD for elasticviscoplastic computational analysis[END_REF] for perfect viscoplasticity problems, and extended for hardening and damage in [START_REF] Bhattacharyya | A LATIN-based model reduction approach for the simulation of cycling damage[END_REF].

The operator H -hence obtained is symmetric, positive and semi-definite. The inverse of this operator is not always defined, e.g. at the spatial and temporal points where the solution is elastic. Therefore it is regularised as

H -= H -+ ̟N -1 , (4.16)
with ̟ being the regularisation coefficient and

N = C 0 0 0 Q 0 0 0 R∞ . (4.17)
With these search direction operators being obtained, the quantities of interest can be calculated.

Internal variables at the global stage

The internal variables being local in space, are calculated without any PGD approximation, but by simply solving first order ODEs in time, locally at each Gauss point (GP). Some early versions of LATIN employed PGD-based technique to solve the internal variables as well (see Ladevèze, 1999, e.g.), however in this thesis PGD on internal variables is not used. The search direction equation for hardening variables eq. (4.11a) combined with the state equation can be written as

-Ẋi+1 -Xi+1/2 = H Z Z i+1 -Ẑi+1/2 = H Z ΠX i+1 -Ẑi+1/2 , =⇒ Ẋi+1 + H Z ΠX i+1 = Xi+1/2 + H Z Ẑi+1/2 . (4.18)
This equation is solved using the same scheme as in eq. (4.6), i.e.

X n i+1 = 1 ∆tn I + H n Z Π -1 Xn i+1/2 + H n Z Ẑn i+1/2 + 1 ∆tn X n-1 i+1 , ∀n = {1, 2, • • • , ñ} . (4.19)
The associated variables can then be simply calculated as

Z n i+1 = ΠX n i+1 , ∀n = {0, 1, • • • , ñ} . (4.20)
As far as the damage variables are concerned, the search direction operator b -being considered to be zero,

Ḋi+1 = Di+1/2 =⇒ D i+1 = Di+1/2 . (4.21)
Finally, knowing σ i+1 and D i+1 , the energy release rate Y i+1 is estimated using

Y n = 1 + ν 2E σ n + ij σ n + ij (1 -D n ) 2 + h σ n - ij σ n - ij (1 -hD n ) 2 - ν 2E σ n kk 2 (1 -D n ) 2 + h -σ n kk 2 (1 -hD n ) 2 , ∀n = {0, 1, • • • , ñ -1} . (4.22)
It can be noticed that the LATIN iteration index i + 1 is dropped in eq. ( 4.22) to avoid confusion with the matrix indices. Although trivial, care should be taken to construct the positive and negative parts of the stress tensor at each time point tn, as explained in section 2.2.5.

PGD formulation of the global stage

The idea here is to solve the global admissibilities using PGD-based model reduction methods, where the quantities of interest are represented in separated variable forms. The boundary conditions being taken into consideration in the elastic initialisation, the solution set which is searched here in terms of corrections has to be kinematically admissible to zero. The static admissibility condition can be expressed ∀δ u which is kinematically admissible to zero as

[0,T ]×Ω ∆σ i+1 : ε (δ u) dV dt = 0, (4.23)
where ∆σ i+1 = σ i+1σ i . In the presence of a linear elastic state law, (4.24) where ∆ε i+1 = ε i+1ε i and ∆ε p i+1 = ε p i+1ε p i . This allows the reformulation of eq. (4.23) into

∆σ i+1 = C ∆ε i+1 -∆ε p i+1 ,
[0,T ]×Ω C∆ε i+1 : ε (δ u) dV dt = [0,T ]×Ω C∆ε p i+1 : ε (δ u) dV dt. (4.25)
This weak form leads to a straightforward PGD representation of the stress and plastic strain corrections.

The quantities of interest are described by a single temporal basis and by separate spatial bases related to each other through constant linear FE operators (see [START_REF] Relun | Multiscale elastic-viscoplastic computational analysis[END_REF][START_REF] Relun | A model reduction technique based on the PGD for elasticviscoplastic computational analysis[END_REF], for details). In the presence of damage however the elastic state law is non-linear and the above formulations are not possible. To circumvent this difficulty, the following formulation is proposed. The stress correction term ∆σ i+1 is severed into ∆σ ′ i+1 that relies on plastic deformation and ∆σ i+1 that is procured from the non-linear elastic state law,

∆σ i+1 = ∆σ ′ i+1 + ∆σ i+1 . (4.26)
The total strain correction obeying the strain partition relation given by eq. (2.18) is also sundered into ∆ε ′ i+1 depending on plasticity and ∆ε i+1 relying on damage,

∆ε i+1 = ∆ε ′ i+1 + ∆ε i+1 . (4.27)
Now, the corrections to the quantities of interest are used to reformulate the search direction equation pertinent to the elastic strain tensor as

∆σ i+1 = C∆ε e i+1 -∆R i+1 , (4.28)
where ∆ε e i+1 = ε e i+1ε e i , and ∆R i+1 given by

∆R i+1 = σ i -σi+1/2 -C ε e i -εe i+1/2 , (4.29)
represents a residual stress term at iteration i + 1. Using eqs. (4.26) to (4.29), one can establish

∆σ ′ i+1 + ∆σ i+1 = C ∆ε ′ i+1 -∆ε p i+1 + C ∆ε i+1 -∆ε R i+1 , (4.30)
where ∆ε R i+1 given by ∆ε R i+1 = C -1 ∆R i+1 , (4.31) can be delineated as a residual strain obtained from non-linear state law.

This additive split of the stress and the total strain leads to separate satisfaction of the static admissibility condition. Equation (4.23) can thereby be split into two separated weak forms, the first one depending on the plastic deformations can be expressed ∀δ u which is kinematically admissible to zero as

[0,T ]×Ω ∆σ ′ i+1 : ε (δ u) dV dt = 0, (4.32) with ∆σ ′ i+1 = C ∆ε ′ i+1 -∆ε p i+1 . (4.33)
The second one that is related to damage can be expressed similarly ∀δ u which is kinematically admissible to zero as As a result of the separation of the stress tensor, eq. (4.37) is rewritten using the quantities depending on plastic deformation, as

∆ εp i+1 -Hσ∆σ ′ i+1 + ∆ε i+1 ≈ 0. (4.38)
The intention thereby is to solve for ∆ εp i+1 , ∆σ i+1 , ∆ε e i+1 with additional satisfaction of the admissibility conditions. It has to be noted, however, that the search direction equation given by eq. (4.38) is an exact equality only in the absence of damage, where ∆σ i+1 = ∆σ ′ i+1 . To obtain physically and mathematically viable solutions, it is necessary that the degree of inequality of eq. (4.38) is not too high, which basically hints that this approach is not feasible for high value of damage. However, for metals, the critical damage level being around 0.2 (see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF], the value of damage in general is low enough to render the equality approximation of eq. (4.38) reasonable.

A different strategy of addressing the non-linear elastic state law is presented in Appendix D by splitting the effective Hooke's tensor using singular value decomposition, which gives a complete satisfaction of the static admissibility for each of the snapshots iteratively at a given global stage. However the major limitation of the strategy as per Appendix D, is the fact that the Hooke's tensor can be decomposed in time and space only if the loading is proportional.

Separable representation of the quantities of interest

The quantities of interest are currently ∆σ ′ i+1 , ∆σ i+1 , ∆ εp i+1 , and ∆ε e i+1 . The plastic strain can be written as ∆ε p i+1 = λ p (t) εp ( x) ⇒ ∆ εp i+1 = λ p (t) εp ( x) .

(4.39)

A similar separation on ∆ε ′ i+1 and the corresponding displacement field lead to

∆ u i+1 = λ u (t) ¯ u ( x) ⇒ ∆ε ′ i+1 = λ u (t) ε ( x) . (4.40)
It needs to be mentioned that although only one space-time pair is used here to approximate each correction, it can easily be extended for more separated terms (see Néron and Ladevèze, 2010, for details).

For the following development, the indices are dropped for simplicity and used if necessary. Introducing the PGD approximation, the variation in the displacement field in eq. ( 4.40) becomes δ (∆ u) = δλ u ¯ u + λ u δ ¯ u with δ ¯ u being kinematically admissible to zero and δλ u does not have any condition. The space-time problem defined by eq. (4.32) is thereby separated into a spatial problem and a temporal problem. The temporal problem consists in calculating λ u such that ∀δλ u ,

[0,T ] λ u δλ u dt Ω Cε : ε ¯ u dV = [0,T ] λ p δλ u dt Ω Cε p : ε ¯ u dV. (4.41)
The proportionality of λ u and λ p is a direct outcome of eq. (4.41), which lets, according to Relun et al. (2011), the most trivial assumption λ u = λ p = λ . This preference of the temporal functions leads to the spatial problem be defined as

[0,T ] λ λ dt Ω Cε : ε δ ¯ u dV = [0,T ] λ λ dt Ω Cε p : ε δ ¯ u dV, (4.42)
∀δ ¯ u which is kinematically admissible to zero. This can be further reduced to

Ω Cε : ε δ ¯ u dV = Ω Cε p : ε δ ¯ u dV, (4.43)
allowing the construction of an operator E in classical FE space such that ε = Eε p . (4.44)

Finally, ∆σ ′ i+1 and ∆ εp i+1 , in a separated variable form can be represented as

∆ εp i+1 = λ (t) εp ( x) , (4.45a) ∆σ ′ i+1 = λ (t) Cε p ( x) . (4.45b)
The construction of the FE operators E and C is explained in Appendix A.

The quantities ∆σ i+1 and ∆ε i+1 depending on damage are solved from the weak form eq. (4.34) together with eq. (4.35) providing

[0,T ]×Ω C ∆ε i+1 -∆ε R i+1 : ε (δ u) dV dt = 0, (4.46)
where ∆ε R i+1 is a quantity already known from the local stage and δ u is kinematically admissible to zero. ∆ε i+1 can thereby be calculated using the same FE operator E as

∆ε i+1 = E∆ε R , (4.47)
and ∆σ i+1 is calculated as

∆σ i+1 = C∆ε R . (4.48)
Lastly the corrections to the stress and elastic strain are given as

∆σ i+1 = ∆σ ′ i+1 + ∆σ i+1 , (4.49a
)

∆ε e i+1 = ∆ε ′ i+1 + ∆ε i+1 -∆ε p i+1 . (4.49b)
4.6.2 Hybrid method to construct the PGD reduced-order basis

As mentioned before, the calculation strategy of the separated forms involves both POD and PGD phases. The POD phase involving the update of the temporal functions by reusing the spatial basis, is relatively less expensive and performed at every global stage of the LATIN method. If the ROM hence obtained is not adequate enough, the ROB is enriched through the PGD phase by the addition of a new space-time mode. A hybrid strategy involving a Galerkin-based formulation to calculate the spatial problem and a formulation based on minimum residual to approximate the temporal problem is used.

Update of the reduced order basis

The idea of update stage at LATIN iteration i + 1 is to reuse the m space functions that have already been generated from LATIN iteration i and update the corresponding time functions. The quantities of interest depending on plastic deformation are written as

∆ εp i+1 = m j=1 ∆ λj (t) εp j ( x) , ∆σ ′ i+1 = m j=1 ∆λ j (t) Cε p j ( x) , (4.50)
where {∆λ j } m j=1 represent the updates of the temporal functions. With the information of eq. (4.50), the constitutive relation error associated with eq. (4.38) is formulated as

e 2 CE = ∆ εp i+1 -Hσ∆σ ′ i+1 + ∆ε i+1 2 H -1 σ , = m j=1 ∆ λj εp j -Hσ m j=1 ∆λ j Cε p j + ∆ε i+1 H -1 σ . (4.51)
Here {∆λ j } m j=1 are the only unknowns and are approximated such that the constitutive relation is minimised, i.e. This results in multi-variable coupled differential equations. Zero-order discontinuous Galerkin method is used to solve this type of equations, as a better approximation compared to shooting method was obtained in [START_REF] Passieux | Time-space radial approximation and multiscale LATIN method[END_REF], [START_REF] Relun | Multiparametric strategy for robust design in fatigue[END_REF]. The method of solving the minimisation problem is detailed in Appendix B. Thereafter, knowing {∆λ j } m j=1 , the total stress and elastic strain are approximated using eq. (4.49). This gives a ROM based on the same ROB. If this approximation is adequate, the particular global step is terminated and the algorithm moves to the next local step. The criterion to enrich the ROB is estimated by a saturation parameter ζ of the error indicator ξ (see section 4.7) which is similar to the criterion proposed in [START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF],

ζ = ξ i -ξ i+1 ξ i + ξ i+1 . (4.53)
If ζ is less than a tolerance level ζ tol at a particular LATIN iteration, a space-time PGD pair needs to be added to enrich the reduced-order basis. On the contrary, if ζ is high enough, the ROM is considered to be effective.

Enrichment of the reduced-order basis

The enrichment of the ROB is performed using a hybrid strategy. The temporal functions are calculated using a minimisation of residual (as mentioned in section 4.6.2.1), and the spatial functions are computed using a Galerkin technique.

The search direction equation during the stage of enrichment becomes such that ∆ εp i+1 is obtained from eq. (4.54) and ∆σ ′ i+1 from eq. (4.56b). For convenience, all the indices indicating the number of modes and number of LATIN iteration are dropped for convenience. The kinematic admissibility condition is written ∀δσ which is statically admissible to zero, as A pseudo-strain ε is defined thereafter such that

∆ εp i+1 -Hσ∆σ ′ i+1 + ∆ε i+1 = 0, (4.54) with ∆ε i+1 = Hσ σi+1/2 -σ up i+1 -εp i+1/2 - εp , up i+1 , ( 4
ε = ∆ ε′ λ = W -1 σ -δ, (4.63)
with W -1 = Hσλ 2 + λ λ C -1 and δ = ∆ε λ . Using the introduced definition of the pseudo- strain, eq. ( 4.62) can be rewritten as ∀δ ¯ u that is kinematically admissible to zero and σ that is statically admissible to zero. This static admissibility condition is re-written with respect to the pseudo-strain, i.e. using eq. (4.59) and by introducing a pseudo-displacement field ¯ u which is kinematically admissible to zero, such that ε

¯ u = ε, Ω Wε ¯ u : ε δ ¯ u dV = - Ω W δ : ε δ ¯ u dV = 0, (4.66)
∀δ ¯ u that is kinematically admissible to zero. This spatial problem is solved using classical FE scheme to obtain the pseudo-displacement ¯ u. The associated pseudo-strain ε can then be easily calculated using the corresponding derivatives of the shape functions. Knowing ε, eqs. (4.59) and (4.63) can be used to calculate the space function εp m+1 as

εp m+1 = 1 λ λ ε -λ λ C -1 W ε + δ . (4.67)
A fixed-point iteration is used to approximate εp m+1 ( x) and λ m+1 (t) which is summarised in algo- rithm 1, with the number of sub-iterations s it taken to be five.

The new space function εp m+1 is orthonormalised with respect to the previously existing spatial basis εp j m j=1 using Gram-Schmidt algorithm (see [START_REF] Ruhe | Numerical aspects of Gram-Schmidt orthogonalization of vectors[END_REF]. In the numerical process all the former time functions {λ j } m j=1 are updated along with the new time function λ m+1 . If the modified time function λ m+1 has an insignificant norm the corresponding space-time pair may be rejected (see Appendix C).

Relaxation of the solution field and convergence criterion

To secure convergence of the LATIN-PGD algorithm, the global step needs to be relaxed. Represented by si+1 , the solution set obtained at the end of the global step i + 1, which was previously being denoted by s i+1 , the new solution set s i+1 can be formulated as (4.70)

s i+1 = ϕs i+1 + (1 -ϕ) s i , ( 4
It has to be noted that the relaxation parameter ϕ, the regularisation coefficient ̟, and the error saturation tolerance ζ tol mostly affect the rate of convergence and not the convergence itself.

Numerical examples

The innovative algorithm has been tested on academic examples in one-dimensional and two-dimensional cases. The relaxation parameter ϕ and regularisation coefficient ̟ are valued to be 0.8 and 0.15 respectively (from [START_REF] Relun | Multiscale elastic-viscoplastic computational analysis[END_REF].

Bar under traction

The one-dimensional test problem considered is a bar with length L = 1000 mm and area of cross section A = 100 mm 2 as depicted in fig. 4.1. The bar is constrained at x = 0 and a sinusoidal prescribed displacement loading of amplitude 1.2 × 10 -3 L is applied at x = L with a time period ∆T = 10 s for 20 cycles. The structure is discretised in space using 90 linear bar elements such that each part has the same number of elements. A time step size of 0.1 s is chosen for temporal discretisation, which results in 100 time elements per cycle and in total 2000 time elements. In the computational model, the criterion ζ tol to enrich the reduced-order basis is taken to be 0.01. The algorithm is stopped if the saturation parameter ζ is lower than 10 -4 . For one-dimensional problems, the search direction operators defined in eq. (4.15) become

Hσ = kn f p n-1 1 (1 -D) 2 , (4.71a) H β = kn f p n-1 -sign σ 1 -D -β + a C β 2 + k f p n a C , (4.71b) H R = kn f p n-1 γ 1 - Rγ 1/2 2R∞ 2 + k f p n γ 2R∞ . (4.71c)
The convergence of the algorithm is measured by the LATIN error indicator ξ depicted in fig. 4.2. Towards the end of the algorithm the LATIN indicator ξ shows a stagnating phenomenon, and neither updating the reduced basis or its enrichment will improve the solution significantly. The orthonormalisation procedure presented in Appendix C results in generation of a total of three PGD pairs due to rejection of insignificant modes. The quantity of interest εp obtained as the sum of products of space-time modes after convergence is shown in fig. 4.3. The internal variables in rate form and the corresponding associated variables for the three sections of the bar are depicted in fig. 4.4. The isotropic hardening variables R and r are referred here as growth stress and growth strain respectively, as they represent the growth of the yield surface. Due to the prescribed boundary conditions, the stress is constant for a given time all over the bar, with decreasing amplitude with time, which basically represents material weakening due to damage. For mat. 1, the plastic strain rate and kinematic hardening variables in rate form show increasing amplitude with respect to time. The amplitudes however decrease for mat. 2 and mat. 3. For mat. 1 there is a monotonic increase in the amplitude of ṙ, and for mat. 2 and mat. 3, the amplitude of ṙ decreases. The associated variable R increases for all the three materials.

The quantities of interest that directly indicate the influence of unilateral condition of micro-defect closure effects are depicted in fig. 4.5. The damage variable D remains practically constant during the compressive part of the loading and increases during the tensile part of the loading. The energy release rate Y also illustrates the difference in behaviour during compression and tension, i.e. very high in tension and extremely low during compression. The elastic strain ε e also shows a higher value during tension than in compression. The values of D at the end of loading for mat. 1, mat. 2 and mat. 3 are 0.22, 0.18 and 0.15 respectively.

General practical engineering problems invariably portray high localisation for both plastic deformation and damage. The presence of stress raisers (e.g. notches, holes, defects and such others) in most engineering components, concentrates the effect of plasticity and damage to a very limited region, while the rest of the structure remains undamaged and elastic. Classical two-dimensional mechanical problems are considered in section 4.8.2 and section 4.8.3 to depict this phenomenon. An "L" shaped structure subjected to a concentrated load represented in fig. 4.6 is considered as a two-dimensional test problem. The structural geometry is defined by length L = 120 mm, and width W = 20 mm. The structure is filleted at the inside corner with fillet radius r = 5 mm to avoid singularity. The thickness of the structure is taken to be 1 mm. A sinusoidal displacement U d (t) of amplitude 1.5 mm is prescribed with a time period ∆T = 10 s, for 5 cycles. The material properties used are the same as given in table 4.1, with σy = 85 MPa. A total of 277 linear two-dimensional quadrilateral isoparametric plane stress elements with 4 Gauss points per element are used to discretise the structure, which generates 337 nodes. A time step size of 0.2 s is chosen for the temporal discretisation. The criterion to enrich the ROB ζ tol is taken to be 10 -2 . Compared to the one-dimensional problem of section 4.8.1, the calculation of the search direction operators for high dimensional cases is considerably expensive. These operators are thereby not calculated at every iteration, but are only updated if there is a saturation of the LATIN error indicator. It is quantitatively measured by the criterion ζ with the tolerance being 10 -2 . The decrease in the LATIN indicator ξ is shown in fig. 4.7. After the first iteration, the indicator reads 5 × 10 -2 . After 21 iterations, the LATIN indicator becomes 9 × 10 -3 with 8 PGD modes being generated. The saturation of the indicator is at ξ = 8 × 10 -4 , with a total of 40 modes being generated. This saturation of LATIN indicator at a relatively high value, compared to fig. 4.2 is due to the fact that a coarse temporal mesh has been used. In case of a finer temporal mesh this LATIN indicator can be lowered even for two-dimensional problems, as shown in the next numerical example (section 4.8.3).

The distribution of the cumulative plastic strain at the end of loading is depicted in fig. 4.8. It is distributed on both sides of the vertical web, with the inner side being predominant and the maximum at the filleted corner of the structure, whereas most part of the structure shows no plastic deformation.

The residual von Mises stress is distributed in a similar manner as the accumulated plastic strain, as shown in fig. 4.9. It is concentrated at the filleted corner with lower valued distribution on both sides of the web.

The distribution of damage variable D, the main quantity of interest, at the end of the loading (t = T ) is represented in fig. 4.10. The distribution is similar to the one of the accumulated plastic strain and residual von Mises stress, i.e. maximum at the filleted corner. The span of damage in the filleted corner, which is the region of interest, after certain load cycles is also represented in fig. 4.10. The damage

Plate with a hole

The goal of this example is to depict a comparative study for virgin and pre-damaged structures and also to present the influence of lower time step size on the LATIN error indicator ξ for two-dimensional structures. The test problem considered is a classic plate with a hole. The rectangular plate of length L = 40 mm, width W = 60 mm and thickness of 1 mm is considered, with a central circular hole of diameter φ = 20 mm. Uniformly distributed sinusoidal displacements of amplitude 0.012 mm and time period ∆T = 10 s are applied on both ends of the plate, as shown in fig. 4.12, for 10 cycles. The symmetry of the structure enables only a quarter of the plate to be considered for the analysis (see fig. 4.12). The material properties are the same as used in section 4.8.2. The structure is discretised using 129 linear two-dimensional quadrilateral isoparametric plane stress elements with 4 Gauss points per element, generating 154 nodes. A finer time step size of 0.1 s is chosen for the temporal discretisation.

The first numerical test is to solve the problem considering all the GPs to be virgin, which is followed by solving the problem considering the GPs which are indicated by * in fig. 4.13 to be pre-damaged and the virgin GPs are marked in •. The decrease in the LATIN indicator ξ for the virgin structure is shown in fig. 4.14. A total of 34 modes are generated to obtain ξ = 9 × 10 -5 in 75 iterations. Expectedly, the indicator is much lower compared to fig. 4.7, because of finer temporal discretisation (see Ladevèze, 1999, for explanation). For the pre-damaged structure, 38 PGD modes are generated in 83 iterations to obtain ξ = 9 × 10 -5 (see fig. 4.15).

The distribution of damage variable D along with the damage spread in the respective regions of interest, for both the cases, at the end of the loading (t = T ) is represented in fig. 4.16. It is obvious from fig. 4.16, that there is a shift of the maximum damaged point in the pre-damaged case, compared to the virgin case. The weakest Gauss point for the virgin structure is GP 1 (see fig. 4.13) and GP 2 (see fig. 4.13) for the pre-damaged structure. The evolution of damage for both of these GPs is represented in fig. 4.17. .17: Damage evolution at the weakest GP of the plate with hole, i.e. GP 1 for the virgin structure, and GP 2 for the pre-damaged structure

Concluding remarks

This is a novel attempt to include unilateral condition of micro-defect closure effect in a LATIN-PGD framework. The algorithm is able to provide appreciable accuracy and robustness as illustrated in the academic examples. In case of high value of damage, on the other hand, this framework will perhaps be not as accurate, however as the focus is on metals, where the maximum damage is not high, this apparent problem is rendered insignificant. This chapter also depicts a hybrid PGD formulation, which according to the knowledge of the author is one of the first in a LATIN framework. However it still does not address the issue at hand which is to tackle large number of cycles. The usage of PGD-based model reduction approach presented here is not enough to reduce numerical cost when large number of cycles is simulated. Therefore a multi-temporal discretisation is needed to effectively simulate fatigue behaviour with less CPU cost which will be detailed in the next chapter.

Chapter 5

Multi-scale temporal discretisation approach

The strategy to compute damage as presented in chapter 4, when applied for cycle by cycle simulations can be extremely time consuming if the problem involves large number of cycles, e.g. 10 5 cycles will have 10 7 time points if 100 time points per cycle are considered and the total computation is unreasonably expensive. To circumvent this problem, a multi-temporal approach is used. The range of problems encompassed here are those where the applied forces or displacements are periodic with a constant time period ∆T , but with amplitudes that may vary with time (see [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF]Ladevèze, 1993, Ladevèze, 1999). To reduce the cost of the mechanical problem, the temporal dependency of the different quantities on the whole time domain [0, T ], which may include several thousands of cycles, is described in parallel by two time scale discretisations: -a coarse temporal discretisation θ defined on the complete temporal domain [0, T ] representing the beginning of every cycle, -a fine temporal discretisation τ i defined within cycle i.

Any instant of interest t can be discretised by the pair θi , τ i as

t = θi + τ i , with i ∈ [0, N -1] and τ i ∈ [0, ∆T ] , ∀t ∈ [0, T ] ,
(5.1)

where N is the number of cycles included in the temporal domain.

The two scales allow to represent in parallel, the slow evolution along the cycles using the coarse temporal discretisation and the rapid evolution within a specific cycle by the fine temporal discretisation. This type of description allows the definition of "nodal cycles" where the quantities of interest are calculated in the fine temporal scale. Two successive "nodal cycles" bound what is called a "time element", where the evolution of the quantities of interest is considered to be slow allowing a finite element like interpolation on the coarse scale of the temporal quantities. The chapter starts with the description of the finite element like interpolation scheme in time. Thereafter, the calculation of a "nodal cycle" is detailed and is followed by various academic test cases, including numerical verification and comparison with mono-scale LATIN-PGD method.

Finite element like time interpolation scheme

Using the two-time scale description, the constitutive relation and the admissibility condition are solved using LATIN-PGD technique (as described in chapter 4) only on the "nodal cycles". The time element m + 1 is bounded by two "nodal cycles" defined by their initial times Θ m and Θ m+1 respectively. For the sake of clarity a change of variables is introduced for each "time element" containing p cycles through θ k , ∀k ∈ [0, p -1] which represents the coarse discretisation for the given "time element" locally. This representation provides the time t to be defined as

t = θ k + τ k , with k ∈ [0, p -1] and τ k ∈ [0, ∆T ] , ∀t ∈ [Θ m, Θ m+1 + ∆T ] .
(5.2)

Once the successive "nodal cycles" m and m+1 have been computed, a one-dimensional interpolation is used to estimate the quantities of interest on [Θ m, Θ m+1 ] based on temporal shape functions (see [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF]. Linear one-dimensional temporal shape functions ν m and ν m+1 for the time element [Θ m, Θ m+1 ] are defined as

ν m (θ k ) = Θ m+1 -θ k Θ m+1 -Θ m , ν m+1 (θ k ) = θ k -Θ m Θ m+1 -Θ m , ∀t ∈ [Θ m, Θ m+1 + ∆T ] .
(5.3) Then, for any temporal quantity of interest χ, its approximation within the time element [Θ m, Θ m+1 ] is written as a polynomial expansion based on pre-defined temporal shape functions {ν m, ν m+1 } and on the knowledge of the quantity for the two "nodal cycles" denoted as χ (Θ m + τ m) and χ (Θ m+1 + τ m+1 ). Similarly to conventional finite elements, the interpolation reads

χ t θ 0 = Θ m τ m θ 1 = θ 0 + ∆T θ 2 = θ 1 + ∆T θ k τ k θ p-1 = Θ m+1 τ m+1 Θ m+1 + ∆T nodal cycle m + 1 nodal cycle m
χ (t = θ k + τ k ) = ν m (θ k ) χ (Θ m + τ m) + ν m+1 (θ k ) χ (Θ m+1 + τ m+1 ) , ∀t ∈ [Θ m, Θ m+1 + ∆T ] ,
(5.4)

with τ m = τ m+1 = τ k belonging to [0, ∆T ] as the period is constant. The interpolation is more tedious than with conventional finite element as the nodes describe no more some single points but some subdomains. It can be noted that the benefit of the fine time discretisation is transferred to the whole element by the interpolation process.

The finite element like description in time of the quantities of interest lies in a set of time elements which spans the time domain (see fig. 5.1). However, it can be noted that this scheme does not include a Galerkin form of the problem on the whole time domain, but only an a posteriori time interpolation scheme is considered.

The idea is first to calculate initial few cycles classically as given in chapter 4, which can be interpreted as a "training stage" as the ROB obtained will be re-used for the two-scale computation. In details, the first element is built from the "nodal cycle" 0, which is the last classically computed cycle 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. 5.2). The computation is pursued for the different "nodal cycles" in succession until the ultimate "nodal cycle" concluding the time domain of interest is achieved, or until the critical damage has been reached. The computation of each "nodal cycle" involves solution of the global admissibility and the local equations using the LATIN-PGD algorithm only in the corresponding temporal domain. As far the material model is considered, the same model as given in table 4.1 is chosen but now only with kinematic hardening. The constitutive relations are rewritten in table 5.1, the global admissibilities are the same as given in table 4.1.

f p = 3 2 [̺ : ̺] + a 2 β : Q -1 β -σy damage function f d = Y -Y th evolution equations εp = k f p n +   3 2 ̺ 3 2 [̺ : ̺]   1 1 -D α = -k f p n +   - 3 2 ̺ 3 2 [̺ : ̺] + aQ -1 β   Ḋ = k d f d n d + state laws ε e ij = 1 + ν E σ + ij 1 -D + σ - ij 1 -hD - ν E σ kk 1 -D + -σ kk 1 -hD δ ij β ij = Qα ij Y = 1 + ν 2E σ + ij σ + ij (1 -D) 2 + h σ - ij σ - ij (1 -hD) 2 - ν 2E σ kk 2 (1 -D) 2 + h -σ kk 2 (1 -hD) 2

Computation of one "nodal cycle"

Each "nodal cycle" is computed using the LATIN-PGD framework as an independent problem. Therefore, a LATIN initialisation is first handled on this sub-domain. Then, the LATIN iterative scheme is processed, i.e. successive computations of a linear global stage and a non-linear local stage are done.

Initialisation

Contrary to classical LATIN strategy in which the initialisation originates from an elastic computation, using the two-time scale approach, the initialisation for each "nodal cycle" benefits from the knowledge of the previous "nodal cycle". Some quantities of interest, such as the stress, the elastic strain, the internal variables for kinematic hardening, and the time functions representing the plastic strain {λ j } µ j=1 exhibit cyclic behaviour. Therefore, their initialisation for "nodal cycle" m + 1 is based on the duplication of their estimation at "nodal cycle" m. Even if the global tendency of these variables is cyclic, due to the influence of increasing damage, they are not perfectly periodic. Thereby, as it is necessary to avoid any non-physical discontinuity at the transition between successive cycles, the initialisation for any cyclic quantity ϑ (Θ m+1 + τ m+1 ) defined ∀τ m+1 ∈ [0, ∆T ] and for a particular Θ m+1 is based on a transformation of the quantity of interest ϑ (Θ m + τ m) defined ∀τ m ∈ [0, ∆T ] and at Θ m, such that the continuity is preserved. A linear transformation between ϑ (Θ m + τ m) and θ (Θ m + τ m) is considered as

θ (Θ m + τ m) = ϑ (Θ m + τ m) + a ϑ • (Θ m + τ m) + b ϑ , ( 5.5)
where a ϑ and b ϑ are some parameters defined for every quantity of interest such that the periodicity of the transformed quantity of interest is guaranteed, i.e.

θ (Θ m) = θ (Θ m + ∆T ) = ϑ (Θ m + ∆T ) .
(5.6)

The initialisation ϑ 0 (τ m+1 ) for the "nodal cycle" m + 1 is thus based on the transformed quantity computed for the "nodal cycle" m as

ϑ 0 (Θ m+1 + τ m+1 ) = θ (Θ m + τ m) = ϑ (Θ m + τ m) + ϑ (Θ m) -ϑ (Θ m + ∆T ) ∆T • (τ m) + ϑ (Θ m + ∆T ) -ϑ (Θ m) .
(5.7) On another hand, the damage variable which is non-cyclic but a non-decreasing function of time, is initialised as constant over the whole "nodal cycle" m + 1, with the magnitude being the one obtained at Θ m + ∆T . The initialisation of the strain energy release rate Y is calculated from the damage and the stress tensor.

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

Local stage

The main difficulty in the local stage is the evaluation of damage at every Gauss point which needs a time integration scheme. To proceed the integration of the damage evolution equation over "nodal cycle" m + 1, the value of the damage at the beginning of this "nodal cycle", i.e. at Θ m+1 is required.

Therefore, it is required to calculate χ (Θ m+1 ) from χ (Θ m). Assume a temporal evolution of any quantity of interest χ defined by a general first order ODE dχ dt + κχ = υ, (5.8) over a time domain with κ and υ being some time-dependent known parameters. For application for damage, these parameters are material-dependent. From the knowledge of an initial value of the quantity of interest χ (t k ), the solution of the ODE at a further time t l can be determined as (see [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF])

χ (t l ) = χ (t l , t k ) + ℜ (t l , t k ) χ (t k ) ,
(5.9) with χ (t k , t l ) being the solution of the ODE at time t l considering zero initial condition at time t k , and ℜ the "resolvent" operator defined as

ℜ (t l , t k ) = exp t l t k -κ (τ ) dτ .
(5.10)

To compute the solution with zero initial condition and the resolvent operator with the minimum numerical cost, the proposed strategy is to compute these quantities only at the "nodal cycle" m + 1, and then to propagate the information using a discretisation based on the instants θ k representing the beginning of the loading cycles within the time element [Θ m, Θ m+1 ]. Considering a time element comprising p cycles, the discretisation is defined such that

θ k = Θ m + k∆T , where k is a natural number comprised in the interval [0, p -2]. Therefore, it is established that θ 0 = Θ m and θ p-1 = Θ m+1 .
The linear interpolation introduced in section 5.1 allows to obtain χ (θ k+1 , θ k ) and ℜ (θ k+1 , θ k ) from the knowledge of the two "nodal cycles" m and m + 1 delimiting the time element as

χ (θ k+1 , θ k ) = ν m (θ k ) χ (Θ m + ∆T, Θ m) + ν m+1 (θ k ) χ (Θ m+1 + ∆T, Θ m+1 ) , (5.11a) ℜ (θ k+1 , θ k ) = ν m (θ k ) ℜ (Θ m + ∆T, Θ m) + ν m+1 (θ k ) ℜ (Θ m+1 + ∆T, Θ m+1 ) .
(

These interpolated values can be used for approximating the solution of the ODE (eq. (5.8)) as

χ (θ k+1 ) = ν m (θ k ) [ χ (Θ m + ∆T, Θ m) + ℜ (Θ m + ∆T, Θ m) χ (θ k )] + ν m+1 (θ k ) [ χ (Θ m+1 + ∆T, Θ m+1 ) + ℜ (Θ m+1 + ∆T, Θ m+1 ) χ (θ k )] .
(5.12)

In the specific case of the evolution equation for damage, the "resolvent" operator ℜ is equal to one. Considering the solution D (τ m) is known ∀τ m ∈ [0, ∆T ], for the first "nodal cycle", Ď may be easily evaluated as D (Θ 0 + ∆T ) -D (Θ 0 ). Otherwise it is computed only for the second "nodal cycle" Ď (Θ m + ∆T, Θ m) by solving the evolution equation using zero initial condition. Thereby, eq. (5.12) is solved at each θ k+1 till D (θ p-1 ) is calculated. This process can also be represented as

D (θ p-1 ) = D (θ 0 ) + ∆D ini + p-2 k=1 ∆D k ,
(5.13) as illustrated in fig. 5.3 Using the interpolation scheme defined by eq. ( 5.12), the initial value of damage for the "nodal cycle" m + 1 can be written as (5.14) where ∆D ini and ∆D f in are the damage increments for the initial and final cycles of the time element respectively.

D (Θ m+1 ) = D (Θ m) + p-2 k=0 θ m+1 -θ k Θ m+1 -Θ m ∆D ini + θ k -Θ m Θ m+1 -Θ m ∆D f in ,
All the other quantities of interest do not require any time integration. They can be calculated from the non-linear material behaviour in rate form independently for every Gauss point before tackling the admissibility conditions at the global stage.

nodal cycle m nodal cycle m + 1 ∆D f in ∆D 3 ∆D 2 ∆D 1 ∆D ini θ 0 = Θ m θ 1 θ 2 θ 3 θ 4 θ p-1 = Θ m+1 Θ m+1 + ∆T t D Figure 5.3: Calculating D (Θ m+1 ), knowing D (Θ m)

Global stage

The global admissibility conditions are calculated for the whole space domain but only for a subdomain in time, which is the "nodal cycle" m + 1. The idea is exactly similar to section 4.6, the only difference being the temporal domain includes only the "nodal cycle" of interest.

Mainly, the global stage is used only for updating the already built ROB, and enrichment is done if there is an absolute necessity. For simple structures with simple boundary conditions that induce no localisation (see e.g. fig. 5.4), the ROB obtained from the "training stage" is considered sufficient, and any new space-time mode added will be rejected after the orthonormalisation procedure as explained in Appendix C. Thereby the enrichment step is skipped completely to save CPU cost. For structures that may induce localisations (see e.g. figs. 5.14 and 5.17), however, a maximum of one space-time pair is added per "nodal cycle" to enrich the reduced order basis. The enrichment condition is the same as explained in section 4.6.

If {λ m+1 j } m j=1 are the time functions used to describe "nodal cycle" m + 1 at LATIN iteration i with the space functions being {ε p, m+1 j } m j=1 , the objective at LATIN iteration i + 1 is to calculate {∆λ m+1 j } m j=1 which are the corrections of {λ m+1 j } m j=1 . The process is exactly as described in section 4.6.2.1, the only difference being the temporal domain, which is now

[Θ m+1 , Θ m+1 + ∆T ].
If enrichment of the ROB is necessary, a new space-time mode εp, m+1 m+1 , λ m+1 m+1 is added according to the process described in section 4.6.2.2.

As far as the kinematic hardening variables are concerned, eq. (4.18) is rewritten over the "nodal cycle" m + 1 as (5.16)

∆ α m+1 i+1 + H m+1 β Π∆α m+1 i+1 + ∆α, m+1 i+1 = 0 (5.15) with ∆α, m+1 i+1 = H m+1 β Πα m+1 i -β m+1 i+1/2 + α m+1 i -α m+1 i+1/2 ,
and ∆α m+1 i+1 = α m+1 i+1 -α m+1 i , ∆ α m+1 i+1 = α m+1 i+1 -α m+1 i .
(5.17) Equation (5.15) is solved to obtain the corrective terms over the "nodal cycle" m + 1 using zero initial condition as described in section 4.5.

The relaxation of the solution and the convergence indicator for "nodal cycle" m + 1 is the same as given in section 4.7.

Numerical examples

Two-dimensional classical mechanical problems subjected to cyclic loadings are considered to investigate the abilities and the robustness of the proposed numerical scheme. All the two-dimensional structures mentioned herein are discretised in space using linear four-noded isoparametric quadrilateral elements with four Gauss points per element. A one-dimensional problem has also been solved to build virtual ε-N curves.

Verification with mono-scale LATIN method

The first set of analysis is performed for a material behaviour characterised by viscoplasticity and kinematic hardening without any contribution of damage. A rectangular plate described by its length L = 50 mm and width W = 100 mm which is constrained on three sides as shown in fig. 5.4 is considered. The unconstrained side of the plate is loaded with a uniformly distributed sinusoidal displacement of the form U d = U 0 sin 2πt ∆T with the time period ∆T = 1 s for 1000 cycles. The goal is to compare the results obtained using two-time scale algorithm with the results obtained from the mono-scale method described in chapter 4. The material considered is Cr-Mo steel at 25 • C with properties given in table 5.2. The yield stress σy considered for the central part is 188 MPa and for the peripheral material is 189 MPa.

The problem is initially solved for U 0 = 0.073 mm. Four separate temporal discretisations are examined using the two-scale computation, with a minimum of 50 cycles to a maximum of 500 cycles per time element (see table 5.3). The result of this analysis, compared to a mono-scale solution is given in fig. 5.5. Plasticity being the defining phenomenon in the process, the relative error is based on the plastic strain rate, i.e.

er 2 = Ω×[0,T ] εp ref -εp ob : εp ref -εp ob dV dt Ω×[0,T ] εp ref : εp ref dV dt + Ω×[0,T ] εp ob : εp ob dV dt , (5.18)
where εp ref is the mono-scale reference solution and εp ob is the solution obtained from the two-time scale computations. Increasing the size of time elements decreases the accuracy of the approximation, i.e. ob is the time required to obtain the solution using the two-time scale computations. This relative CPU time expectedly decreases with increase in number of cycles per time element. It should also be noted from fig. 5.5 that there is a drastic drop in the calculation time for two-time scale computation, compared to the mono-scale computational time, even if 50 cycles per time element are used.

Similar tests are performed for material that is susceptible to damage along with visco-plasticity and kinematic hardening with the damage threshold considered to be (1.4σy) 2

2E

. The problem is solved for U 0 = 0.070 mm using the four temporal discretisations explained previously. The relative error estimation with respect to a reference mono-scale solution is now based on the damage variable and is given by

er 2 d = Ω×[0,T ] D ref -D ob • D ref -D ob dV dt Ω×[0,T ] D ref • D ref dV dt + Ω×[0,T ] D ob • D ob dV dt , ( 5.19)
where D ref and D ob are obtained from a reference mono-scale and a two-time scale LATIN-PGD algorithms respectively. The relative error, which is plotted in fig. 5.6, strongly increases with respect to the size of time elements. Although it is much higher than the one of the visco-plastic case shown in fig. 5.5, it is still in the acceptable range, i.e. a minimum of 2.6 % and a maximum of 4.2 % is observed. On an other hand, the reduction in CPU time is here even more drastic compared to the case without damage.

Further tests are conducted to analyse the effect of the load level on the accuracy of different temporal discretisations represented in fig. 5.7a. These comparisons are with respect to the results obtained using It is seen that with increase in load level, the coarser temporal discretisation results in higher inaccuracy as the damage behaviour becomes non-linear (see fig. 5.8). However for low damage level, the damage evolution is comparatively linear and coarse temporal mesh provides relatively acceptable accuracy. As far as the computational time is concerned (see fig. 5.7b), there are two aspects that should be explained. First of all, for a given load level, finer temporal mesh will be computationally expensive as the number of "nodal cycles" computed is higher. Secondly, for a given temporal discretisation, increasing the load level will increase the CPU cost as higher load induces higher degree of non-linearity and more LATIN iterations are needed at each "nodal cycle" for convergence. It should be noted that the comparisons depicted in fig. 5.7 are with respect to a two-time scale solution using 50 cycles per time element and not a mono-scale solution.

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 adaptive temporal discretisation scheme based on the damage level should be used .

Influence of the "training stage"

Before using the adaptive scheme, it should also be checked if the number of cycles those are initially calculated for the generation of the ROB has any influence in the solution. Four separate tests are conducted on the same structure given in fig. 5.4 using U 0 = 0.063 mm for 1 × 10 5 cycles with 50, 20, 10, and 5 cycles for the "training stage". A uniform temporal mesh of 100 cycles per time element is considered thereafter. Considering the solution obtained using 50 cycles in the "training stage" as the reference, the solutions of the other cases are compared using the error criterion of eq. (5.19). The error criterion, given in fig. 5.9, shows that with the decrease in number of cycles of the training phase, the final solution deviates further from the reference solution. However the amount of inaccuracy is very low and it is not necessary to include high number of cycles in the "training stage", as it will drastically increase the numerical cost. Hence for the next analyses described henceforth will include a "training phase" of 10 cycles.

Simulation of large number of cycles

Now the aim is to observe the damage behaviour involving large number of cycles. An adaptive temporal meshing scheme given in table 5.4 is tested on the simple structure given in fig. 5.4, which is chosen to avoid any localisation and stress raisers which may lead to regional concentration of damage.

D max m < 0.001 500 0.001 < D max m ≤ 0.01 200 0.01 < D max m ≤ 0.1 100 D max m > 0.1 50
Using term D max m that describes the maximum damage value at the end of "nodal cycle" m, the size of the time element [Θ m, Θ m+1 ] i.e. the value of Θ m+1 is determined based on table 5.4.

A total of 5 × 10 5 cycles are simulated with the load amplitude U 0 = 0.060 mm. The "training stage" is based on the calculation of initial 10 cycles, after which the reduced order basis is reused to calculate every "nodal cycle". The adaptive numerical scheme computes only 2185 "nodal cycles" and the rest is interpolated. The search direction operators are calculated only once for each "nodal cycle", and then reused during the subsequent LATIN iterations. A total of three PGD modes are obtained during the training phase of the problem, and these basis functions are not enriched during the calculation of the "nodal cycles".

The distribution of damage at the end of loading is depicted in fig. 5.12. The central region of the plate consists of a weaker material, hence the damage evolution is much higher as displayed in fig. 5.13a. For the peripheral region, the damage evolution is much lower and it also shows saturation after certain load cycles (see fig. 5.13b).

The tests performed till now is based on fig. 5.4 where the basis generated from the "training stage" is considered to be optimal and is later on re-used for the "nodal cycles" without enrichments. However the examples provided in sections 5.3.4 to 5.3.5 are on structures that induce localisation, hence the basis .13: Damage evolution in the plate for a total of 5 × 10 5 cycles obtained from "training stage" is not considered to be optimal, and the ROB is enriched by a maximum of one PGD mode per "nodal cycle".

Pre-damaged structure

The idea of this analysis is to simulate and compare the behaviour of virgin and pre-damaged materials, for which a two-dimensional notched plate is considered. The structure consists of a rectangular plate of length L = 20 mm, width W = 10 mm, and two semicircular notches of equal radii φ = 4 mm as illustrated in fig. 5.14. The structure is loaded with uniformly distributed sinusoidal displacements of the form U d = U 0 sin 2πt ∆T with time period ∆T = 1 s on the longitudinal ends for 20000 cycles. The quarter of the plate generated out of symmetry for the numerical analysis is also shown in fig. 5.14, along with the region of interest where damage localisation is anticipated.

U d U d L φ W U d
Figure 5.14: A plate with semicircular notches subjected to distributed loads For a given load amplitude U 0 = 0.0057 mm, the first analysis considered is to assume the complete structure to be virgin. A total of 7 PGD modes are calculated during the "training stage" which are then used to approximate the "nodal cycles". A total of 171 "nodal cycles" and 11 PGD modes are calculated to simulate 20000 cycles. Thereafter the Gauss points marked in red in fig. 5.15 are assigned an initial damage value of 0.005, and the problem is recalculated. A total of 186 "nodal cycles" and 13 PGD modes are needed to simulate the total behaviour, with a "training stage" of 10 cycles and 7 PGD modes.

Variable amplitude loading

The intent herein is to simulate the difference in behaviour of a given structure subjected to constant amplitude (CA) load and variable amplitude (VA) load. The effect of sequence in case VA load, i.e. lowto-high (L-H) or high-to-low (H-L) is also studied. The structure used for this analysis is a square plate of length L = 40 mm. This plate consists of an elliptical hole with semi-major axis a = 10 mm and semiminor axis b = 5 mm. The plate is subjected to uniformly distributed sinusoidal displacements of the form U d = U 0 sin 2πt ∆T with time period ∆T = 1 s on the ends, as shown in fig. 5.18 for 10000 cycles. The quarter of the plate, with the symmetric boundary conditions, along with the anticipated region of interest, is also shown in fig. 5.18.

U d U d L 2a 2b U d Figure 5
.18: A square plate with an elliptical hole subjected to distributed loads With the same material properties as before, the first part of the analysis is to simulate the complete loading using constant amplitude. The "training stage" consists of 10 cycles, which yields 8 PGD pairs. Using the adaptive temporal scheme, a total of 95 "nodal cycles" are simulated by generating 19 PGD pairs.

In the next part of the analysis, the load amplitude U 0 is considered to be uniformly increasing as shown in fig. 5.19. The "training stage" similar to the constant amplitude case consists of 10 load cycles generating 7 PGD modes. A total of 90 "nodal cycles" are generated with a total of 18 PGD modes at the end of the last "nodal cycle". It has to be mentioned that the mean value of the load amplitude for the three cases is considered to be the same. The spread of damage for these three loading cases, at certain cycles, in the region of interest is shown in fig. 5.20. Damage is obtained to be lower for constant amplitude loading as compared to the H-L loading but higher that L-H loading. The effect of sequence can be observed if the evolution of damage is plotted for the weakest GP for the three load cases (fig. 5.21). The evolution of damage for H-L loading is found to be higher than for L-H loading.

Virtual ε-N curves

To develop virtual ε-N curves, a one-dimensional structure as shown in fig. 4.1 is considered. The bar is considered to be made of a single material with the material properties given in table 5.2. A verification study is done at the beginning, which involves comparing the solutions obtained using different uniform temporal meshes with the mono-scale reference solution. A total of 2000 cycles are simulated using the mono-scale LATIN-PGD method. Thereafter the same problem is solved using different uniform temporal meshes and the comparative solution is depicted in fig. 5.22a. As the temporal mesh gets finer the solution approaches the reference solution. A more quantitative description of this fact is presented in fig. 5.22b, where the error, defined by eq. ( 5.18), increases with increase in number of cycles per time element. The CPU time however shows a drastic decrease as compared to the reference CPU time, and the relative CPU cost decreases with increase in number of cycles per time element. To build an ε-N curve, a critical damage level of 0.2 is chosen and the structure is loaded at a given strain amplitude εa. The structure is assumed to have failed when the damage value reaches 0.2 and the number of cycles needed is noted. This ordeal gives one point of the ε-N curve, which is thereby repeated for approximately 20 times with different strain amplitudes to obtain other points, which are then joined to obtain one particular ε-N curve. As far as the temporal discretisation is considered, each virtual point is obtained using a uniform temporal mesh which varies for different strain amplitudes, the lowest being 10 cycles per time element and the highest being 200 cycles per time element. The mean-strain effect on the ε-N curve is depicted in fig. 5.23a. For a negative mean strain εm, the evolution of damage will be less as the compressive part will be more than the tensile part, thereby for a given amplitude εa, more cycles will be needed to achieve a critical damage level, as compared to a zero mean-strain case. For a positive mean-strain, however, the effect is reversed, with a higher tensile part, the number of cycles needed to reach the critical damage value is less compared to a zero mean-straincase. The effect of different yield stresses is shown in fig. 5.23b, where, with decrease in yield stress the material becomes more susceptible to damage and less cycles are needed to reach the critical damage value for a given strain amplitude. Also the damage threshold considered here is directly proportional to the yield stress, hence higher σy induces lower damage. The effect of different values of initial damage D in is shown in fig. 5.23c, where with higher initial damage value, the number of cycles needed to reach the critical damage becomes less for a given strain amplitude. The effect of various moduli of elasticity is displayed in fig. 5.23d. For a higher value of the modulus of elasticity, the stress becomes higher at a given εa and generates lesser number of cycles to obtain the critical damage. Also the damage threshold chosen here is inversely proportional to the modulus of elasticity, hence higher E induces higher damage.

Concluding remark

An approach has been proposed in this chapter that can effectively handle large number of cycles. Verification studies that were performed, showed confirmation with the mono-scale LATIN method, with finer temporal discretisation providing more accuracy but consuming more CPU time. Adaptive temporal discretisation has thereby been used to obtain optimal compromise between numerical cost and accuracy. This method is effective to deal with various types of loading and initial conditions and also been used to simulate virtual ε-N curves.

In a way this method can be compared to the classical jump-cycle approach. Indeed, the "nodal cycles" used in this non-incremental framework can be compared to the cycle computed after the jumps used in an incremental framework. The finite element time description is however novel in this approach and can effectively approximate the intermediate cycles in a time element, i.e. the jumped cycles. The choice of the initial conditions for the cyclic quantities at each "nodal cycle" is similar to the jump cycle method. However, the initial condition for damage is linearly approximated in jump cycle, whereas here, a "resolvent" technique was used, which is more accurate.

Chapter 6 Conclusion and future perspective

A novel technique has been developed to include damage as an internal variable in a classical LATIN-PGD framework. The main difficulty that arose herein was the elastic state law which was no longer linear. To circumvent this difficulty, the elastic state law was solved in the local stage of the LATIN method and the stress and total strain tensor were separated into parts depending on plasticity and parts depending on damage. The separated parts should also be admissible and PGD was applied to the plastic part only. The unilateral condition of micro-defect closure effect was used in the damage description, which provided the difference between tension and compression. An on the fly PGD technique based on a greedy algorithm was used to build the ROM which satisfied the global admissibilities. A hybrid PGD method was used to solve the global step where the temporal approximations were calculated though the minimisation of a mechanical residual and the spatial approximations were obtained from classical Galerkin method. PGD-based approximation was used to approximate the quantities depending on plastic deformation and quantities depending on damage were already known from the local stage of the LATIN algorithm. For the hardening terms, no model reduction approach was utilised as they are obtained from cheap local first order ODEs. This technique was tested in simple academic examples subjected to cyclic loading and the results were qualitatively satisfactory. The test results also showed that the technique was able to handle non-proportional loading quite efficiently.

This approach was thereafter extended so that a large number of cycles could be simulated. A twoscale temporal discretisation approach was adopted, which took into account the slow evolution of the quantities of interest along the load cycles and their fast evolution within a cycle. The quantities of interest were calculated only at certain "nodal cycles" and finite element like temporal interpolation allowed to approximate the intermediate cycles. The first few cycles were calculated using mono-scale LATIN technique to obtain the ROB ("training stage" of the two-scale method), which were then utilised to obtain the approximations at each "nodal cycle". Every "nodal cycle" was subjected to a separate LATIN algorithm where the ROB of the previous "nodal cycle" was updated and enriched only if necessary. A "resolvent" technique introduced in Cognard and [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF] for visco-plasticity was used to calculate the initial conditions for damage at each "nodal cycle". This method was compared with the solution of mono-scale LATIN-PGD method and acceptable accuracy was obtained with drastic reduction in numerical expense. It was thereafter used to simulate large number of cycles and also to obtain virtual ε-N curves. An adaptive temporal scheme was also used where the sizes of the time elements are not considered to be equal but adapted according to the damage evolution. This adaptive temporal scheme was used to simulate pre-damaged structure and structure subjected to variable amplitude load.

Although, non-proportional loading can be handled, the limitation of the proposed method to include damage in a LATIN-PGD clearly is the fact that high value of damage will lead to inaccuracy. Although for metals where the critical damage is relatively low, this is not a problem, for damage computations in materials like rubber, the damage may be much higher compared to metals. Therefore, a formulation that is able to tackle high values of damage with acceptable accuracy should be developed in the future. Other 93 corresponding rows of O red tot are thereby filled with zeros to obtain Otot. Equation (A.4) can thereby be written as {ūtot} = Otot εp tot .

(A.6)

The total assembled matrix containing the derivatives of the shape functions can then be used to rewrite eq. (A.6) as

B T tot {ūtot} = B T tot Otot εp tot =⇒ {εtot} = E εp tot , (A.7)
where E is the FE operator relating the space function of plastic strain εp tot and total strain {εtot}. Now to find the FE operator C, the strain partition relationship is used, i.e.

{ε e tot } = {εtot} -εp tot =⇒ {σtot} = Ctot {εtot} -εp tot =⇒ {σtot} = Ctot (E -I) εp tot =⇒ {σtot} = C εp tot , (A.8)
where Ctot is the complete elasticity matrix for all the Gauss points, and C is the FE operator which relates the space function of plastic strain εp tot and stress {σtot}.

Appendix B Solution technique of the temporal problem

The purpose herein is to elaborate the methodology of solving the temporal problem that is defined as a minimisation of a mechanical residual in chapter 4.

The residual is defined in terms of the constitutive relation error defined in eq. (3.31). As the separated representation is only for the stress and plastic strain, only the first part of eq. (3.31) is rewritten as

e 2 CE = ∆ εp -Hσ∆σ ′ + ∆ε 2 M , (B.1) with M = H -1 σ . (B.2)
The LATIN indices i + 1 are dropped for clarity and will be sustained for the complete appendix. Considering the case where a new mode is added, i.e. Equation (B.5) leads to a strong form consisting of a second order ODE with mixed boundary conditions (see [START_REF] Passieux | Time-space radial approximation and multiscale LATIN method[END_REF]. However, solution of eq. (B.5) is calculated in a weak sense. A discontinuous Galerkin scheme of order zero is used to solve eq. (B.5), the detail of which can be found in [START_REF] Passieux | Time-space radial approximation and multiscale LATIN method[END_REF], [START_REF] Relun | Multiparametric strategy for robust design in fatigue[END_REF]. The main idea is to have piecewise constant functions for every time element and at each node there is a jump or discontinuity in the nodal unknowns (see fig.

B.1 for a particular temporal quantity λ ).

λ t t 0 t 1 t 2 t n-1 t n λ (t 0 ) λ (t - 1 ) λ (t + 1 ) λ (t - 2 ) λ (t + 2 ) λ t + n-1 λ (t - n ) λ (t + n ) Figure B
.1: Discontinuous Galerkin scheme of order zero For any time element [t n-1 , tn], λ t + n-1 = λ t - n , and the derivative λ is approximated as

λ = λ t + n = λ t + n -λ t - n ∆tn , (B.7)
with ∆tn = tn -t n-1 .

To derive the elemental equation for a particular time element [t n-1 , tn], the domain of integration is only for [t n-1 , tn] in eq. (B.5), i.e. It is needless to say that this formulation is exactly similar to first-order finite difference method, where the nodal unknowns can be expressed as

q el = λ n-1 λn , (B.14)
where λ n-1 and λn represent the nodal unknowns at time points (or nodes) t n-1 and tn respectively. The elemental equation, i.e. eq. (B.10) must be assembled for all the time elements resulting in the system of linear equations Ttotqtot = btot, (B.15)

where, Ttot is obtained from the assembly of T el , btot from the assembly of b el , and qtot will contain all the nodal unknowns. Generally, this problem is solved with the knowledge of λ 0 , i.e. λ at t 0 , which basically shows that the boundary value problem has been transformed into an initial value problem. Now, if the minimisation problem is for the update of the time functions, the linearised form invariably includes more than one degree of freedom per node. The number of degree of freedom per node is exactly the same as the number of modes that are being updated. If m represents the number of modes that are present, the corrections of the stress and the plastic strain rate can be written as,

∆ εp = m i=1 ∆ λi εp i , (B.16a) ∆σ ′ = m i=1 ∆λ i Cε p i . (B.16b)
The constitutive relation error will thereby be written as

e 2 CE = m i=1 ∆ λi εp i -Hσ m i=1 ∆λ i Cε p i + ∆ε 2 H -1 σ . (B.17)
The minimisation will henceforth be on all the m time functions, which generates the element matrix of the form

b el =                           d 1 1 ∆tn d 2 1 ∆tn . . . d m 1 ∆tn - d 1 1 ∆tn + d 1 0 - d 2 1 ∆tn + d 2 0 . . . - d m 1 ∆tn + d m 0                           , with d i 1 = Ω εp i : H -1 σ ∆ε dV,
q el =            ∆λ 1, n-1 ∆λ 2, n-1 . . . ∆λ m, n-1 ∆λ 1, n ∆λ 2, n . . . ∆λm, n            , (B.19)
and

T el =             P 1,1 P 1,2 . . . P 1,m Q 1,1 Q 1,2 . . . Q 1,m P 2,1 P 2,2 . . . P 2,m Q 2,1 Q 2,2 . . . Q 2,
P m,1 P m,2 . . . P m,m Q m,1 Q m,2 . . . Q m,m Q 1,1 Q 1,2 . . . Q 1,m R 1,1 R 1,2 . . . R 1,m Q 2,1 Q 2,2 . . . Q 2,m R 2,1 R 2,2 . . . R 2,m . . . . . . . . . . . . . . . . . . . . . . . . Q m,1 Q m,2 . . . Q m,m R m,1 R m,2 . . . R m,m             , with P i,j = a i,j 11 (∆tn) 2 , Q i,j = - a i,j 11 (∆tn) 2 + a i,j 10 ∆tn , R i,j = a i,j 11 (∆tn) 2 - a i,
λ 2 m+1 dt/ [0,T ] λ 2 1 dt > µ tol then Update m = m + 1 end Result: εp i = εp 0 + m j=1 λj εp j
Hooke's tensor in compliance form ( Ki+1 = C-1 i+1 ) at each space-time point is calculated as

Ki+1 (t, x) = Gi+1 (t, x)           1 -ν -ν 0 0 0 -ν 1 -ν 0 0 0 -ν -ν 1 0 0 0 0 0 0 1 + ν 0 0 0 0 0 0 1 + ν 0 0 0 0 0 0 1 + ν           , (D.5) with Gi+1 (t, x) = 1 E (1 -D i+1 (t, x)) if sg (t, x) ≥ 0 = 1 E (1 -hD i+1 (t, x))
if sg (t, x) < 0.

(D.6)

A singular value decomposition is performed on Gi+1 such that it is represented as

Gi+1 (t, x) ≈ q l=1 g l i+1 (t) Ḡl i+1 ( x) , (D.7)
where q is the number of pre-selected modes. From this representation of Gi+1 (t, x), the Hooke's tensor can be written as

Ki+1 (t, x) ≈ q l=1 k l i+1 (t) Kl i+1 ( x) , (D.8)
where

Kl i+1 ( x) = Ḡl i+1 ( x)           1 -ν -ν 0 0 0 -ν 1 -ν 0 0 0 -ν -ν 1 0 0 0 0 0 0 1 + ν 0 0 0 0 0 0 1 + ν 0 0 0 0 0 0 1 + ν           , (D.9)
and k l i+1 (t) = g l i+1 (t). As far as the Hooke tensor in stiffness form Ci+1 (t, x) is concerned, it is easily represented in a similar space-time separable form as

Ci+1 (t, x) ≈ q l=1 c l i+1 (t) Cl i+1 ( x) , (D.10)
where

c l i+1 (t) = 1 k l i+1 (t) , Cl i+1 ( x) = Kl i+1 ( x) -1 , ∀ l ≤ q.
(D.11) e.g. en fatigue à faible nombre de cycles, il est plutôt recommendé de considérer une approche de type ≪ strain-life ≫ comme proposé dans [START_REF] Coffin | A study of the effects of cyclic thermal stresses on a ductile metal[END_REF] et [START_REF] Manson | Behaviour of materials under conditions of thermal stress[END_REF] pour des matériaux métalliques. Similairement aux courbes S-N qui sont utilisées pour des approches de type ≪ stress-life ≫, les courbes ∆ε-N ou εa-N utilisées pour les approches ≪ strain-life ≫ (voir [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF] sont des représentations des déformations totales ∆ε ou de leur amplitude εa ou en fonction du nombre de cycles de survie à la fatigue. Les méthodes empiriques peuvent être étendues pour traiter des cas plus sophistiqués. Par exemple, l'effet dû à la contrainte moyenne peut être pris en compte (voir Gerber, 1874[START_REF] Goodman | Mechanics Applied to Engineering[END_REF][START_REF] Soderberg | Factor of safety and working stress[END_REF], ainsi que celui lié aux concentrations de contraintes dues à des entailles en utilisant la règle de Neuber (voir Neuber, 1937). La durée de vie en fatigue d'une structure soumise à des blocs de chargements de différentes amplitudes peut être étudiée en utilisant une loi d'endommagement cumulative couplée avec une courbe S-N . La loi d'accumulation peut être linéaire comme celle de Palmgren-Miner (voir Miner, 1945[START_REF] Palmgren | Die Lebensdauer von Kugellagern[END_REF] ou non-linéaires comme la loi d'endommagement proposée par [START_REF] Marco | A concept of fatigue damage[END_REF].

Les méthodes pré-citées sont phénoménologiques par nature et basées sur l'utilisation de relations empiriques provenant de résultats expérimentaux. Par exemple, elles n'offrent pas la flexibilité nécessaire pour analyser les effets dus à l'ordre dans lequel les différents niveaux de chargement sont appliqués à la structure. Pour dépasser cette limitation, les changements observés dans la structure durant le chargement cyclique peuvent être décrits par des variables internes dans le cadre de la mécanique des milieux continus. L'endommagement défini en tant que variable interne est utilisé pour quantifier la phase d'initiation de fissures macroscopiques qui visent à réduire la capacité portante du matériau (voir [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]. L'utilisation de cette approche pour prédire l'espérance de vie à la fatigue d'une structure ou d'un matériau fut introduite par [START_REF] Chaboche | A non-linear continuous fatigue damage model[END_REF], ils utilisèrent ainsi un modèle d'endommagement continu non-linéaire pour décrire les différentes phases du procédé de détérioration du matériau. Diverses modifications et développements ont été réalisés au fil des années pour améliorer la modélisation de l'endommagement par fatigue et incorporer les différents phénomènes physiques aussi précisément que possible. Par exemple, quand le chargement est important, la structure subit d'importantes déformations plastiques, menant à une durée de vie restreinte, en général inférieure à 10 5 cycles, ce cas est dénommé fatigue oligocycliques ou à faible nombre de cycles. Au contraire, lors de cas de fatigue à grand nombre de cycles, le chargement est beaucoup plus faible que la limite d'élasticité. Donc, aucune déformation plastique macroscopique n'est mise en jeu, la structure peut supporter un très grand nombre de cycles. Le modèle d'endommagement à deux échelles (voir [START_REF] Lemaitre | A two scale damage concept applied to fatigue[END_REF] est un développement majeur pour la modélisation de la fatigue à grand nombre de cycles, il représente le comportement macroscopique élastique alors que l'endommagement est évalué seulement à l'échelle microscopique.

Pour étudier la propagation d'une macro-fissure, le cadre de la mécanique de la fracture est utilisé connaissant la pré-existence d'une fissure dans le matériau. La loi la plus traditionnelle pour décrire la croissance des fissures est la loi de Paris-Erdogan (voir Paris and Erdogan, 1963).

Plus d'informations sur les techniques classiques et les développements récents au sujet de la modélisation et des simulations pour des cas de fatigue peuvent être trouvées par exemple dans [START_REF] Cui | A state-of-the-art review on fatigue life prediction methods for metal structures[END_REF]. Une approche de description continue de l'endommagement est ici utilisée pour modéliser le comportement à la fatigue de manière flexible. Ceci permet de considérer la chronologie des différents cycles ou les effets d'inertie dus à des cycles à hautes fréquences. Cependant, cette approche peut mener à des coûts de calcul prohibitifs. Des développements ont été réalisés dans le cadre de techniques numériques nouvelles et robustes pour réduire le coût de calcul. Par exemple, ce problème peut être surmonté en utilisant des techniques de réduction d'ordre de modèles.

E.1.2 Techniques de réduction de modèles pour des calculs en fatigue

La réduction d'ordre de modèle est une famille de stratégies numériques qui a montré son efficacité pour de nombreux problèmes à grande dimension tels que les études paramétriques ou les calculs en temps réel (voir Cueto et al., 2014). La solution est recherchée en résolvant le problème de Galerkin dans une base d'ordre réduit, c'est-à-dire dont la dimension est beaucoup plus faible que la taille du modèle originel à grande dimension. L'utilisation d'une décomposition orthogonale aux valeurs propres (Proper orthogonal decomposition : POD) dans le champ de la mécanique est basée sur la construction initiale d'une base réduite désignée ici comme base POD, à partir des snapshots de la solution de l'étape d'apprentissage. Ensuite, cette base est utilisée pour résoudre le problème d'intérêt à coût réduit (Néron andLadevèze, 2010, Ryckelynck, 2009). La POD fournit une base réduite optimale au sens des moindres carrés pour la solution du problème d'ordre complet. En utilisant une décomposition propre généralisée (Proper Generalised Decomposition : PGD), le problème est également résolu sur une base réduite, mais cette base est définie au cours du calcul en fonction du problème considéré par un algorithme glouton [START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF]Chinesta et al. ( , 2014b)), [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF]. Cette technique ne requiert aucune phase d'apprentissage et évalue la qualité du modèle réduit au cours des itérations enrichissant la base réduite si besoin.

La méthode LATIN (voir Ladevèze, 1999) est une approche numérique dans laquelle une approximation de la solution est recherchée sur l'ensemble du domaine spatio-temporel à chaque itération, l'équilibre global traité en tant que problème linéarisé et le comportement non-linéaire du matériau sont considérés séparément. Donc, cette approche offre un cadre approprié pour inclure des techniques de réduction d'ordre de modèle même pour des problèmes a priori non linéaires. Elle a été développée initialement pour résoudre des problèmes plastiques et visco-plastiques et a montré une réduction drastique du coût numérique en comparaison aux approches traditionnelles aussi bien dans son cadre initial que dans de nombreuses applications dans les dernières décennies (voir [START_REF] Ladevèze | PGD in linear and nonlinear computational solid mechanics[END_REF][START_REF] Relun | A model reduction technique based on the PGD for elasticviscoplastic computational analysis[END_REF]. Une nouvelle version est proposée dans ces travaux de thèse pour étudier des problèmes (visco-)plastiques avec endommagement par fatigue prenant en compte l'effet de fermeture des micro-défauts.

Pour des simulations des fatigue appliquées à un grand nombre de cycles, des stratégies numéri-ques particulières sont nécessaires pour éviter le coût numérique dû à l'étendue du domaine temporel, il peut être remarqué que dans le cadre de simulations de fatigue, ce n'est pas un calcul en temps réel qui est recherché, mais un calcul à temps réduit.

E.1.3 Schémas temporaux efficaces pour des calculs de fatigue

Pour des cas de chargements de fatigue à grand ou très grand nombre de cycles ou dans des cas de fatigue combinée, c'est-à-dire pour lesquels le chargement est une combinaison de chargement à large amplitude-basse fréquence et faible amplitude-haute fréquence (voir Suresh, 2001), le coût numérique dû aux schémas d'intégration temporelle peut être extrêmement élevé. Ce coût peut être réduit en bénéficiant de la nature cyclique du chargement pour ne pas calculer explicitement l'ensemble des cycles.

La procédure par saut de cycles est une technique très robuste (voir [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF] qui évite le calcul de blocs entiers de cycles. Après avoir calculé en détail un ensemble de cycles, le résultat est utilisé pour établir une tendance d'évolution et extrapoler les quantités d'intérêt durant la suite du chargement. Ainsi, un état initial extrapolé est utilisé pour calculer un cycle futur et extrapoler un bloc de cycles suivants, ainsi de suite jusqu'à ce que l'ensemble du domaine temporel soit simulé. Une fonction de contrôle a été proposée de sorte que la longueur des sauts de cycles soit adaptée pour garantir une qualité d'approximation acceptable (voir [START_REF] Cojocaru | A simple numerical method of cycle jumps for cyclically loaded structures[END_REF]. Cette méthode est très appropriée et effiace pour des systèmes aux comportement quasi-linéaires. Cependant, la fonction de contrôle permet également de considérer des comportements non-linéaires, les sauts de cycles sont alors automatiquement raccourcis voire abandonnés.

Pour les cas de fatigue combinée, les techniques d'homogénéisation temporelles peuvent être utilisées. Elles sont basées sur une hypothèse de séparation d'échelle de temps entre une échelle de temps large associée avec le chargement à basse fréquence et une échelle de temps court associée avec le chargement à haute fréquence. Le rapport entre ces deux temps est supposé suffisamment faible pour considérer que les deux échelles sont indépendantes. Ainsi, le comportement dû à la fréquence élevée est homogénéisé et le problème est résolu pour une discrétisation temporelle adaptée au chargement à basses fréquences. Développée initialement pour des cas quasi-statiques et des comportement plastiques (voir [START_REF] Guennouni | Réponse homogénéisée en temps de structures sous chargements cycliques[END_REF], cette méthode a été étendue plus récemment pour des cas avec endommagement ou effets dynamiques (voir Oskay and Fish, 2004b).

Dans le cadre de la méthode LATIN, comme l'ensemble du domaine temporel est étudié à chaque itération, une double échelle temporelle dans le cadre d'une discrétisation de type ≪ éléments finis ≫ a été proposée pour des problèmes de visco-plasticité (voir [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF]Ladevèze, 1993, Ladevèze, 1999). Seulement quelques cycles d'intérêt, appelés ≪ noeuds temporels ≫ ou ≪ cycles nodaux ≫, sont calculés en utilisant un schéma temporel traditionnel. L'évolution des quantités d'intérêt entre deux noeuds est interpolée en utilisant des fonctions de forme de manière similaire à la méthode des éléments finis spaciaux.

Les méthodes empiriques étant basées sur des expériences complexes et coûteuses, le but de ce projet est de développer une approche numérique permettant de prévoir le comportement en fatigue basé sur la mécanique des milieux continus. Pour faire face au coût numérique induit par ces modèles, un schéma numérique basé sur une réduction d'ordre de modèle approprié est proposé qui permet de simuler des tests virtuels pour de très grands nombres de cycles avec un coût raisonnable. Le schéma numérique proposé est flexible quant au modèle phénoménologique d'endommagement considéré. Le cadre général ainsi que les équations considérées sont résumées dans la section E.2. Ensuite, la réduction d'ordre de modèle basée sur la méthode LATIN-PGD est introduite dans la section E.3. Finalement, le traitement numérique du problème temporel qui permet de réduire le coût numérique est détaillé dans la section E.4.

E.2 Approche basée sur la mécanique des milieux continus

f d • δ ˙ u dV dt + ∂Ω 2 ×[0,T ] F d • δ ˙ u dSdt = 0, (
δσ : ε dV dt + ∂Ω 1 ×[0,T ] δσ • n • ˙ U d dSdt = 0, (E.2)
où ε est le tenseur de déformations totales qui est cinématiquement admissible et peut être séparé en une partie élastique ε e et une partie plastique ε p additivement. δσ est un champ statiquement admissible.

Les propriétés mécaniques des matériaux sont décrites par un ensemble de relations constitutives. Les équations d'état pour les matériaux élasto-plastiques soumis à un endommagement unilatéral sont obtenues à partir de la fonction d'énergie libre (voir [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF]) Une réduction de modèle est utilisée par l'intermédiaire d'une séparation des dépendances en temps et en espace. L'approche LATIN-PGD permet pour les structures ayant un comportement visco-plastique de définir les deux quantités d'intérêt, qui sont les contraintes et les déformations plastiques, en utilisant une unique base temporelle. Pour examiner des cas incluant de l'endommagement, une nouvelle quantité d'intérêt, la déformation élastique, est introduite. Le champ solution est noté s = εp , ε e , Ẋ, Ḋ, σ, Z, Y où X représente l'ensemble des variables d'écrouissage (cinématique et isotropique) et Z est la variable conjuguée à X. La loi d'état non-linéaire peut être incluse dans l'étape globale conduisant à une séparation de variables utilisant différentes fonctions temporelles pour les contraintes et les déformations. Sinon, comme présenté ici, la non-linéarité due à la loi d'état peut être considérée à l'étape locale et les contraintes peuvent être séparées en deux contributions, l'une définie à l'étape locale et l'autre écrite en tant que décomposition en formes dépendant indépendemment du temps et de l'espace utilisant les mêmes fonctions temporelles que pour la déformation plastique.

ε e ij = 1 + ν E σ + ij 1 -D + σ - ij 1 -hD - ν E σ kk 1 -D + -σ kk 1 -hD δ ij , (E.3a) β ij = Qα ij , (E.3b) R = R∞ (1 -exp (-γr)) , (E.3c) Y = 1 + ν 2E σ + ij σ + ij (1 -D) 2 + h σ - ij σ - ij (1 -hD) 2 - ν 2E σ kk 2 (1 -D) 2 + h -σ kk 2 (1 -hD) 2 . (E.
εp = k f p n +   3 2 ̺ 3 2 [̺ : ̺]   1 1 -D , (E.4a) α = -k f p n +   - 3 2 ̺ 3 2 [̺ : ̺] + aQ -1 β   , (E.4b) ṙ = k f p n + , (E.4c) Ḋ = k d f d n d + , (E.4d) avec ̺ = σ D 1 -D -β et σ D étant
f p = 3 2 [̺ : ̺] + a 2 β : Q -1 β -R -
L'algorithme est initialisé par la résolution d'un problème considérant les conditions aux limites du problème d'intérêt mais un comportement du matériau parfaitement élastique quelles que soient les conditions de chargement. Puis, des corrections liées à la plasticité et à l'endommagement sont ajoutées successivement à la solution élastique à chaque itération. L'ensemble des équations est divisé en deux sousgroupes, l'un comprenant les équations globales et linéaires alors que le second comprend les équations locales et non-linéaires. Une itération LATIN se compose de deux étapes :

- L'estimation de la solution est recherchée alternativement dans les deux variétés jusqu'à atteindre la convergence. À partir de l'estimation de la solution dans une étape, l'approximation est recherchée dans la variété suivante en utilisant certains opérateurs linéaires appelés opérateurs de direction de recherche.

E.3.1 Étape locale

À l'étape locale, les équations d'évolution pour les variables internes sont résolues, elles sont locales en espace et non-linéaires. La loi d'état élastique étant non-linéaire est également traitée durant cette étape. À partir de la solution s i ∈ A d à l'itération LATIN i, l'estimation ŝi+1/2 ∈ Γ est réalisée de telle sorte que les directions de recherche locale soient satisfaites

     εp i+1/2 -εp i -Xi+1/2 -Ẋi εe i+1/2 -ε e i Di+1/2 -Ḋi      + B +    σi+1/2 -σ i Ẑi+1/2 -Z i σi+1/2 -σ i Ŷi+1/2 -Y i    = 0. (E.8)
Ici, la direction de descente est notée B + . En suivant l'idée proposée par Ladevèze (1999), la direction de recherche est considérée verticale de telle manière que

B + -1 = 0. (E.9)
La solution de l'équation définissant la direction de recherche éq. (E.8) ainsi que des équations d'évolution éq. (E.4) et de la loi élastique non-linéaire éq. (E.3a) constitue l'estimation ŝi+1/2 . À partir de cette approximation, la solution de l'étape globale s i+1 est recherchée.

   εp i+1 -εp i+1/2 -Ẋi+1 -Xi+1/2 ε e i+1 -εe i+1/2    -B - σ i+1 -σi+1/2 Z i+1 -Ẑi+1/2 σ i+1 -σi+1/2 = 0 (E.10a) Ḋi+1 -Di+1/2 -b -Y i -Ŷi+1/2 = 0, (E.10b) où B -= H -0 0 C -1 . L'
∆σ i+1 = σ i+1 -σ i , ∆ε e i+1 = ε e i+1 -ε e i et ∆ εp i+1 = εp i+1 -εp i . (E.13)
Les corrections en termes de contraintes et de déformations totales à l'étape globale de l'itération i+1 sont ainsi séparées en parties dépendant de la déformation plastique (∆σ ′ i+1 , ∆ε ′ i+1 ) et de l'endommagement (∆σ i+1 , ∆ε i+1 ) de telle sorte que

∆σ i+1 = ∆σ ′ i+1 + ∆σ i+1 , (E.14a) ∆ε i+1 = ∆ε ′ i+1 + ∆ε i+1 . (E.14b)
À partir de ces séparations et de l'équation de direction de recherche ainsi que de la relation additive de décomposition des déformations, il peut être établi que

∆σ ′ i+1 + ∆σ i+1 = C ∆ε ′ i+1 -∆ε p i+1 + C ∆ε i+1 -∆ε R i+1 , (E.15)
où ∆ε R i+1 peut être interprété comme une déformation résiduelle provenant de la loi d'état non-linéaire à l'itération i + 1. Les termes correctifs ∆σ i+1 et ∆ε i+1 sont ainsi obtenus à partir de l'équation d'équilibre, directement. D'autre part, si seulement la partie plastique est considérée, la direction de recherche peut être écrite ∆ εp i+1 -Hσ∆σ ′ i+1 + ∆i+1 = 0, (E.16) où ∆i+1 est un terme correctif plastique provenant de l'étape locale et Hσ représente la partie découplée de H -qui rellie la contrainte aux taux de déformation plastique.

Les termes de correction relatifs au comportement plastique ∆σ ′ i+1 et ∆ε ′ i+1 sont alors écrits dans une forme séparable en utilisant la PGD.

E.3.2.1 Séparation de variables

La Décomposition Propre Généralisée (PGD) est une technique de réduction d'ordre de modèles flexible, qui n'est pas basée sur une phase d'entraînement. Étant donné qu'à chaque itération LATIN, les quantités d'intérêt sont approximées sur l'ensemble du domaine spatio-temporel par une forme linéaire de l'équilibre mécanique, l'utilisation de la PGD couplée avec la méthode LATIN est opportune. Toute fonction dépendant de plusieurs variables indépendantes peut s'écrire sous la forme d'une somme infinie de produits de fonctions dépendantes d'une seule variable (voir Ladevèze, 1999(voir Ladevèze, , 2014)), la PGD recherche des approximations des quantités d'intérêt sous la forme de sommes finies de produits de fonctions définies dans des espaces de plus petite dimension que la fonction originelle en utilisant un algorithme glouton. Donc, cette approximation inclut une erreur due à la troncature de la série de formes séparables. Ici la déformation plastique et la part des contraintes relatives à la déformation plastique sont approximées par Cet algorithme par point fixe entre les problèmes spaciaux et temporels converge rapidement. Une fois le tenseur des contraintes connu à l'itération i + 1, le taux de libération d'énergie de déformation élastique relative à l'endommagement est calculée à partir de l'éq (E.3d).

E.3.3 Relaxation du champ solution et critère de convergence

Pour garantir la convergence de l'algorithme LATIN-PGD, l'étape globale requiert l'utilisation d'une relaxation. Soit si+1 la solution obtenue à la fin de l'étape globale i + 1, préalablement notée s i+1 , une nouvelle solution s i+1 est calculée

s i+1 = ϕs i+1 + (1 -ϕ) s i (E.
χ (t = θ k + τ k ) = Θ m+1 -θ k Θ m+1 -Θ m χ (Θ m + τ m) + θ k -Θ m Θ m+1 -Θ m χ (Θ m+1 + τ m+1 ) , ∀t ∈ [Θ m, Θ m+1 + ∆T ] ,
(E.28) avec τ m = τ m+1 = τ k appartenant à [0, ∆T ], la période étant constante.

Les quelques premiers cycles sont calculés classiquement pour construire une base réduite initiale (cette phase peut également être interprétée comme une ≪ étape d'apprentissage ≫). Ensuite, les cycles nodaux sont calculés progressivement. Le dernier cycle calculé classiquement défini sur le domaine [Θ 0 , Θ 0 + ∆T ] devient le cycle nodal 0 et à partir de sa connaissance, l'idée est de calculer le cycle nodal 1 défini sur [Θ 1 , Θ 1 + ∆T ]. Ensuite, connaissant le cycle nodal 1, le cycle nodal 2 est calculé, et ainsi de suite. Cette procédure est suivie jusqu'à ce que le dernier cycle nodal soit calculé.

E.4.1 Initialisation

L'initialisation des quantités d'intérêt au cycle nodal m + 1 à partir de leurs valeurs au cycle nodal m dépend des quantités qui sont considérées. La contributon de l'écrouissage isotropique est ici négligée. Les quantités qui sont cycliques, à savoir les contraintes, les déformations élastiques et les variables cinématiques, sont dupliquées à partir du cycle nodal m. Une transformation est considérée de telle sorte qu'elles deviennent périodiques et que la continuité avec le cycle m soit ainsi préservée. Les fonctions temporelles représentant les déformations plastiques {λ j } m j=1 sont également dupliquées de manière similaire. Les quantités non-cycliques, à savoir ici uniquement l'endommagement, sont initialisées comme constantes sur le cycle nodal m+1 avec la valeur obtenue à Θ m +∆T . Le taux de restitution de l'énergie de déformation liée à l'endommagement est calculé à partir de l'endommagement et du tenseur des contraintes. Le champ solution pour le cycle nodal m + 1 est calculé ultérieurement en utilisant les deux étapes de l'algorithme LATIN-PGD jusqu'à convergence.

E.4.2 Étape locale

E.4.3 Étape globale

Les modes spatiaux calculés pour les cycles initiaux sont réutilisés alors que les fonctions temporelles des modes PGD sont recalculées pour le cycle nodal m + 1. L'initialisation des fonctions temporelles {λ j } m j=1 est telle qu'elle assure la continuité des quantités à la transition entre les cycles nodaux m et m + 1. Puis, les corrections des fonctions temporelles {∆λ j } m j=1 sont calculées en résolvant l'éq. (E.19) avec des conditions initiales nulles. Si nécessaire, les bases PGD peuvent également être enrichies comme détaillé dans la section E.3.2.3. Concernant les variables cinématiques, qui sont également cycliques, elles peuvent être traitées de manière similaire. Leur initialisation est pareillement réalisée de manière à garantir la continuité vis-à-vis du cycle nodal m, et l'éq. (E.12) est résolue sur le cycle nodal m + 1. Ainsi les quantités sont calculées en termes de correction en résolvant l'équation différentielle ordinaire éq. (E.11) avec des conditions initiales nulles.

Cette approche LATIN-PGD à deux échelles de temps permet de bénéficier d'une réduction drastique de coût numérique en comparaison avec la technique LATIN-PGD traditionnelle. Plus le niveau de chargement est élevé, plus l'endommagement est important ainsi que l'erreur relative due aux maillages grossiers. Ainsi, il apparaît qu'un maillage adaptatif serait approprié pour optimiser le temps de calcul sans sacrifier la qualité de l'approximation numérique, comme proposé dans la table E.2.

E.4.4 Exemples numériques

E.5 Conclusion

Une approche numérique innovante a été introduite ici pour le calcul de l'endommagement par fatigue. Une technique non-incrémentale a été utilisée comme cadre numérique et le coût numérique a été réduit par l'utilisation de la PGD qui transforme le problème à grande dimension en un problème à dimension réduite. Une approche multi-échelle en temps a été proposée pour étendre la méthode LATIN-PGD à la simulation de l'évolution de l'endommagement lors de chargement de fatigue. Son efficacité numérique a été démontrée et ouvre une voie pour des simulations complexes basées sur la mécanique des milieux continus pour la fatigue à très grand nombre de cycles.

Stichworte:

  LATIN Methode, Schädigungsmechanik, (Visko)Plastizität, Metallermüdung, Modellreduktion, PGD, Mehrskalenansatz in der Zeit E.5 Discrétisation du domaine temporel utilisant deux discrétisations en temps : θ le long des cycles et τ au cours d'un cycle, p étant le nombre de cycles pour l'élément temporel [Θ m, Θ m+1 ] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.6 Plaque rectangulaire avec différentes limites d'élasticité . . . . . . . . . . . . . . . . . E.7 Erreur et temps de calculs relatifs pour différents maillages temporels en comparaison avec un calcul LATIN mono-échelle . . . . . . . . . . . . . . . . . . . . . . . . . . . E.8 Calculs pour différents niveaux de chargement . . . . . . . . . . . . . . . . . . . . . . E.9 Distribution et évolution de l'endommagement . . . . . . . . . . . . . . . . . . . . . . E.10 Une plaque carrée avec un trou elliptique soumise à des chargements répartis . . . . . . E.11 Distribution de l'endommagement pour certains cycles de chargement dans la région d'intérêt pour la plaque avec un trou elliptique . . . . . . . . . . . . . . . . . . . . . . E.12 Évolution de l'endommagement en fonction du nombre de cycles de chargement pour le point de Gauss le plus vulnérable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.13 Courbe ε-N virtuelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of tables 4.1 The set of constitutive relations and admissibility conditions that need to be satisfied . . 4.2 Material properties for Cr-Mo steel at 580 • C . . . . . . . . . . . . . . . . . . . . . . 5.1 Constitutive relations for only kinematic hardening . . . . . . . . . . . . . . . . . . . 5.2 Material properties for Cr-Mo steel at 25 • C . . . . . . . . . . . . . . . . . . . . . . . 5.3 Different temporal discretisations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Adaptive temporal discretisation based on damage level . . . . . . . . . . . . . . . . . E.1 Différentes discrétisations temporelles . . . . . . . . . . . . . . . . . . . . . . . . . . E.2 Discrétisation temporelle adaptative basée sur le niveau d'endommagement . . . . . . . vi List of abbreviations

  Figure 1.2: Different types of constant amplitude loading

  Figure 1.6: Constant life diagram

Figure 1 . 8 :

 18 Figure 1.8: Comparison between LDR and NLDR

Figure 1

 1 Figure 1.10: Crack growth rate curve showing different regions

  Figure 1.11: Cyclic elasto-(visco)plastictic behaviour

  Figure 2.2: Stress-strain relationship for hardening and perfect plasticity

Figure 2

 2 Figure 2.4: RVE showing virgin, damaged, and equivalent states

Figure 2 . 5 :

 25 Figure 2.5: Closure of micro-defects during compression

Figure 3

 3 Figure 3.1: Newton-Raphson algorithm at time step n + 1 (de Souza Neto et al., 2011)

  Figure 3.2: SVD approximations of a given solution field

  Figure 3.3: Decrease of error and singular values with SVD modes

  is suggested for more information and comparison between SVD, PCA, and KLE. It can be depicted that the spatial basis vectors obtained through SVD and PCA are identical, as for example the first four spatial basis vectors of the reference solution field U ref of fig. 3.2 are shown in fig. 3.4.

Figure 3

 3 Figure 3.4: First four spatial basis vectors using SVD (-) and KLE (--)

  (3.21) Appropriate separated weak forms similar to eqs. (3.18) and (3.19) can be obtained. Interested readers can consult Chinesta et al. (2014b) for an introduction to PGD technology, and Chinesta et al. (

  development thereafter is straightforward, i.e. this separable representation when introduced in the weak form eq. (3.28), yields separate equations for space and time which are solved iteratively.

  40) is sought at time t n+1 , the goal is achieved iteratively similar to LATIN method. The iterations can be divided into, -local stage: Seeking ŝ k n+1 from s k n+1 such that constitutive material behaviour i.e. σ k n+1 = C ε k n+1 is satisfied. The solution set must also satisfy the "search direction" equation, given by ε k n+1 = ε k n+1 .

Figure 3

 3 Figure 3.6: Similar iterative behaviour of Newton-Raphson algorithm as compared to LATIN method.

  ∆σ i+1 : ε (δ u) dV dt = 0, (4.34) with ∆σ i+1 = C ∆ε i+1 -∆ε R i+1 . (4.35) Furthermore, considering only the plastic part of the search direction, eq. (4.11a) is re-formulated as ∆ εp i+1 -Hσ∆σ i+1 + ∆ε i+1 = 0, (4.36) with ∆ε i+1 = Hσ σi+1/2σ i -εp i+1/2 -εp i . (4.37)

  ∆ ε′ : δσ dV dt = 0. (4.60) Considering δσ = λ δ σ, eq. (4.60) can be written as Ω ∆ ε′ λ : δ σ dV = 0, (4.61) ∀ δ σ which is statically admissible to zero, where = [0,T ] dt and λ being the known function of time. Using eqs. (4.54) and (4.59), eq. (4.61) can be written as Ω Hσλ 2 σ + λ λ C -1 σ -∆ε λ : δ σ dV = 0, (4.62)

  .68)where ϕ is a relaxation parameter. The convergence of the iterative algorithm is determined by a relative LATIN indicator. This indicator is basically the distance between the local solution and the global solution

Figure 4

 4 Figure 4.1: A bar in traction

  of pairs m / number of iterations i LATIN indicator ξ indicator w.r.t. number of pairs indicator w.r.t. number of iterations

Figure 4 . 2 :

 42 Figure 4.2: Evolution of the LATIN indicator with respect to the number of PGD pairs or LATIN iterations for the bar problem

Figure 4

 4 Figure 4.3: Space-time modes needed to approximate εp (x, t) in the bar under cyclic loading at convergence

  Figure 4.6: An "L" shaped structure subjected to a concentrated load

  Figure 4.7: Evolution of the LATIN indicator with respect to the number of PGD pairs or LATIN iterations for the "L" shaped structure

Figure 4 . 9 :Figure 4 D

 494 Figure 4.9: Distribution of residual von Mises stress in the "L" shaped structure at t = T

  Figure 4.12: A plate with a hole subjected to distributed loads and the symmetric part considered for analysis

Figure 4

 4 Figure 4.13: FE mesh depicting the virgin (•) and pre-damaged ( * ) Gauss points

Figure 4 .Figure 4

 44 Figure 4.14: Evolution of the LATIN indicator with respect to the number of PGD pairs or LATIN iterations for the virgin structure

Figure 5 . 1 :

 51 Figure 5.1: Discretisation of the time domain using two time scale discretisations: θ along the cycles and τ within a cycle, p being the number of cycles for the time element [Θ m, Θ m+1 ]

Figure 5

 5 Figure 5.2: Pictorial representation of the two-time scale scheme

Figure 5

 5 Figure 5.4: Rectangular plate with different yield stresses

Figure 5 Figure 5

 55 Figure 5.5: Relative accuracy and saved time for different temporal meshes for the structure without damage

  Figure 5.7: Comparison with the size of time elements relative to 50 cycles per time element solution

  Figure 5.9: Variation of accuracy with respect to the number of cycles of the "training stage"

Figure 5 .

 5 Figure 5.10: Distribution of damage after end of loading, i.e. after 1 × 10 5 cycles

Figure 5 .

 5 Figure 5.12: Damage distribution at the end of loading, i.e. after 5 × 10 5 cycles

Figure 5 .

 5 Figure 5.15: FE discretisation showing virgin (•) and pre-damaged ( * ) Gauss points

Figure 5

 5 Figure 5.19: Variation of amplitude with number of load cycles for three different load cases

Figure 5

 5 Figure 5.21: Evolution of damage with respect to number of load cycles for the weakest GP

  Figure 5.22: Comparison with different sizes of time elements relative to mono-scale solution

  Figure 5.23: ε-N curves

∆

  εp = λ εp , (B.3a) ∆σ ′ = λ Cε p , (B.3b) with the index of the number of modes being dropped for convenience. Equations (B.2) and (B.3) can be used to rewrite eq. (B.1) as e 2 CE = λ εp -Hσλ Cε p + ∆ε 2 Hσλ Cε p + ∆ε : H -1 σ λ εp -Hσλ Cε p + ∆ε dV dt, (B.4) with εp being the known space function. Equation (B.4) needs to be minimised with respect to the time function λ , which gives the following weak form, [0,T ] δ λ a 11 λ -a 10 λ + d 1 dt -[0,T ]δλ a 10 λ -a 00 λ + d 0 dt = 0,

LFigure

  Figure E.1: Problème de référence dans le domaine Ω

E

  .3.2 Étape globale incluant une réduction d'ordre de modèle À l'étape globale, la solution s i+1 ∈ A d satisfait les lois d'état, les conditions d'admissibilité et les directions de recherche descendantes :

  ∆ ε′ : δσ dΩ dt = 0, (E.21) et la condition d'admissibilité statique définie par l'éq. (E.1) telle que ∀δ u cinématiquement admissible à temporelle λ m+1 est résolue similairement à l'étape d'apprentissage par minimisation d'un résidu mécanique λ m+1 = arg min

  Figure E.3: Distribution de l'endommagement dans la structure en L et évolution de l'endommagement au point de Gauss le plus vulnérable

  Durant l'étape locale, toutes les quantités d'intérêt excepté l'endommagement ne nécessitent aucune intégration temporelle et peuvent être estimées directement. L'unique difficulé réside dans l'intégration de Ḋ pour obtenir la variable d'endommagement. Pour intégrer l'éq. (E.4d) sur le cycle nodal m + 1, il est nécessaire de connaître la condition initiale à θm. Considérons une équation générale ordinaire de premier ordre (voir[START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF] dχ dt + κχ = υ (E.29) définie sur l'ensemble du domaine temporel avec κ et υ des quantités temporelles connues. L'idée est de calculerχ (Θ m+1 ) à partir de χ (Θ m). L'élément temporel [Θ m, θ m+1 ] est discrétisé en certaines instances θ k de telle sorte que θ k = Θ m + k∆T , avec k = 0, 1, 2, • • • , p -2, où p est le nombre de cycles dans l'élément temporel [θ m-1 , θm]. Ainsi, on obtient θ 0 = Θ m et θ p-1 = Θ m+1 . Connaissant χ (θ k ), l'éq. (E.29) peut être résolue pour obtenir χ (θ k+1 ) sous la formeχ (θ k+1 ) = χ (θ k+1 , θ k ) + ℜ (θ k+1 , θ k ) χ (θ k ) (E.30)où χ représente la solution homogène de l'ODE et ℜ représente l'opérateur ≪ résolvant ≫ (voir[START_REF] Bhattacharyya | A model reduction technique in space and time for fatigue simulation[END_REF][START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -new approaches and non-incremental methods of calculation[END_REF]. Le défi désormais est de calculer χ et ℜ avec un coût numérique réduit. La manière la plus simple consiste à évaluer les quantités seulement aux cycles nodaux m et m + 1 puis à utiliser une interpolation linéaire pour estimerχ (θ k+1 , θ k ) et ℜ (θ k+1 , θ k ), c'est-à-dire χ (θ k+1 , θ k ) = Θ m+1 -θ k Θ m+1 -Θ m χ (Θ m + ∆T, Θ m) + θ k -Θ m Θ m+1 -Θ m χ (Θ m+1 + ∆T, Θ m+1 ) , (E.31a) ℜ (θ k+1 , θ k ) = Θ m+1 -θ k Θ m+1 -Θ m ℜ (Θ m + ∆T, Θ m) + θ k -Θ m Θ m+1 -Θ m ℜ (Θ m+1 + ∆T.Θ m+1 ) . (E.31b)Ces valeurs interpolées peuvent être utilisées pour écrire l'éq. (E.30) commeχ (θ k+1 ) = ν m (θ k ) [ χ (Θ m + ∆T, Θ m) + ℜ (Θ m + ∆T, Θ m) χ (θ k )] + ν m+1 (θ k ) [ χ (Θ m+1 + ∆T, Θ m+1 ) + ℜ (Θ m+1 + ∆T, Θ m+1 ) χ (θ k )] . (E.32) À partir de χ (θ 0 = Θ m), χ (θ p-1 = Θ m+1 ) estcalculé progressivement, et l'éq. (E.29) est résolue pour le cycle nodal m + 1. L'évolution de l'endommagement définie par l'éq. (E.4d) est estimée pour le cycle m + 1 en utilisant la technique précédemment définie.

LFigure E. 6 :

 6 Figure E.6: Plaque rectangulaire avec différentes limites d'élasticité

Figure E. 7 :

 7 Figure E.7: Erreur et temps de calculs relatifs pour différents maillages temporels en comparaison avec un calcul LATIN mono-échelle

  Figure E.8: Calculs pour différents niveaux de chargement

Table 4 .

 4 2: Material properties for Cr-Mo steel at 580 • C

	Material parameter		Value
	E		134000 MPa
	ν		0.3
	R∞		30 MPa
	γ		2 MPa
	Q		5500 MPa
	a		250 MPa
	k D	2.778 MPa -n d s -1
	n D		2
	K		1220 MPa s 1/n
	n		2.5
	k		K -n
	h		0.2
	σy	mat. 1	mat. 2	mat. 3
		80 MPa	82.5 MPa	85 MPa

Table 5

 5

	.1: Constitutive relations for only kinematic hardening
	viscoplastic
	yield function

Table 5 .

 5 2: Material properties for Cr-Mo steel at 25 • C

	Material parameter	Value
	E	199740 MPa
	ν		0.3
	Q	20602 MPa
	a	141.4 MPa
	k D	0.128 MPa -n d s -1
	n D		2
	K	195 MPa s 1/n
	n		12.5
	k		K -n
	h		0.2
	σy	189 MPa
	Table 5.3: Different temporal discretisations
	Temporal	Number of cycles	Total number
	discretisation	per time element	of time elements
	I	50	20
	II	100	10
	III	200	5
	IV	500	2
	increases er, although this increase is very less and even 500 cycles per time element provides results of
	appreciable accuracy (see fig. 5.5). The relative CPU time is estimated using t cal ob /t cal ref , where t cal ref is
	the reference mono-scale computation time and t cal

Table 5

 5

	.4: Adaptive temporal discretisation based on damage level
	Damage range	Number of cycles per time element

  λ a 11 λ -a 10 λ + d 1 dt -It must be mentioned that, the quantities a 11 , a 10 , a 00 can be time-dependent and the quantities d 1 , d 0 are most definitely time-dependent. Therefore, in eq. (B.10) average values of these quantities are used, which obviously vary for different time elements.

	Using the discontinuous Galerkin formulation, eq. (B.8) is rewritten as
	δλ t + n	a 11 (∆tn) 2 -	a 10 ∆tn	-	a 10 ∆tn	+ a 00 λ t + n + -	a 11 (∆tn) 2 +	a 10 ∆tn	λ t -n	(B.9)
	+	d 1 ∆tn	-d 0 -δλ t -n		a 11 (∆tn) 2 -	a 10 ∆tn	λ t + n -	a 11 (∆tn) 2 λ t -n +	d 1 ∆tn	= 0
	which gives the following linear system of equations
										T el q el = b el ,	(B.10)
	with
						q el =	λ t -n λ t + n	=	λ t + n-1 n λ t +	,	(B.11)
												d 1	
								b el =	 -	∆tn ∆tn + d 0 d 1	 ,	(B.12)
	and			T el =	  -	a 11 (∆tn) 2 a 11 (∆tn) 2 +	a 10 ∆tn	-a 11 (∆tn) 2 -a 11 (∆tn) 2 + a 10 ∆tn -∆tn a 10 ∆tn a 10	+ a 00	  .	(B.13)
					tn						tn
											δλ a 10 λ -a 00 λ + d 0 dt.	(B.8)
				t n-1						t n-1

δ

  m

	. . .	. . .	. . .	. . .	. . .	. . .	. . .	. . .

  This elemental form needs to be assembled to achieve the form of eq. (B.15) which can then be solved using the initial values {∆λ 1, 0 , ∆λ 2, 0 , • • • ∆λ m, 0 }, to obtain qtot.

	λj	εp j + λm+1	εp m+1 -	m j=1	εp m+1 , εp j εp j , εp j	εp j + λm+1	m j=1	εp m+1 , εp j εp j j , εp	(B.20) εp j ′
	with
	j 10 ∆tn j dV, a 00 = -a i,j 10 a, b = + a i,j 00 , where a i,j 11 = a : C b dV ∆tn Ω i : Cε p εp Ω HσCε p i : Cε p for j = {1, . . . , m} do a i,j 10 = Ω Calculate j dV, ∀i, j ∈ {1, • • • , m} . εp i : H -1 σ εp j dV, Ω εp m+1 , εp j λj = λj + λm+1 εp j , εp j
	end
	Calculate
	εp m+1 =	εp m+1 -	m j=1	εp m+1 -m+1 , εp εp j εp j , εp j εp m j=1 j , εp εp m+1 , εp j εp j , εp j m+1 -j=1 εp j m	εp εp j , εp j m+1 , εp j	j εp	1/2
	and
	λm+1 = λm+1		εp m+1 -	m j=1	εp m+1 , εp j εp j , εp j	εp j , εp m+1 -	m j=1	εp εp j , εp j m+1 , εp j	j εp	1/2
	if
	[0,T ]

  établie à partir de la condition cinématique d'admissibilité, qui est définie telle que

	-
	Ω×[0,T ]

E.1) avec le tenseur des contraintes σ statiquement admissible, et δ u qui est un champ cinématiquement admissible. D'autre part, la compatibilité des déplacement prescrits avec la déformation générée dans la structure est

  la partie déviatorique du tenseur des contraintes. La limite d'élasticité est définie également par le modèle de Marquis-Chaboche :

  th la valeur limite, n d et k d des paramètres spécifiques au matériau considéré. L'initiation d'une fissure macroscopique est indiquée par l'endommagement critique Dc, valeur pour laquelle la rupture du matériau est considérée atteinte. E.3 Approche par réduction d'ordre de modèle basée sur la méthode LATIN pour les calculs d'endommagement L'approche LATIN cherche à résoudre les équations définissant le problème sur l'ensemble du domaine spatio-temporel à chaque itération. Les équations sont traitées itérativement, à savoir l'équilibre global de la structure d'une part, la loi élastique et les équations d'évolution non linéaires d'autre part.

	σy.	(E.5)
	Les paramètres relatifs aux matériaux k et n décrivent le comportement visqueux. La fonction limite
	relative à l'endommagement est donnée par
	f d = Y -Y th	(E.6)
	avec Y

  le problème global et linéaire est résolu dans l'espace A Il peut être remarqué que l'éq. (E.3d) bien que non linéaire est abordée avec le groupe des équations linéaires comme une étape de post-traitement à partir de la connaissance du tenseur des contraintes et de la variable d'endommagement à la fin de chaque itération. La solution exacte du problème sex est définie comme l'intersection des deux variétés par

	d qui appartient à la variété des conditions
	d'admissibilité éq. (E.1) et éq. (E.2), des lois d'état linéaires éq. (E.3b) et éq. (E.3c), et de la loi
	d'état non-linéaire pour l'endommagement éq. (E.3d) ;
	-le problème local et non-linéaire est résolu dans l'espace Γ qui appartient à la variété des équations
	d'évolution éq. (E.4) et de la loi d'état élastique éq. (E.3a) qui n'est pas linéarisable à cause de
	l'endommagement.

sex ∈ A d ∩ Γ.

(E.7)

  opérateur H -appartient à l'espace tangent associé à la solution ŝi+1/2 dans la variété Γ et C est le tenseur de Hooke pour le matériau non endommagé. En considérant que la variable d'endommagement n'est pas actualisée durant l'étape linéaire, l'opérateur de direction de recherche b - Z étant la part découplée de H -qui relie les variables internes à leurs variables associées. Les variables d'écrouissage sont ainsi obtenues en résolvant l'éq. (E.10a) en temps à chaque point de Gauss. Il peut être mentionné que la validité de l'éq. (E.11) présuppose que les variables d'écrouissage isotropique aient été transformées de telle sorte que l'équation d'état devienne linéaire, c'est-à-dire qu'elle suivent une ≪ formulation normale ≫ (voirBhattacharyya et al., 2017, Cognard and[START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF], pour plus de détails).La difficulté pour calculer les contraintes et déformations par rapport aux travaux précédents également basés sur la méthode LATIN est que loi d'état élastique éq. (E.3a) est ici potentiellement non-linéaire à cause de la possibilité d'endommagement. La présence d'un problème non-linéaire s'oppose a priori à l'utilisation d'une technique de réduction de modèle. L'idée proposée ici vise à transformer ce problème non-linéaire en deux équations indépendantes en décomposant les contraintes et les déformations totales en un partie dépendant de la déformation plastique d'une part, et une partie dépendant du niveau d'endommagement d'autre part.

	Les quantités d'intérêt à ce point σ i+1 , ε e i+1 et εp i+1 sont représentées en forme corrective à l'itération i + 1, c'est-à-dire
	est défini comme nul.
	La première étape étant le calcul des variables d'écrouissage, les équations d'état sont combinées sous
	la forme
	Z i+1 = ΠX i+1	(E.11)

où Π est un opérateur linéaire incluant les paramètres relatifs aux lois d'état. L'équation de direction de recherche pour les variables d'écrouissage éq. (E.10a) combinée avec l'équation d'état éq. (E.11) peut être écrite -( Ẋi+1 -Xi+1/2 ) = H Z (ΠX i+1 -Ẑi+1/2 ), (E.12) avec H

  est le nombre de paires considérées dans la décomposition, et C est un opérateur linéaire qui lie les fonctions spatiales relatives aux contraintes et aux déformations plastiques.E.3.2.2 Mise à jour de la base réduiteL'algorithme glouton est tel qu'après avoir défini une première paire de fonctions d'espace et de temps à la volée, à chaque itération une première décomposition est recherchée en réutilisant les fonctions spatiales prédéfinies et en actualisant les fonctions temporelles. Cette étape est équivalente à une décomposition orthogonale aux valeurs propres (POD) sur la base spatiale établie. En considérant que m modes spatiaux-temporels ont été générés pour approximer les contraintes et les déformations plastiques à l'itération i, les termes correctifs des contraintes et des taux de déformation plastique à l'itération i + 1 sont donnés par Les fonctions temporelles mises à jour sont calculées par minimisation d'un résidu mécanique qui est défini comme la norme de l'opérateur définissant la direction de recherche, c'est-à-dire Puis, si l'amélioration de l'approximnation n'est pas suffisamment importante, une nouvelle paire est ajoutée à la décomposition.E.3.2.3 Enrichissement de la base spatio-temporelleL'objectif de la phase d'enrichissement est d'ajouter une nouvelle paire de fonctions spatio-temporelles. Les corrections des contraintes et déformation plastiques sont écrites comme intention de calculer les quantités séparables λm+1 et εp m+1 . Une stratégie hybride est utilisée, la fonction spatiale εp m+1 est calculée à partir d'une formulation de Galerkin, en utilisant la condition d'admissibilité cinématique éq.E.2, tel que ∀δσ statiquement admissible à zéro

			m
	εp ( x, t) =	λj (t) εp j ( x) ,
			j=1 m						(E.17)
	σ ′ ( x, t) =	λ j (t) Cε p j ( x) ,
			j=1
	où m ∆ εp i+1 ( x, t) =	m	∆ λj (t) εp j ( x) ,
			j=1 m						(E.18)
	∆σ ′ i+1 ( x, t) =	∆λ j (t) Cε p j ( x) .
			j=1
		m				m
	{∆λ j } m j=1 = arg min {∆λ i } m j=1	j=1	∆ λj	εp j -Hσ	j=1	∆λ j Cε p j + ∆	H -1 σ	.	(E.19)
	∆ εp i+1 ( x, t) = λm+1 (t) εp m+1 ( x) , ∆σ ′ i+1 ( x, t) = λ m+1 (t) Cε p m+1 ( x) ,			(E.20)
	avec l'

Table E

 E Le problème est résolu pour U 0 = 0.070 mm en considérant d'une part l'approche mono-échelle et d'autre part l'approche bi-échelle avec les différents maillages temporels détaillés dans la table E.1. L'erreur relative est basée sur l'endommagement et est définie comme suit (Dms -Dts) • (Dms -Dts) dΩ dt Ω×[0,T ] Dms • Dms dΩ dt + Ω×[0,T ] Dts • Dts dΩ dt , (E.33) où Dms et Dts sont les endommagements obtenus en utilisant les algorithmes LATIN-PGD mono-échelle et bi-échelle respectivement. L'erreur relative qui est représentée dans la fig. E.7, croît grandement avec la taille des éléments. Le temps de calcul relatif est estimé par la formule t cal ts /t cal ms où t cal ms est le temps de calcul en utilisant l'algorithme mono-échelle et t cal ts le temps requis pour le calcul bi-échelle. Ce temps relatif, comme prévu, diminue lorsque le nombre de cycles par élément temporel augmente.

	er 2 d =	Ω×[0,T ]
		.1: Différentes discrétisations temporelles
	Discrétisation	Nombre de cycles	Nombre total
	temporelle	par élément temporel	d'éléments temporels
	I	50	20
	II	100	10
	III	200	5
	IV	500	2

Table E .

 E 2: Discrétisation temporelle adaptative basée sur le niveau d'endommagement

	Niveau d'endommagement	Nombre de cycles par élément temporel

Chapter 6-Conclusions and Future Perspective sophisticated damage models, for instance the 2-scale damage model can also be included, if the objective is to perform HCF simulations. Inertia effects, although excluded in the current work, can be introduced if the intent is to calculate dynamic fatigue.

As far as the two-scale discretisation is concerned, a better criterion for the adaptive temporal meshes may also be included in the future. Although, simple cyclic loading has been tested, the challenge of simulation of block loadings, or narrow-band random loading should also be explored.

Lastly, all these problems were tested in simple academic problems, and trials should be done if this concept can be utilised in practical engineering problems.

Appendix A Calculation of the finite element operators

The objective herein is provide a brief methodology for the construction of the FE operators E and C defined in section 4.6.1.

To derive the finite element operator E introduced in eq. (4.44), consider the spatial problem eq. (4.43), i.e. where Ktot and Ltot are the assembled matrices, {ūtot} and εp tot are the total number of nodal unknowns and Gauss point quantities respectively. It is needless to mention that although the stiffness matrix Ktot depends on the total number of nodal degrees of freedom, Ltot depends also on the total number of Gauss points. The right hand side of eq. (A.4) can be viewed as a "pseudo nodal force vector", hence care should be taken to include contribution of all the Gauss points to the "pseudo nodal force vector" at each element while building Ltot. Thereafter the procedure is quite straightforward. An operator

is calculated, where removal of the corresponding rows and columns of Ktot and the corresponding rows of Ltot where the nodal degrees of freedom are known, results in K red tot and L red tot respectively. The

Orthonormalisation of the space functions

The spatial basis functions obtained through the PGD-based technique as described in section 4.6.2.2 must be orthonormal to each other. The most common method to orthonormalise the spatial basis vectors is to use a Gram-Schmidt algorithm.

Considering εp 0 be the elastic initialisation, at any LATIN iteration i, m is the number of modes used to approximate the plastic strain rate, i.e.

Considering the case where a space-time mode is added to describe the plastic strain rate

where {ε p } m j=1 are orthonormal basis vectors but εp m+1 is not orthonormal to {ε p } m j=1 . The spatial vector εp m+1 is orthonormalised by algorithm 2, with respect to the existing basis. To maintain equality of the solution, the temporal functions are also modified. The number of modes m is updated depending on the contribution of the new but modified temporal function and the plastic strain rate can be written as

This process allows to compare the contribution of the different modes. All the modes that are generated are not considered, depending on the norm of the time functions corresponding to the orthonormal spatial basis. However, because of the Gram-Schmidt projection, useful information pertaining to the rejected modes is still taken into account.

103

Algorithm 2: Orthonormalisation of space functions Write

Alternative method of incorporating non-linear elastic state law

In chapter 4, the incorporation of damage in LATIN-PGD has been performed by including the nonlinear elastic state equation in the local stage. However it is also possible to incorporate the state law in the global stage.

In the local stage, the search direction equations are

Along with these, the evolution equations are solved.

In the global stage, the static admissibility condition can be written ∀δ u which is kinematically admissible to zero, as

C is the effective Hooke's tensor that depends on the damage variable. To separate the static admissibility condition into space and time, the effective Hooke's tensor needs to be represented in a separable form.

D.1 Series expansion of non-linear elastic law

At any LATIN iteration i + 1, the effective Hooke's tensor Ci+1 is expressed in the form of series of linear state laws. As per the definition of proportional loading

D.2 Separable representation of the quantities of interest

The corrective terms at LATIN iteration i + 1 are now written as

where q is the number of SVD modes used to describe the non-linear state equation as series of linear state laws. Then each sub-corrective term is written in terms of separable form of space and time.

The development of the PGD formulation starts from the static admissibility condition defined by eq. (D.2) and is similar to the development in section 4.6.1. The only thing to be noted is because of eq. (D.10), the Hooke tensor is also represented in a separable form.

For any sub-iteration l, using the same approximation as that of section section 4.6.1 gives

εp . The linear operators E l i+1 and C l i+1 are calculated the same way as described in Appendix A. The only difference is they will not be constant but change at each sub-iteration l. Subsequently, the space-time modes are calculated using the methods described in section 4.6.2. It is evident that the time functions for space and plastic strain are no longer equal. As far as the two spatial bases are concerned, the relation between them is not consistent, and orthonormalisation of one basis does not necessarily guarantee orthonormality of the other, hence both bases need separate orthonormalisation.

The limitation of this method lies on the fact that the Hooke's tensor needs to be separated, which is not possible for non-proportional loading. However, for the case of proportional loading this method shows appreciable accuracy compared to the method described in chapter 4. This method has been tested in the same one-dimensional problem (classic example of proportional loading) as given in section 4.8.1 for 5 cycles. For simplicity, the bar has been considered to be made of a single material with σy = 85 MPa. Considering the solution obtained using the method described in chapter 4 as reference, the solution obtained through the method involving separation of Hooke's tensor furnished an error of 1%, where the error is defined by eq. (5.19)). This particular method, although opens a new portal in solving damage mechanics using LATIN-PGD algorithm, is not advisable. The reason is not only the limitation in solving non-proportional loading, but also the higher numerical expense due to sub-iterations as compared to the method described in chapter 4.

Appendix E Extended summary in French E.1 Introduction

Le phénomène de fatigue est d'une grande importance dans le dimensionnement des pièces mécaniques ainsi que des structures de génie civil. La fatigue dans un sens général peut être définie comme le changement de propriétés d'une structure soumise à un chargement répété. Ces travaux s'intéressent à la fatigue mécanique, c'est-à-dire aux fluctuations créant des contraintes dans le maté-riau. La nature du chargement peut être parfaitement périodique, non-périodique ou aléatoire. Quelle que soit la nature de la fatigue, la rupture de la structure a lieu suite à trois phases (voir Suresh, 2001) :

-la nucléation et croissance de micro-vides qui amorce une fissure macroscopique, -la propagation stable de la fissure macroscopique, -la propagation instable de la fissure qui mène à la rupture totale du matériau.

Trois approches sont principalement utilisées par les ingénieurs pour garantir la sécurité d'une structure (voir Suresh, 2001). L'approche dite ≪ safe-life ≫ vise à assure que la structure demeure sûre durant un certain nombre de cycles de chargement. L'approche dite ≪ fail-safe ≫ requiert que la structure soit capable de supporter un certain niveau d'endommagement sans mener immédiatement à sa rupture totale, qui serait catastrophique. Le dimensionnement en ≪ tolérance à l'endommagement ≫ s'intéresse à la capacité d'une structure à survivre à la présence d'endommagement dans l'attente d'une réparation. La vie totale en fatigue est la somme du nombre de cycles correspondant à l'initialisation de la macro-fissure et du nombre de cycles correspondant à la propagation menant à la rupture catastrophique. Les philosophies de dimensionnement parmi les ingénieurs varient selon le stade de la rupture qui est étudié, mais également selon les propriétés du matériau, le chargement, la géométrie de la structure, etc.

Quel que soit l'objectif de l'ingénieur, de nombreux développements en recherche ont lieu pour représenter les phénomènes complexes mis en jeu durant la fatigue cyclique ou aléatoire avec des modèles flexibles et pertinents.

E.1.1 Différentes modélisations proposées pour l'analyse de la fatigue

Les approches historiques pour l'analyse de la fatigue sont basées sur l'exploitation de courbes empiriques introduites par Wöhler (voir Wöhler, 1860). La représentation du niveau de contraintes imposées à la structure en fonction du nombre de cycles durant lesquels la structure survit, observé durant les expériences, produit des courbes dites courbes S-N . Ceci peut être un outil efficace pour une analyse directe de type ≪ safe-life ≫ si un chargement parfaitement périodique et de même contrainte moyenne que celle de la courbe S-N est considéré. À partir du comportement asymptotique de la courbe, une limite d'endurance, niveau de chargement en dessous duquel le matériau n'atteindra théoriquement jamais la rupture, peut être définie. Pour des chargements conduisant à des déformations plastiques conséquentes Title: A model reduction approach in space and time for fatigue damage simulation

An innovative numerical scheme is proposed to predict the life time of mechanical components under fatigue loading, using classical continuum damage mechanics based models. The large time increment (LATIN) method, which is a non-incremental solution technique and builds the solution iteratively for the complete space-time domain, along with proper generalised decomposition (PGD) that separates the quantities of interest with respect to space and time, is used as a numerical framework. The first goal of the research is to introduce damage as an internal variable in the LATIN-PGD framework in a reasonable computational time.

The PGD-based model reduction is however not enough to solve problems involving large number of load cycles. The second part of the research is to overcome this bottleneck by using a multi-time scale approach, that takes into account the rapid evolution of the quantities of interest within a load cycle and their slow evolution along the load cycles.