Federica

Martin

Renaud Pierre

Keywords: Turbomachinery, DIC, Rupture, Triaxiality, Mesh dependency, Mass scaling

This thesis

addresses the amelioration of the burst prediction of turboengines rotating parts. The principal axis of this study was the identification of the material behavior up fracture under multi-axial loading and the identification of the local condition at failure. To accomplish this, an integrated digital image correlation approach was employed. The other thematic axis of this thesis was the objective and robust prediction of rupture using the damage mechanics in the framework of explicit simulations and the study of the influence, on the prediction, of various scaling techniques used to decrease the computation time. Two materials were studied, the titanium alloy TA6V and the nickel alloy Udimet 500. The samples were defined to guarantee similar loading conditions as those in the disks in service and to allow the use of field measurement techniques. A commercial code (ABAQUS) is used in an integrated approach to Digital Image Correlation (DIC). The principle of this method is to optimize the correlation of the test directly with the model parameters without passing by an intermediate procedure requiring the reconstruction of the deformation field beforehand. This technique, which can be seen as highly regularizing, grants a considerable robustness to the method, which allows obtaining information from degraded speckle Ces remerciements viennent conclure en peu plus que trois années de thèse passées au LMT Cachan et chez Safran Turbomeca à Bordes à cote de Pau. Je tiens tout d'abord à remercier les membres du jury, notamment Khemais Saanouni pour avoir accepté d'en être le président ainsi que Pierre-Olivier Bouchard et Jean-Charles Passieux pour avoir pris les rôles des rapporteurs. Un grand merci à Olivier Paulien-Camy, Hervé Chalons et Jean-Phillip Ousty de m'avoir accueillir à Bordes et de m'avoir guidé au cours de la thèse malgré la grande distance. Un grande merci aussi à Louis Roux pour la discussion et d'être mon homme à Pau. Je remercie Olivier Allix de m'avoir proposé cette thèse, qui m'a permis de travailler à l'interconnexion du monde numérique et expérimental. Il m'a tant laissé une grande liberté dans mes recherches que me guidé dans les sujets difficiles. Egalement merci à Cuong Ha-Minh d'être le co-encadrant de la thèse et de m'aider dans l'amélioration des mes présentations scientifiques. Je remercie également Francois Hild pour avoir encadré la partie expérimentale de cette thèse, même si originalement pas prévu. Il s'est toujours pris le temps de répondre à mes questions concernant la corrélation d'images.

Concernant mes rencontres au LMT, je dois d'abord citer les gents qui font fonctionner le LMT, comme Xavier qui ma donne un autre regard sous la conception mécanique, mais aussi Lydia qui est toujours là pour organiser tous. Egalement merci à Boubou, Rémi l'autre Xavier et Frisou, sans eux le LMT ne serait pas une place si unique dans le monde de recherche.

Un grand merci à ceux qui ont assisté à mes

Résumé : Cette thèse s'intéresse à l'amélioration de la prévision de l'éclatement des pièces tournantes des turbomachines. L'axe de recherche principal a été celui de l'identification du comportement du matériau jusqu'à rupture pour des sollicitations multi-axiales représentatives et l'identification des conditions de rupture locale elles mêmes. Pour cela une approche basée sur la corrélation d'image numérique intégrée a été suivie. L'autre axe abordé dans cette thèse a été celui de la prévision objective et robuste de la rupture par la mécanique de l'endommagement dans le cadre de simulations explicites utilisant des techniques de scaling pour diminuer les temps de calcul. Deux matériaux ont été étudiés, un alliage titane TA6V et un alliage nickel Udimet 500. Des éprouvettes ont été définies pour permettre d'avoir des sollicitations proches de celles en service et pour permettre d'utiliser des techniques de mesure de champs. Un code commercial (ABAQUS) est utilisé dans une approche intégrée de la corrélation des images numériques (CIN). Le principe de ces méthodes est d'optimiser la corrélation directement à partir des paramètres du modèle sans passer par une procédure intermédiaire demandant de reconstruire préalablement les champs de déformations. Cette technique, qui peut-être vu comme très régularisante, confère à la méthode une grande robustesse, ce qui permet d'obtenir des informations même en présence de mouchetis dégradés lors des phases ultimes de chargement. L'ensemble des paramètres constitutifs des modèles peut alors être identifié avec un seul essai hétérogène. Dans l'étude un accent particulier a été mis sur l'analyse des champs de triaxialité, dont l'importance sur les conditions de rupture est suspectée, en utilisant des éprouvettes fines et des éprouvettes épaisses. Des lois de comportement ont pu être identifiées pour des niveaux de déformations plastiques jusqu'à 3 fois supérieurs à ceux atteints dans des essais de traction uni-axiaux et ceci en présence d'adoucissement global. L'intégration étroite des essais et des simulations a également permis d'analyser l'état de déformation et contraintes des éléments proches de la surface de rupture juste avant celle-ci et ainsi de tester des critères de rupture potentiels avec un grand nombre de données. Pour la seconde partie de la thèse, plus exploratoire faute de temps, une étude préliminaire a porté sur la possibilité d'utiliser des modèles à taux d'endommagement limités combinés avec des techniques de scaling pour utiliser des simulations explicites dans le cadre de sollicitations centrifuges quasi-statiques. Les premiers résultats montrent qu'une accélération est possible avec des niveaux d'erreurs relativement bas sur les vitesses de rotation conduisant à la rupture et ceci dans des temps comparables avec des simulations implicites quasi-statique, l'avantage étant la robustesse des analyses en dynamique explicite. patterns encountered at the end of the test. The entire set of constitutive parameters can be identified with only one heterogeneous test. Within this study, a particular accent was laid on the analysis of triaxiality fields by using thin and thick samples. The importance of these fields for the fracture conditions is suspected. Constitutive models were identified for plastic strains that were about three times higher than the ones achievable in uniaxial tensile test on smooth sample. The close integration of the experiments and simulations allowed the analysis of deformation and stresses of the elements near the failure surface at the instant just before failure. This allowed the testing of different criteria with a large amount of data.

For the second axis of the thesis, preliminary studies examined the possibility to use limited damage rate models combined with scaling techniques. These latter allow performing explicit simulations in context of a centrifugal quasi-static loading. The results show that acceleration is possible. The errors of the rotational speeds, which lead to disk fracture, are relatively low. The simulation times are comparable to implicit quasi-static simulations, while the main advantage is the maintained robustness in explicit simulations.

List of Figures

CHAPTER 1

Introduction

In its most general definition, a turbomachine transfers energy between a fluid and a rotor. Different classes of turbomachines exist: compressors, pumps or turboengines. Turboengines or gas turbines are rotating internal combustion engines, with a continuous combustion process. Depending on the application, several designs of turboengines exist. Among them, turboprops and turboshafts make use of an external rotor. The schematic design of a helicopter turboshaft engine is presented in figure 1.1. [START_REF] Mazière | Eclatement des disques de Turbomachines -Burst of turboengine disks[END_REF]: the combustion gas flows out of the combustion chamber in the high pressure turbine, which is part of the first shaft and drives the compressor. The power turbine stage of the second shaft transforms the remaining stream energy into mechanical energy. This shaft also ensures the connection to the gearbox.

Today's research in turboengines is driven by three main factors: fuel-efficiency, weightgains and cost-reduction. Therefore, manufacturers constantly seek for design improvements. These circumstances lead to more and more realistic computations, which require a better understanding and characterization of the constitutive behavior of the materials.

In this context, this thesis focuses on possible improvements in the methods of material modeling and characterization. In addition, initial work regarding possible improvements of the burst computation was conducted.

The design of turbomachines for aeronautical applications has to comply with strict regulation rules. The European Certification Specification for Engines (CS-E) requires the testing of every engine component for a variety of loading conditions. For rotating disks and impellers, the main tests concern low cycle fatigue life, disk integrity, modal analysis and crack tolerance. While every test adds different requirements to the design, the integrity assessment is maybe the most challenging one. The test assumes a hypothetical failure event: the gearbox connected to the power turbine may break or a combustion malfunction may occur. In this case, the rotational speed abruptly increases, exceeding the maximum design speed and entering what is called the overspeed regime. The CS-E demands electronic overspeed protection, which initiates the shutdown of the turbine. In addition, the certification requires a proof of the disk stability in the overspeed regime. To efficiently design the disk, the knowledge of the ultimate stability speed, the burst speed, is required.

The estimation of the burst speed can rely on global criteria. The average hoop stress criterion proposed by Robinson is one of the most used [START_REF] Robinson | Bursting tests of steam-turbine disk wheels[END_REF]). In the case of complex geometries and refined material models, this criterion is not accurate enough.

Recent attempts rely instead on local burst predictions. Finite element simulations were performed, which precisely take into account both the complex geometry of actual disks and the constitutive model of the material. Many studies have been carried out on this subject, notably the one of Mazière, in which many aspects were precisely studied, e.g. the effect of viscoplasticity or possible Portevin-Le Chatelier instabilities. From the work of Mazière, it also appears that burst prediction through local approaches strongly depends on the modeling of the constitutive behavior [START_REF] Mazière | Overspeed burst of elastoviscoplastic rotating disks -Part I: Analytical and numerical stability analyses[END_REF]. For the studied material, the Hosford yield criterion [START_REF] Hosford | A generalized isotropic yield criterion[END_REF], was found to be more appropriate than the von Mises criterion.

However, no clear improvement of the burst prediction was obtained. Therefore, we seek for other possible explanations. A first one is that the modeling and the identification of the used material have mainly been performed based on uniaxial tests and with a von Mises elasto-plastic model. The stress state in these tests is characterized by a triaxiality of about 0.33. In contrast, in the disk computations the triaxiality ranges from 0.4 to 0.8. A second possible explanation consists in the large values of equivalent plastic strain, which are reached in the burst calculations. The values lie in the post-peak regime of the typical strain-stress curves. In this regime, strain and stress heterogeneities occur in the tensile characterization test. This makes a classical identification, which assumes test homogeneity, inappropriate.

The last issue concerns the condition of failure itself. For the materials used, no clear information is available either about the possible occurrence of damage or the appropriate form of the local criterion.

From these considerations, we decided that progress should be possible on all these questions by applying full-field measurement techniques to heterogeneous tensile experiments. We primarily focused on stress triaxiality, as it is a key parameter in damage and failure [START_REF] Mcclintock | A Criterion for Ductile Fracture by the Growth of Holes[END_REF][START_REF] Gurson | Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I-Yield Criteria and Flow Rules for Porous Ductile Media[END_REF][START_REF] Lemaitre | A Continuous Damage Mechanics Model for Ductile Fracture[END_REF][START_REF] Malcher | An assessment of isotropic constitutive models for ductile fracture under high and low stress triaxiality[END_REF]. We also tried to investigate the influence of the choice of the plasticity criterion, which has been the subject of many studies in recent years (Bao et al., 2004a;[START_REF] Wierzbicki | Calibration and evaluation of seven fracture models[END_REF][START_REF] Malcher | An improved damage evolution law based on continuum damage mechanics and its dependence on both stress triaxiality and the third invariant[END_REF]. The use of flat notched tensile geometries gives rise to levels of stress triaxiality, which are more representative to the ones encountered in the disks. The fact that the studied titanium alloy exhibits an important tension-compression asymmetry was disregarded [START_REF] Tuninetti | Anisotropy and tension-compression asymmetry modeling of the room temperature plastic response of Ti-6Al-4V[END_REF]. In fact, for the applications of interest, compression does not occur.

The considered sample geometries induce non-uniform strain fields. In this case, the measurement of load-dependent displacement fields and their subsequent exploitation are necessary [START_REF] Wattrisse | Analysis of strain localisation during tensile test by digital image correlation[END_REF]. Different optical techniques can be used to measure displacement fields (Rastogi et al., 2012). Among them, Digital Image Correlation (DIC) is increasingly employed due to its versatility and applicability to "any" scale of observation (Rastogi et al., 2012).

Even though very powerful, DIC suffers from the 'resolution/spatial resolution' curse since it is an inverse problem (Hild et al., 2012a). In particular, it does not allow elements of arbitrarily small sizes to be used. For the analyzed tests, this limitation occurs for a minimum element size of 60 pixels, due to the poor random texture. To correct for this shortcoming, it was decided to rely on Integrated DIC (I-DIC) approaches, which allow linking even more closely measurements and simulations. This technique was first applied with closed-form (i.e. elastic) solutions [START_REF] Hild | Digital image correlation: From measurement to identification of elastic properties -A review[END_REF] and then generalized to numerically generated solutions [START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF].

The use of these techniques allowed for the qualification of the constitutive models for levels of equivalent plastic strain of up to 40%. These levels are far beyond the beginning of the post-peak regime in tensile test on un-notched specimens, which is about 10 % for the studied titanium alloy. Moreover, these levels of equivalent plastic strain were obtained for a range of stress triaxiality ranging from 0.4 to 0.7.

Similarly to existing studies, the presented analysis of failure relies on the identified elasto-plastic model [START_REF] Diiorio | Plastic deformation, damage and rupture of PM Ti-6Al-4V at 20K under monotonic loading[END_REF][START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF][START_REF] Papasidero | Ductile fracture of aluminum 2024-T351 under proportional and non-proportional multi-axial loading: Bao-Wierzbicki results revisited[END_REF]. In addition, we propose an image-based analysis of the state of stresses and strains in the vicinity of the failure surface at the instant just before failure. This analysis allowed testing different criteria with a large amount of data.

The experimental part of this thesis and the associated software developments did not allow us to go as far as we expected regarding the last aspect, namely, the simulation of failure itself. This part, which was conducted in parallel, is at present relatively decoupled from the experimental finding.

Today, realistic simulations of complex industrial problems including all their technological complexity can only be made using commercial codes. In the case of the simulation of failure, two main difficulties arise. The first one concerns only quasi-static or slow dynamic problems. In this case, implicit schemes are the dominant approaches. The numerous instabilities associated with damage and cracking make the convergence problematic and sometimes impossible. Here, path following techniques are mandatory and are often used in combination with viscosity. In order to ensure convergence and robustness, explicit algorithms can be used. In this case, mass scaling techniques allow for the increase of the critical time step. It is unclear though how the use of these techniques affects the prediction of failure.

The second difficulty is the occurrence of spurious mesh dependence in the failure prediction. The most widespread technique to obtain mesh independent results is the one of non-local models [START_REF] Besson | Continuum models of ductile fracture: a review[END_REF]. The implementation of such models in industrial software is more than cumbersome and the identification of these models remains an issue.

A much easier possibility is the development of a bounded rate version of the model [START_REF] Allix | The bounded rate concept: A framework to deal with objective failure predictions in dynamic within a local constitutive model[END_REF]. This approach requires only classical, i.e. local, modifications of the constitutive law. Nevertheless, such approaches have been much less studied and in any case have to be performed in combination with dynamic analyses. Considering these two issues, the proposed path is the combination of explicit simulations and a bounded rate model. This thesis explores the first steps of such an approach.

Outline

This document is structured in six chapters. A literature review on the topics concerned follows this introduction. The identification of the plastic behavior constitutes the main part of this thesis. Two separate chapters address possible solutions and perspectives in the identification of fracture models and the successful simulation of quasi-static mesh independent fracture. A conclusion ends the investigation.

Chapter II: State of the art This chapter introduces the topics of burst prediction, elasto-plastic modeling and the digital image correlation methods. The chapter begins with the description of the studied titanium and nickel alloys as well as the experimental phenomena, which occur during deformation. The history of disk burst prediction and the encountered issues are presented. Plasticity and rupture modeling are presented in light of the stress state influence. Finally, DIC and I-DIC methods are explained.

Chapter III: Experimental study of plasticity using DIC and I-DIC The chapter presents the conducted experiments. The choice of the experimental geometries for high triaxialities is explained and the procedures that are necessary to obtain the model parameters are detailed. Different aspects of the material model identification via I-DIC are addressed. Finally, different models are discussed for the studied titanium and nickel alloys.

Chapter IV: Towards the identification of material failure criteria The identification of damage and failure induces additional challenges to the identification procedure. This chapter discusses the experimental results. A characterization of the failure modes of the two materials is presented. The analysis of the fracture surfaces in connection with the I-DIC approach allows the definition of different fracture loci.

Chapter IV: Perspectives for the simulation of material failure The introduction of damage and failure leads to a mesh dependence and a possible nonconvergence of implicit simulations. The influence of the centrifugal forces on the mesh dependence is explored. Explicit dynamic simulations are applied to simple cases and the resulting errors are analyzed. In this context, a literature review on the current state of mass and time scaling techniques is provided.

CHAPTER 2

State of the art

This chapter focuses on establishing the state of the art of the main aspects studied in this thesis, namely, the disk burst prediction, the modeling of metals including failure especially with respect to multi-axial loading, the identification of such models and in particular the use of Digital Image Correlation (DIC).

These aspects are as much as possible oriented towards the case of TA6V and Udimet 500 alloys and are put in relation with the context of the thesis, namely, the prediction of the burst speed limit in turbomachines.

Description of TA6V and Udimet 500 alloys

In turboengines, some components undergo large rotational speeds. Moreover, the internal combustion leads to severe gradients in both temperature and mechanical fields. These operating conditions require the use of materials with an ultimate strength larger than the one of standard engineering steels. High creep resistance, for lifetime safety, and preferable low weight, for aeronautical applications, are further requirements. The compression of the intake air and the subsequent combustion raise the engine temperature to a minimum of 300 ∘ C. The air pollution at airports leads to the phenomenon of hot corrosion as the particles enter the compressor. This requires a good corrosion resistance, not only in the turbine section, but also in the compressor [START_REF] Ejaz | Failure due to structural degradation in turbine blades[END_REF].

Due to the strength requirements, the disks have a bulky structure. Forging provides a cost-effective production process. The process-induced anisotropy is removed through intensive successive heat-treatments. This allows the assumption of a high level of isotropy in the disks.

As this can be visualized in figure 2.1, lightweight titanium-based alloys are among the best candidates to meet these requirements. Due to the reduced resistance at high temperatures, titanium alloy applications are restricted to compressor parts. There, the temperatures rarely exceed 600 ∘ C [START_REF] Rao | Materials for Gas Turbines -An Overview[END_REF]. In the turbine section the combustion gas can exceed 700 ∘ C. Nickel-based alloys combine high strength with a significant higher temperature tolerance. Although they are almost twice as heavy as titanium alloys, they are frequently used in the turbine section.

Due to their application in helicopter engines, this study examines the titanium alloy Ti6Al4V or TA6V and the nickel-based super alloy NiCr19Co18Mo4Ti3Al3 or Udimet 500 [START_REF] Lovatt | Material selection and processing[END_REF].

(U500). Their actual microstructure depends strongly on the heat-treatment. Examples of the microstructures of the two alloys are displayed in figure 2.2. The alloy composition and the final morphology of the grains determine the mechanical properties.

The titanium-based alloy TA6V

The titanium alloy Ti6Al4V is a standard aeronautical alloy with widespread use in the structural parts of planes and engines. TA6V alloy combines several titanium phases. The two main phases are a hexagonal phase, 𝛼, and a body centered cubic phase, 𝛽, cf. figure 2.2(a). Additional elements provide stability and corrosion-resistance to these phases. Table 2.1 presents an overview of the alloy elemental composition. The approximate density is 𝜌 = 4.42 𝑔/𝑐𝑚 3 , close to the density of pure titanium 𝜌 = 4.51 𝑔/𝑐𝑚 3 . Depending on the heat treatment, the typical mechanical properties are characterized by a Young's modulus ranging from 𝐸 = 110 to 140 𝐺𝑃 𝑎, a Poisson's ratio of 𝜈 = 0.32, an ultimate stress 𝜎 𝑚𝑎𝑥 ranging from 900 to 1200 𝑀 𝑃 𝑎 and a final "macroscopic" elongation measured on tensile test 𝜀 𝑚𝑎𝑥 ranging from 13 to 16 % [START_REF] Abouridouane | Bruchverhalten von Leichtmetallen unter Impact-Beanspruchung[END_REF].

The nickel-based superalloy Udimet 500

The increased temperature and corrosion resistance of the nickel-based alloy Udimet 500 is associated with the higher amount of chromium and cobalt. In a typical composition, these sum up to 35 % of the total weight, cf. table 2.1. The high amount of chromium has a significant impact on the density, 𝜌 = 8.08 𝑔/𝑐𝑚 3 , which is lighter than pure nickel, 𝜌 = 8.9 𝑔/𝑐𝑚 3 . The material consists of a primary body centered cubic phase, 𝛾, and a strengthening face centered cubic phase, 𝛾 ′ . The large quantity of additional alloying elements and higher amounts of carbon lead to a significant formation of carbides [START_REF] Davis | Nickel, cobalt, and their alloys[END_REF], cf. figure 2.2(b). The high content of alloying addition increases the probability of large inclusions. Nevertheless, the carbides are well dispersed and agglomerate mainly at the grain boundaries [START_REF] Ejaz | Failure due to structural degradation in turbine blades[END_REF]. Similar to TA6V the large variety in heat-treatments leads to a large range of mechanical properties: an elastic behavior with a Young's modulus around 𝐸 = 200 𝐺𝑃 𝑎 and a Poisson's ratio of 𝜈 = 0.3. The ultimate stress 𝜎 𝑚𝑎𝑥 ranges from 690 to 1300 𝑀 𝑃 𝑎 and a final "macroscopic elongation" 𝜀 𝑚𝑎𝑥 from 4 to 18 % [START_REF] Lambert | Structural stability of Udimet 500 a nickel base superalloy[END_REF].

(a) TA6V microstructure: (1) 𝛼-phase (dark structures) is filled with (2) 𝛽-phase (white patches), after [START_REF] Breutinger | Verformungsverhalten und Verformungskinetik von Titan technischer Reinheit und der Titanlegierung TiAl6V4 im Bereich niedriger homologer Temperaturen von 0[END_REF].

(b) U500 microstructure: (1) 𝛾-phase (light gray) is filled with (2) 𝛾 ′ -phase (black), (3) carbides exist (dark gray), after [START_REF] Sajjadi | Effects of cooling rate on the microstructure and mechanical properties of the Ni-base superalloy UDIMET 500[END_REF].

Microstructural processes of plastic deformation

In the process of metallic deformation, one distinguishes three phases: elastic, elasto-plastic and material failure. During elasticity, the reversible extension of the atomic bonds allows the storage of external energy. For most metals, the amount of energy stored in elastic deformation is small. Therefore, it remains macroscopically linear [START_REF] Chaboche | Mécanique des matériaux solides[END_REF]. Due to the resistance of the atomic bonds, the internal stresses increase. The high macroscopic stiffness is a result of the large increase in stress for small strains. As the external energy grows above a certain threshold, the yield stress 𝜎 𝑦 , a second mechanism of energy storage becomes active, namely, the reorganization of the crystal lattice through the nucleation of dislocations and dislocation movement [START_REF] Besson | Mécanique non linéaire des matériaux[END_REF]. Dislocations are line defects in the crystal lattice. These processes allow the absorption of large amounts of energy. The dislocation motion and creation is irreversible, and so is the macroscopic phenomenon of plasticity. This enables for the possibility of large permanent strains.

Other phenomena characterize plasticity. The creation of new dislocations hinders the movement of other dislocations. An additional stress increase is necessary to overcome this locking. Macroscopically, a further increase in internal stresses is necessary for an increase in deformation; the material hardens. This effect diminishes with increasing deformation; the hardening saturates, cf. figure 2.3. Dislocation movement happens without volume change.

A more complex phenomenon is observed for cyclic behavior associated to difference with the yield stress in tension and compression [START_REF] Bauschinger | Über die Veränderungen der Elasticitätsgrenze von Eisen und Stahl[END_REF]. This Bauschinger effect does not occur under purely tensile load states encountered during disk burst. Failure of the material can occur in different manners, brittle or ductile, which is further discussed in section 2.5. Figure 2.3 shows a schematic representation of typical stress-strain for a ductile and a brittle material.

Burst speed calculation methods for turbomachine disk

Background of burst speed predictions

The design requirements of all components do not permit for irreversible plastic deformation even at the highest rotational speed, which occurs in normal operation mode. As explained earlier, it is possible that the engine accelerates to speeds greater than this highest design speed, entering the overspeed regime. To ensure flight safety, Airworthiness Requirements (CS-E ) demand the proof of the structural integrity at these higher speeds. In practice, the integrity speed is 120 % of the maximum design rotational speed. To that end, turboengine manufacturers have to conduct disk integrity analyses and disk integrity tests. Structural failure of the disk occurs after an additional speed margin. This final speed is referred to as the burst or limit speed.

Due to high costs of disk burst experiments, the design process relies considerably on the predictions of the burst speed. Predictions of lower precision have to be balanced with higher design security factors. Consequently, the precision of the burst speed predictions is a key aspect in the design process.

The research conducted on disks bursts focuses on three main areas, namely, the correct prediction of the limit speed, the prediction of the number of burst fragments and the prediction of fatigue cracking [START_REF] Mazière | Eclatement des disques de Turbomachines -Burst of turboengine disks[END_REF]. The current study is related to the first point.

Experimental study of burst speed

The experimental overspeed test takes place in a closed pit with massive walls, which allows scarce possibility for observation and measurements. Therefore, most of the experimental findings concentrate on a post-mortem examination of the disk fragments. Based on the fracture surface [START_REF] Servetnik | Energy-based method for gas turbine engine disk burst speed calculation[END_REF], disk fracture may be classified into two basic types. In the hoop mode burst, the entire disk disintegrates along radial planes in two or more parts (cases one and two in figure 2.4). The second type is the rim peel burst. A portion of the outer diameter breaks away from the intact central hub section. Frequently, the burst fragments present a combination of these two types (cases three and four in figure 2.4). The post-failure analysis remains complicated, due to secondary fracture as the disk fragments impact against the pit wall.

Figure 2.4: Disk fragments after burst: the disk fracture can happen along radial planes or be more complex [START_REF] Nozhnitsky | Numerical Simulation of Spin Testing for Turbo Machine Disks Using Energy-Based Fracture Criteria[END_REF][START_REF] Kuzmin | Yield Surface Investigation of Alloys During Model Disk Spin Tests[END_REF][START_REF] Forest | Mécanique des milieux continus[END_REF].

While the experimental burst speed is externally recorded, the experimental deformation of the disk remains mostly unknown. A counter-example was reported [START_REF] Kuzmin | Yield Surface Investigation of Alloys During Model Disk Spin Tests[END_REF]. The rim displacement was measured with the eddy current technique. In most cases, the comparisons between experiments and simulations can only be made on the limit speed. This requires a proper definition of a burst criterion and its application in the simulations.

Prediction of the burst speed with global criteria

The calculation of the disk burst speed has attracted the attention of researchers since a long time. Amongst the first studies are the solutions of elastic strain fields in simple disk geometries [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF][START_REF] Timoshenko | Theory of elasticity[END_REF][START_REF] Roark | Formulas for stress and strain[END_REF]. These analytical solutions allowed the definition of the main criterion for the estimation of the burst speed. Burst occurs when the mean hoop stress reaches the maximum tensile strength of the material [START_REF] Robinson | Bursting tests of steam-turbine disk wheels[END_REF][START_REF] Laszlo | Rotating disks in the region of permanent deformation[END_REF][START_REF] Waldren | Burst Strength of Rotating Discs[END_REF][START_REF] Percy | An experimental study of the burst strength of rotating disks[END_REF]. This criterion is considered as appropriate for the burst prediction of rotating rings, but not for large bulk disk [START_REF] Tvergaard | On the burst strength and necking behaviour of rotating disks[END_REF]. For non-axisymmetric disks, studies were carried out to demonstrate its limits [START_REF] Mazière | Eclatement des disques de Turbomachines -Burst of turboengine disks[END_REF]. A correction factor was proposed for the application of this criterion to disks [START_REF] Manavi | Centrifugal Rotor Tri-Hub Burst For Containment System Validation[END_REF].

Recent developments emphasize two issues: the influence of the geometric non-linearity and alternative burst criteria. In fact, even though geometrical effects could be thought of as negligible they play an important role in the mechanism of the instability as is explained in what follows. Let us consider, in a Finite Element (FE) framework, the expression of the centrifugal forces 𝐹 𝑐𝑒𝑛𝑡 ,

𝐹 𝑐𝑒𝑛𝑡 = 𝑀 (𝑥 0 + 𝑢) θ2 , (2.1)
where 𝑀 is the mass matrix and θ denotes the rotational speed. Centrifugal forces evolve in the disk depending on the initial position 𝑥 0 but also on the displacement 𝑢 (geometrical effect). Therefore, the load also increases even if the angular velocity remains constant. This effect is known as spin-softening or centrifugal softening. The equilibrium equation reads

𝐹 𝑖𝑛𝑡 = 𝐾 𝐸𝑃 .𝑢 = 𝐹 𝑐𝑒𝑛𝑡 , (2.2)
where 𝐾 𝐸𝑃 is the elasto-plastic secant stiffness matrix. This equation can be rewritten as

(𝐾 𝐸𝑃 -𝑀 θ2 )𝑢 = 𝐾 𝑡 .𝑢 = 𝑀 .𝑥 0 θ2 . (2.3)
An increase in the rotational speed leads to a decrease in the total stiffness 𝐾 𝑡 . This analysis can be conducted further in order to derive instability criteria. This was done in simple cases [START_REF] Hong | A New Method for Calculating Burst Speed of Aeroengine Disks[END_REF][START_REF] Baddour | A revisit of spinning disk models. Part I: derivation of equations of motion[END_REF][START_REF] Lila | Stability Loss of Rotating Elastoplastic Discs of the Specific Form[END_REF] and in a more systematic manner [START_REF] Mazière | Overspeed burst of elastoviscoplastic rotating disks -Part I: Analytical and numerical stability analyses[END_REF].

In [START_REF] Mazière | Eclatement des disques de Turbomachines -Burst of turboengine disks[END_REF], FE elasto-plastic computations were performed taking into account spin-softening by means of large displacement and large strains. To illustrate such computations, figure 2.5(a) presents an axisymmetric disk model with a linear plastic hardening and a constant rotational acceleration. A global indicator for the total amount of deformation is the radial rim displacement, 𝑢 𝑟 , normalized with the initial position 𝑟 0 . In what follows, this external radial deformation is denoted by 𝜀 𝑟 such that

𝜀 𝑟 = 𝑢 𝑟 𝑟 0 .
(2.4)

Figure 2.5(b) shows the evolution of this global deformation and of the equivalent plastic strain in an external element denoted by 𝑁 . As long as no plasticity takes place in this element, and only a few areas of the disk have reached the yield stress, cf. figure 2.5(c); the external deformation 𝜀 𝑟 remains low. The instant the element enters in the plastic phase is marked with a dashed line. From this instant a large amount of plastic strain exists in several areas of the disk, cf. figure 2.5(d). The effect of the spin softening becomes then significant and the external radial deformation increases very fast, which characterizes the surge of instability, cf. figure 2.5(b).

Prediction with local approaches

The failure prediction of the disk may also be based on local criteria involving either the local strain energy density, the equivalent plastic strain, the principal stress, the principal strain or more sophisticated quantities [START_REF] Mazière | Eclatement des disques de Turbomachines -Burst of turboengine disks[END_REF][START_REF] Nozhnitsky | Numerical Simulation of Spin Testing for Turbo Machine Disks Using Energy-Based Fracture Criteria[END_REF][START_REF] Servetnik | Energy-based method for gas turbine engine disk burst speed calculation[END_REF]. Those criteria can be simply applied when post-processing the computation. A second option is the enrichment of the material behavior with a damage model affecting the constitutive law. This leads to much more complex issues, further discussed in chapter 5.

Apart from the definition of the criterion itself, its application raises the question of the sensitivity of the result with respect to the constitutive law. It seems that this aspect is much more crucial for local criteria than for global ones, especially for complex structures. For this reason, a parameter study of a plasticity model was performed for a disk [START_REF] Squarcella | The importance of the material properties on the burst speed of turbine disks for aeronautical applications[END_REF].

Influence of the equivalent stress on burst predictions Figure 2.6 shows the calculations for two disks. They were performed with von Mises yield criterion and Hosford yield criteria.

In the first case, only the experimental speed was recorded. Calculations were carried out until non-convergence. The prediction with the von Mises criterion overestimated the maximum experimental speed. In contrast, using the Hosford criterion led to a correct prediction. The optimum exponent for the criterion was 𝑛 = 80, which corresponds roughly to the Tresca yield criterion [START_REF] Mazière | Eclatement des disques de Turbomachines -Burst of turboengine disks[END_REF].

In the second case, a measurement of the experimental deformation was performed via eddy current technique. The level of deformation is lower, due to the different shape of the disks. The calculations were only carried out up to the maximum experimental speed. The prediction with the von Mises criterion is closer to the experimental curve up to the experimental burst speed. The authors report an optimum exponent of 𝑛 = 6 [START_REF] Kuzmin | Yield Surface Investigation of Alloys During Model Disk Spin Tests[END_REF]. However, it must be noted that the computation could have been prolongated after the experimental burst speed, which raises once again the question of the burst criterion.

Local criteria for burst predictions

A local criterion based on burst experiments was proposed in [START_REF] Nozhnitsky | Numerical Simulation of Spin Testing for Turbo Machine Disks Using Energy-Based Fracture Criteria[END_REF]. It introduces a pseudo equivalent plastic strain as follows,

𝜀 𝑝 𝑚𝑎𝑥 = 𝑝 𝑚𝑎𝑥 ( 1 - 𝜎 𝑚 𝜎 𝑓 ) 𝛼 .
(2.5)

The strain at fracture 𝑝 𝑚𝑎𝑥 and the stress at fracture 𝜎 𝑓 are the one, measured from uniaxial experiments which for the strain is always debatable. As in many criteria, the mean stress 𝜎 𝑚 affects the failure, while the parameter 𝛼 describes the embrittlement at high mean stress. For a nickel alloy, it is in the range of 𝛼 = [1..2]. These results indicate that the stress state can have an effect on the prediction. Let us note another aspect of the failure prediction, the one of defects. Even though a low-defect tolerance is in operation in the disk production, defects may still exist and possible effects were studied [START_REF] Scolavino | Critical Speed of Flawed Rotors: Global vs. Local Approach[END_REF].

Constitutive modeling of elasto-plasticity

The research fields that deal with the complete modeling of elasto-plasticity at a low scale are dislocation dynamics and crystal plasticity. With today's processing power, the application of these models to large-scale applications is not feasible in acceptable computation times. Instead, the phenomena are described in the theory of continuum mechanics using a homogenization approach. This is reasonable, considering the relation of their scale, 10 -7 𝑚𝑚, to the engineering application in tensile tests or large parts, 10 -1 -10 3 𝑚𝑚 [START_REF] Chaboche | Mécanique des matériaux solides[END_REF].

The field of plasticity modeling comprises many aspects. The following paragraphs are limited to topics of interest for this study. Notably, plasticity is considered associative and rate-independent. Moreover, we do not consider kinematic hardening, as this aspect has to be considered principally for cyclic loading.

For more details on these topics see [START_REF] Malvern | Introduction to the mechanics of a continuous medium[END_REF][START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF][START_REF] Besson | Mécanique non linéaire des matériaux[END_REF].

Description of the stress state Principal stresses and stress invariants

The Cauchy stress tensor 𝜎 can be characterized by three invariants the high, the intermediate and the low principal stress, 𝜎 1 ≥ 𝜎 2 ≥ 𝜎 3 . The maximum shear stress can then be expressed as,

𝜏 𝑚𝑎𝑥 = 1 2 |𝜎 1 -𝜎 3 |. (2.6)
For any symmetric tensor, the first invariant is its trace,

𝐼 1 = Tr { 𝜎 } = 𝜎 11 + 𝜎 22 + 𝜎 33 = 𝜎 1 + 𝜎 2 + 𝜎 3 , (2.7)
where 𝜎 11 , 𝜎 22 , 𝜎 33 are the normal stresses. The mean normal stress is 𝜎 𝑚 = 1 3 𝐼 1 . The stress tensor 𝜎 can be decomposed into the sum of two other tensors. The shear-free hydrostatic stress-tensor, 𝜎 ℎ , and the stress deviator, 𝑠,

𝜎 = -𝜎 ℎ + 𝑠 = -𝜎 ℎ 𝐼 + 𝑠 = 1 3 𝐼 1 𝐼 + 𝑠. (2.8)
For an isotropic material, the hydrostatic stress, 𝜎 ℎ = -𝜎 𝑚 , only produces volume changes.

In contrast, the deviatoric stress tensor is associated with isochoric deformation.

The tensor 𝑠 is analogously characterized by its invariants 𝐽 1 , 𝐽 2 , 𝐽 3 .

𝐽 1 = 0, 𝐽 2 = 1 2 (𝑠 2 1 + 𝑠 2 2 + 𝑠 2 3 ) = 1 6 [ (𝜎 11 -𝜎 22 ) 2 + (𝜎 22 -𝜎 33 ) 2 + (𝜎 33 -𝜎 11 ) 2 ] + 𝜎 2 12 + 𝜎 2 23 + 𝜎 2 31 = 1 6 [ (𝜎 1 -𝜎 2 ) 2 + (𝜎 2 -𝜎 3 ) 2 + (𝜎 3 -𝜎 1 ) 2 ] , 𝐽 3 = det(𝑠).
(2.9)

Definition of triaxiality and Lode angle

The triaxiality ratio 𝜂 is defined from the second invariant 𝐽 2 and the hydrostatic stress as

𝜂 = -𝜎 ℎ √ 3𝐽 2 = 𝜎 𝑚 𝜎 𝑒𝑞 , (2.10)
where 𝜎 𝑒𝑞 = √ 3𝐽 2 is the von Mises equivalent yield stress. Negative triaxiality values indicate compressive loading, positive values tensional loading. Other triaxiality indicators have been defined for example in [START_REF] Schafer | Triaxiality and Fracture of Steel Moment Connections[END_REF] as

𝜂 2 = 𝜎 1 𝜎 𝑒𝑞 . (2.11)
In recent years, the effect of the third invariant on fracture has been the subject of many studies. A convenient way of application is the use of the Lode angle that introduces the normalized third invariant 𝜉 [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF],

𝜉 = ( ( 27 2 𝐽 3 ) 1/3 √ 3𝐽 2 ) 3 = 𝑐𝑜𝑠(3𝜃), θ = 1 - 6𝜃 𝜋 = 1 - 2 𝜋 𝑎𝑟𝑐𝑐𝑜𝑠𝜉.
(2.12)

From the definition of θ it is recognized that θ = -1 corresponds to an axial symmetry loading in compression, θ = 0 to a plane strain state or a generalized shear one and θ = 1 represents an axial symmetry loading in tension [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF].

Isotropic hardening laws

The calculation of stresses results from 𝜎 = C 𝑡𝑎𝑛 𝜀 𝑒 using the tangent stiffness tensor, C 𝑡𝑎𝑛 .

The elastic strain 𝜀 𝑒 can be obtained from the total strain 𝜀 𝑡 via 𝜀 𝑒 = 𝜀 𝑡 -𝜀 𝑝 . In this case, the aim of every plasticity model is the description of the plastic strain 𝜀 𝑝 . A basic plasticity model for associative plastic flow consists of the normality rule, the Kuhn-Tucker conditions and a yield function 𝑓 ,

d𝜀 𝑝 = d𝑝 𝜕 𝑓 𝜕 𝜎 , ṗ ≥ 0, 𝑓 ≤ 0, ṗ𝑓 = 0, 𝑓 = 𝜎 𝑒𝑞 -𝑅(𝑝),
(2.13)

where 𝜎 𝑒𝑞 denotes the equivalent stress and 𝑅 the current yield stress. The equivalent plastic strain, 𝑝, is an internal state variable, which defines the level of the current yield stress. The von Mises stress is the standard equivalent stress and is a function of the second invariant of the stress deviator,

𝜎 𝑒𝑞 = 𝜎 𝑣𝑚 = √ 3𝐽 2 = 𝑌. (2.14)
where 𝑌 is the yield stress, achieved under uniaxial stress. Two other frequently employed isotropic yield criteria are the maximum shear and the maximum normal stress criterion,

𝜎 𝑡𝑟 = 𝜏 𝑚𝑎𝑥 = 1 2 |𝜎 1 -𝜎 3 |, 𝜎 𝑟𝑎 = 𝜎 1 .
(2.15)

In reference to their inventors, the criteria bear the names Tresca and Rankine [START_REF] Pilkey | Peterson's stress concentration factors[END_REF]. The current yield stress is defined as the sum of the initial yield stress and the strain hardening. Various formulations of the strain hardening have been proposed [START_REF] Larour | Strain rate sensitivity of automotive sheet steels: influence of plastic strain, strain rate, temperature, microstructure, bake hardening and pre-strain[END_REF].

The power hardening law (Ludwik's law) [START_REF] Ludwik | Elemente der technologischen Mechanik[END_REF] describes the hardening with

𝑅 = 𝜎 𝑦 + 𝐾𝑝 1/𝑀 , (2.16)
where 𝜎 𝑦 is the yield stress, 𝐾 the hardening modulus in MPa and 𝑚 the hardening exponent. The exponential hardening (Voce's law) [START_REF] Voce | The relationship between stress and strain for homogeneous deformation[END_REF] uses yield stress 𝜎 𝑦 , hardening stress 𝑅 ∞ and hardening strain coefficient 𝛾,

𝑅 = 𝜎 𝑦 + 𝑅 ∞ (1 -𝑒𝑥𝑝(-𝛾𝑝)).
(2.17)

The third law is a combined linear-exponential model, the extended Voce law with an additional hardening modulus 𝐶,

𝑅 = 𝜎 𝑦 + 𝑅 ∞ (1 -𝑒𝑥𝑝(-𝛾𝑝)) + 𝐶𝑝. (2.18)
More advanced multi-stage hardening laws exist, which emphasize the complex nature of strain hardening [START_REF] Lecarme | Void growth and coalescence in ductile solids with stage III and stage IV strain hardening[END_REF].

Influence of the triaxiality and of the Lode angle

The relationship between the macroscopic state of stress and the strain hardening is indirect and strongly depends on the microstructure. Therefore, various formulations for the equivalent stresses and the yield functions have been proposed depending on the material. For example some proposals introduce the influence of the hydrostatic pressure, 𝜎 ℎ , and/or the third deviatoric invariant, 𝐽 3 . Here, we examine the Hosford criterion and the linear Drucker-Prager criterion. More complex propositions comprise a sophisticated parameterization of the equivalent stress with triaxiality and Lode angle [START_REF] Bigoni | Yield criteria for quasibrittle and frictional materials[END_REF][START_REF] Wierzbicki | Calibration and evaluation of seven fracture models[END_REF][START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF][START_REF] Gao | On stress-state dependent plasticity modeling: Significance of the hydrostatic stress, the third invariant of stress deviator and the non-associated flow rule[END_REF]. Hosford criterion is a modification of the von Mises criterion [START_REF] Hosford | A generalized isotropic yield criterion[END_REF]. Only the exponent 𝑛 varies the form of the yield surface,

𝜎 ℎ𝑜 = [ (𝜎 1 -𝜎 2 ) 𝑛 + (𝜎 2 -𝜎 3 ) 𝑛 + (𝜎 1 -𝜎 3 ) 𝑛 2 ] 1 𝑛 . (2.19)
Compared to the von Mises criterion, it can be seen that the Hosford criterion includes a dependence on the Lode angle, cf. appendix A. The yield surface is equal to the von Mises criterion for 𝑛 = [2, 4]. For 𝑛 = 1 and 𝑛 → ∞ it approaches the Tresca criterion (Hosford, 1996). For 1 < 𝑛 < 4 and 𝑛 > 4 it lies between those two criteria and for 2 ≤ 𝑛 ≤ 4 the surface lies outside the von Mises cylinder [START_REF] Banabic | Sheet metal forming processes: constitutive modelling and numerical simulation[END_REF]. A relation of the parameter 𝑛 to the crystallographic structure of the material was suggested [START_REF] Logan | Upper-bound anisotropic yield locus calculations assuming <111>-pencil glide[END_REF]. For body centered cubic materials 𝑛 = 6 was recommended.

In its extended form, the Drucker-Prager criterion accounts for the influence of both triaxiality and third invariant. Neglecting the latter, the parameter 𝛽 introduces the sensitivity to pressure through a modified yield function (ABAQUS, 2011),

𝑓 = 𝜎 𝑣𝑚 -𝜎 ℎ 𝑡𝑎𝑛(𝛽) - 1 3 𝑡𝑎𝑛(𝛽)𝑅.
(2.20)

Anisotropic equivalent stress

The yield criteria discussed so far are applicable to isotropic materials. However, polycrystalline materials can have anisotropic properties. For TA6V such anisotropy was examined in [START_REF] Khan | Deformation induced anisotropic responses of Ti-6Al-4V alloy Part II: A strain rate and temperature dependent anisotropic yield criterion[END_REF]. In addition, an important tension-compression asymmetry was reported in [START_REF] Tuninetti | Anisotropy and tension-compression asymmetry modeling of the room temperature plastic response of Ti-6Al-4V[END_REF]. An anisotropic yield criterion was introduced by [START_REF] Hill | The mathematical theory of plasticity[END_REF],

𝜎 ℎ𝑖 = √ 𝐹 (𝜎 22 -𝜎 33 ) 2 + 𝐺(𝜎 33 -𝜎 11 ) 2 + 𝐻(𝜎 11 -𝜎 22 ) 2 + 2𝐿𝜎 2 23 + 2𝑀 𝜎 2 31 + 2𝑁 𝜎 2 12 , (2.21)
where, 𝐹, 𝐺, 𝐻, 𝐿, 𝑀, 𝑁 are coeffcients detailed in appendix H. In [START_REF] Karafillis | A general anisotropic yield criterion using bounds and a transformation weighting tensor[END_REF][START_REF] Bron | A yield function for anisotropic materials application to aluminum alloys[END_REF][START_REF] Wang | The equivalent plastic strain-dependent Yld2000-2d yield function and the experimental verification[END_REF][START_REF] Soare | Convex polynomial yield functions[END_REF][START_REF] Dunand | Experiments and modeling of anisotropic aluminum extrusions under multi-axial loading -Part I: Plasticity[END_REF] some modifications of the Hill criterion were proposed to be more suited for the materials that were studied.

Failure of metals

Fracture mechanics is amongst the earliest theories of failure. It assumes an existing crack and studies the conditions of propagation. It searches to calculate the driving force on a crack and to compare it to the material resistance to fracture [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF][START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF]. Continuous developments led to perception of the importance of notches in this context [START_REF] Irwin | Relation of stresses near a crack to the crack extension force[END_REF]. We will not discuss the development of this theory as no visible crack is accepted in the studied turbomachine.

The literature on metallic materials classifies metals into two groups with regard to the fracture mode: brittle and ductile. In the extreme, the term brittle specifies materials that fracture without any measurable plastic deformation. In what follows the term brittle refers to materials that fail without any previous strain softening or variation of the elastic modulus. In contrast, the term ductile refers to a material exhibiting macroscopic strain softening and a variation of the elastic modulus [START_REF] Hancock | On the mechanisms of ductile failure in highstrength steels subjected to multi-axial stress-states[END_REF], cf. figure 2.3.

Microstructural processes during metal failure Brittle fracture

At the microscopic level, one distinguishes between intra-and inter-granular mechanisms of failure. The dislocations accumulate at the grain boundaries, forming small voids. If a certain energy storage threshold is reached, these grain boundary voids grow and fracture is immediate by cleavage [START_REF] Besson | Mécanique non linéaire des matériaux[END_REF]. Frequently, the fracture runs along the grain boundaries, it is intergranular, cf. figure 2.7(a). However, a trans-granular fracture passing through the grains is possible as well, cf. figure 2.7(b). The surfaces are orientated normal to the maximum loading fracture first. In consequence, in a smooth tensile sample, perfect uniaxial loading leads to a macroscopic angle of the fracture surface of 0 ∘ .

Ductile fracture

The distortion of the crystal lattice around dislocations and other crystallographic defects invokes local stress concentrations. If the energy required to fracture the atomic bonds is beneath the energy of dislocation creation and movement, the atomic bonds breach. A micro void comes into existence. Consequently, the energy is "stored" in the void [START_REF] Besson | Mécanique non linéaire des matériaux[END_REF].

Figure 2.8 presents the three stages of void growth up to the establishment of a micro crack [START_REF] Garrison | Ductile fracture[END_REF]. This process is called rupture and can take place over a considerable part of the entire deformation range. There is no concentration to any particular location and nucleation takes place in the entire grain. Therefore, the fracture [START_REF] Flatten | Lokale und nicht-lokale Modellierung und Simulation thermomechanischer Lokalisierung mit Schädigung für metallische Werkstoffe unter Hochgeschwindigkeitsbeanspruchungen[END_REF].

surface is always transgranular, showing small dimples, which result from void growth. The grains whose primary glide planes coincide with the direction of maximum shear undergo plastic deformation first. Naturally, these grains are the first to undergo rupture and eventually fracture. In a uniaxial tensile test, many materials fracture in direction of the maximum shear (45 ∘ ).

The advances in X-Ray tomography allow for the observation on much smaller scales [START_REF] Maire | On the Application of X-ray Microtomography in the Field of Materials Science[END_REF][START_REF] Morgeneyer | Evolution of voids during ductile crack propagation in an aluminium alloy sheet toughness test studied by synchrotron radiation computed tomography[END_REF][START_REF] Benzerga | Ductile fracture by void growth to coalescence[END_REF]. The direct observation of void growth being possible, a growing research focuses on a better understanding of the causes that lead to the perceived macroscopic effects (i.e. the softening in the stress strain curve) in order to provide better models of the softening phase.

Constitutive modeling of material fracture Influence of the stress state

The final fracture strain differs, depending on the local stress state [START_REF] Bridgman | Studies in large plastic flow and fracture with special emphasis on the effects of hydrostatic pressure[END_REF][START_REF] Johnson | Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures[END_REF]. The earliest works describe the evolution of cylindrical and spherical holes in tensile test geometries [START_REF] Mcclintock | A Criterion for Ductile Fracture by the Growth of Holes[END_REF][START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF]. Similar observations exist for bars of increased thickness. Most authors report a decrease in the fracture strain with increasing triaxiality. For compression, a limit of triaxiality was reported below which fracture never occurs [START_REF] Bao | On the cut-off value of negative triaxiality for fracture[END_REF]. Recently conducted research, on the aluminum alloy 2024-T351, suggests that the phenomenon may be even more complex. At high triaxiality, the critical strain slightly increased [START_REF] Papasidero | Ductile fracture of aluminum 2024-T351 under proportional and non-proportional multi-axial loading: Bao-Wierzbicki results revisited[END_REF]. The influence of second characteristic invariant of the stress state, the Lode angle, has attracted attention [START_REF] Lode | The influence of the intermediate principal stress on yielding and failure of iron, copper and nickel[END_REF]. However, its influence seems restricted to low triaxialities [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF].

Studies of notch-dependent fracture of brittle materials address the stress sensitivity of failure not to a specific mechanism, but to the increased maximum principal stress, which acts in the notch vicinity [START_REF] Berto | Recent developments in brittle and quasi-brittle failure assessment of engineering materials by means of local approaches[END_REF].

The influence of the stress state on ductile fracture is the focus of a tremendous number of scientific publications [START_REF] Barsoum | Rupture mechanisms in combined tension and shear-micromechanics[END_REF][START_REF] Dunand | Hybrid experimental-numerical analysis of basic ductile fracture experiments for sheet metals[END_REF][START_REF] Danas | Influence of the Lode parameter and the stress triaxiality on the failure of elasto-plastic porous materials[END_REF][START_REF] Brünig | Micro-mechanical studies on the effect of the stress triaxiality and the Lode parameter on ductile damage[END_REF]. Under shear (i.e. low triaxiality), voids form but do not enlarge. The coalescence of two voids is the result of the nucleation of smaller intermediate voids.

In contrast, under tension (i.e. high triaxiality), the voids enlarge and eventually merge with neighboring voids [START_REF] Besson | Continuum models of ductile fracture: a review[END_REF]. High triaxiality may lead to the effect of transition between brittle and ductile fracture or a combined brittle/ductile failure (Lemaitre, 2005).

Damage mechanics

The concept of continuum damage mechanics was first introduced by [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF]. The main assumption is that the effect of voids and cracks can be averaged and characterized at the macroscopic level by damage indicators [START_REF] Besson | Mécanique non linéaire des matériaux[END_REF]. This approach has been notably developed for metal by Lemaitre and coworkers [START_REF] Chaboche | A CDM Approach of Ductile Damage with Plastic Compressibility[END_REF] and are often referred to as macroscopic damage models. A second class was proposed by Rice together with Tracey and Gurson. This approach aims at the description of the growth of voids in the elasto-plastic regime. The original Gurson model [START_REF] Gurson | Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I-Yield Criteria and Flow Rules for Porous Ductile Media[END_REF] is still the focus of many improvements.

Comparisons between the Gurson and Lemaitre models can be found in the literature [START_REF] Malcher | An assessment of isotropic constitutive models for ductile fracture under high and low stress triaxiality[END_REF]. Moreover, in [START_REF] Cao | A comparative study of three ductile damage approaches for fracture prediction in cold forming processes[END_REF] the different formulations of these models are also studied.

The Lemaitre model

The model formulation originally proposed by Lemaitre does not consider any microscopic effects (pores, inclusions) [START_REF] Chaboche | Anisotropic creep damage in the framework of continuum damage mechanics[END_REF][START_REF] Lemaitre | A Continuous Damage Mechanics Model for Ductile Fracture[END_REF]Lemaitre, 2005). The macroscopic stress and the modified yield function are given by

σ = 𝜎 1 -𝐷 , 𝑓 = 𝜎 𝑒𝑞 1 -𝐷 -𝑅(p), (2.22)
where σ, p are the effective stress and effective equivalent plastic strain, which are supposed to correspond to the stress and real equivalent plastic strain in a hypothetical undamaged reference configuration of the material. The damage, 𝐷, is a state variable that denotes the ratio of damaged to undamaged volume,

𝐷 = 𝑉 𝑑𝑎𝑚𝑎𝑔𝑒𝑑 𝑉 𝑢𝑛𝑑𝑎𝑚𝑎𝑔𝑒𝑑 .
(2.23)

The original proposition of the Lemaitre model describes the damage rate as

Ḋ = ṗ 1 1 -𝐷 ( -𝑌 𝑆 ) 𝑠 , (2.24)
where 𝑠, 𝑆 are material parameters and 𝑌 is the energy release rate density,

-𝑌 = 𝜎 2 𝑒𝑞 6𝐺(1 -𝐷) 2 + 𝜎 2 ℎ 2𝐾(1 -𝐷) 2 , (2.25)
and 𝐺, 𝐾 are the shear and the bulk modulus. A summary of possible improvements of the Lemaitre model concerning anisotropy, compressive loading and adaptation of the damage potential can be found in [START_REF] Bouchard | An enhanced Lemaitre model formulation for materials processing damage computation[END_REF]. Recently, a proposition with Lode angle influence was made (Cao et al., 2014a). A simplified form, which removes the coupling from the damage rate expression ex-ists [START_REF] Allix | The bounded rate concept: A framework to deal with objective failure predictions in dynamic within a local constitutive model[END_REF], (2.26) where 𝑝 0 denotes a threshold, which marks the onset of damage [START_REF] Lemaitre | Anisotropic damage law of evolution[END_REF].

𝐷 = ⟨ p -𝑝 0 𝑝 𝑐 -𝑝 0 𝑅 𝜈 ⟩ 𝛼 + , 𝐷 < 1,
The non-linearity parameter, 𝛼, reflects the effect of the damage rate increase towards fracture. The function 𝑅 𝜈 introduces the dependence on triaxiality. For a linear triaxiality independent damage law (𝛼 = 1, 𝑅 𝜈 = 1), the maximum equivalent plastic strain at fracture is 𝑝 𝑐 .

Formulations of the triaxiality function

The original form of the triaxiality function follows from derivations in a continuum thermodynamics framework [START_REF] Germain | Continuum Thermodynamics[END_REF]Lemaitre, 2005),

𝑅 𝜈 = 𝑅 𝜈 1 = 2 3 (1 + 𝜈) + 3(1 -2𝜈) (𝜂) 2 .
(2.27)

For an almost incompressible material, 𝜈 ≈ 0.5, the triaxiality dependence vanishes. Several authors have proposed corrections to improve the consistency with experimental results. For instance for a 16MnCr5 steel, propositions with differences between tension and compression exist [START_REF] Behrens | Verification of the damage model of effective stresses in cold and warm forging operations by experimental testing and FE simulations[END_REF],

𝑅 𝜈 2 = ⟨ 𝜂 + 2 3 ⟩ , 𝑅 𝜈 3 = 𝑒 2/3(𝜂-(1/3)) ,
(2.28)

In addition, proposition 𝑅 𝜈 2 has a negative cut-off value at 𝜂 = -2 3 , below which no damage occurs.

To account for more freedom in the choice of the function, a generalized form was proposed for an application to the high strength steel 18MnD5 [START_REF] Oueslati | Etude physique et modélisation numérique de l'endommagement d'alliages d'aluminium : application au sertissage[END_REF][START_REF] Winter | Identification of anisotropic elasto-plastic and damage behaviour from Arcan tests instrumented with Digital Image Correlation[END_REF],

𝑅 𝜈 = 𝐴 + 𝐵𝜂 2 .
(2.29) Similar to the other propositions, this function may be modified so that uniaxiality (𝜂 = 1/3) results in 𝑅 𝜈 = 1,

𝑅 𝜈 4 = (1 + 𝛽 𝐷 [𝜂 2 -(1/3) 2 ]).
(2.30)

For the steel 18MnD5, 𝛽 𝐷 = 6 was reported. Figure 2.9 presents a comparison of these propositions. The original formulation 𝑅 𝜈 1 is presented for U500 (𝜈 = 0.3) and TA6V (𝜈 = 0.32).

The Gurson-Tvergaard-Needleman model

The description of void growth and coalescence goes back to [START_REF] Brown | The Initiation and Growth of Voids at Second Phase Particles[END_REF]. A first theory was established by [START_REF] Gurson | Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I-Yield Criteria and Flow Rules for Porous Ductile Media[END_REF] and enriched with more features (pores, inclusions, triaxiality) by [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF][START_REF] Needleman | Mesh effects in the analysis of dynamic ductile crack growth[END_REF]. Hence, its

𝑓 = 𝐽 2 - 1 3 [ 1 + 𝑞 3 𝑓 *2 -2𝑞 1 𝑓 * 𝑐𝑜𝑠ℎ ( 𝑞 2 3𝜎 ℎ 2𝜎 𝑦 )] 𝜎 2 𝑦 .
(2.31)

The porosity 𝑓 * is the volume of the voids divided by the reference volume

𝑓 * = 𝑉 𝑣𝑜𝑖𝑑𝑠 𝑉 𝑟𝑒𝑓 .
(2.32)

For more details of the calculation of 𝑓 * cf. [START_REF] Malcher | An assessment of isotropic constitutive models for ductile fracture under high and low stress triaxiality[END_REF]. The analysis of the stress state relied for a long time on analytical solutions. For axisymmetric samples, formulas predict the strain fields in elasticity [START_REF] Neuber | Theory of Stress Concentration for Shear-Strained Prismatical Bodies With Arbitrary Nonlinear Stress-Strain Law[END_REF]. To correct the measured data, other propositions use stress concentration factors [START_REF] Pilkey | Peterson's stress concentration factors[END_REF] or geometrical considerations to define a global triaxiality value [START_REF] Bridgman | Studies in large plastic flow and fracture with special emphasis on the effects of hydrostatic pressure[END_REF]. Using FE calculations, a correction of these latter formulas was proposed for round notched bars and flat grooved plates [START_REF] Bao | On the cut-off value of negative triaxiality for fracture[END_REF][START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF],

Existing tests to examine the stress state influence

𝜂 𝑟𝑜𝑢𝑛𝑑 = 1 3 + 1.4𝑙𝑛 ( 1 + 𝑎 2𝑅 ) , 𝜂 𝑓 𝑙𝑎𝑡 = + √ 3 3 [ 1 + 2𝑙𝑛 ( 1 + 𝑡 4𝑅 )] . (2.33)
Where, 𝑅 is the radius of the notch and 𝑎 the axisymmetric sample's thinnest diameter and 𝑡 the thickness in the case of plates To describe the stress state in any arbitrary shape, a general method based on FE calculations was proposed (Bao et al., 2004b;[START_REF] Luo | Experiments and modeling of anisotropic aluminum extrusions under multi-axial loading -Part II: Ductile fracture[END_REF]. An average value of the triaxiality is defined over the time,

𝜂 𝑎𝑣 = 1 𝑝 𝑚𝑎𝑥 ˆ𝑝𝑚𝑎𝑥 0 𝜂 d𝑝.
(2.34)

Other quantities, as for example the average Lode angle, are obtained analogously. Table 2.2 presents an overview of the average triaxiality values, reported for different geometries.

Average values allow a straightforward reduction of huge amounts of data. On the other hand, they conceal the complex shape of the triaxiality field. The variation in the field can be analyzed by tracing single elements. For simple geometries, such as the round notched, the variation depends largely on the radius [START_REF] Dunand | Hybrid experimental-numerical analysis of basic ductile fracture experiments for sheet metals[END_REF][START_REF] Erice | Flow and fracture behaviour of FV535 steel at different triaxialities, strain rates and temperatures[END_REF][START_REF] Fourmeau | Anisotropic failure modes of high-strength aluminium alloy under various stress states[END_REF]. In addition, the representation of the entire stress-space becomes more frequent [START_REF] Dunand | Hybrid experimental-numerical analysis of basic ductile fracture experiments for sheet metals[END_REF]Kim et al., 2014a). A discussion of the stress state analysis methods is given in the next chapter, cf. section 3.1. [START_REF] Wierzbicki | Calibration and evaluation of seven fracture models[END_REF][START_REF] Bao | On the cut-off value of negative triaxiality for fracture[END_REF] Flat, hole 0.31-0.56 [START_REF] Wierzbicki | Calibration and evaluation of seven fracture models[END_REF][START_REF] Luo | Experiments and modeling of anisotropic aluminum extrusions under multi-axial loading -Part II: Ductile fracture[END_REF][START_REF] Lou | Prediction of ductile fracture for advanced high strength steel with a new criterion: Experiments and simulation[END_REF] Flat, notch 0.4-0.76 [START_REF] Wierzbicki | Calibration and evaluation of seven fracture models[END_REF][START_REF] Luo | Experiments and modeling of anisotropic aluminum extrusions under multi-axial loading -Part II: Ductile fracture[END_REF][START_REF] Hammer | Plastic Deformation and Ductile Fracture of Ti-6Al-4V under Various Loading Conditions[END_REF] 2.7 Graphical representation of the fracture locus For increasing triaxiality, the equivalent plastic strain at fracture decreases. Earlier findings reported a cut off value at 𝜂 = -1 3 [START_REF] Bao | On the cut-off value of negative triaxiality for fracture[END_REF], below which fracture never occurs. In contrast, more recent findings provide a fracture locus in the range -0.7 ≤ 𝜂 ≤ 0.7 [START_REF] Papasidero | Ductile fracture of aluminum 2024-T351 under proportional and non-proportional multi-axial loading: Bao-Wierzbicki results revisited[END_REF]. In figure 2.11(b) the Lode angle dependence was included. The locus for the fracture strain 𝜀 𝑓 is then three-dimensional.

Tensional tests of TA6V at low temperatures led to brittle fracture [START_REF] Diiorio | Plastic deformation, damage and rupture of PM Ti-6Al-4V at 20K under monotonic loading[END_REF]. Figure 2.11(c) presents the identified fracture locus. The valid domain is limited by the maximum principal stress and the maximum equivalent plastic strain. The interpretation of this locus is more complicated, namely, fracture at low strains depends only on the principal stress, while for higher strains it depends also on the maximum equivalent plastic strain.

The central issue about these loci is how the data are obtained from the experiments. For the presented triaxiality -equivalent plastic strain locus, the FE averaging method from equation (2.34) was used. For the creation of the maximum principal stress -equivalent plastic strain locus, the elemental FE values at the centre of the sample were used.

Identification of the constitutive law

Standard approach

The constitutive model, that is the relationship between stress and strain, needs to be identified from experiments. Moreover, a fracture model requires the knowledge of the maximum experimental strain or stress.

The standard technique to measure the relationship between stress and strain consists of the use of a smooth sample under uniaxial tension. The actual measurement of the strains in the tensile test in many cases still relies on strain gages or extensometers.

Assuming a uniform stress distribution in a sample of initial cross-section, 𝑆 0 , and together with the measured force 𝐹 , one obtains the engineering stress, 𝜎 𝐸 = 𝐹 𝑆 0 . For extensometers, the engineering strain follows from the measured sample elongation 𝛥𝑙 and the initial length 𝑙 0 to 𝜀 𝐸 = 𝛥𝑙 𝑙 0 . Assuming incompressibility, the true stress 𝜎 𝑇 and true strain 𝜀 𝑇 are then,

𝜀 𝑇 = ln (1 + 𝜀 𝐸 ) , 𝜎 𝑇 = 𝜎 𝐸 (1 + 𝜀 𝐸 ).
(2.35)

Measurement of heterogeneous strain fields

Most strain gages are only certified up to 10 % level (Micro-Measurements, 2007). In addition, strain gage and extensometer do not deliver representative strain data if the strain field is not homogeneous. In particular, no reliable information regarding the maximum strain state is in general available. For a smooth sample, heterogeneous strain fields occur as soon as a non-uniform reduction happens corresponding to the beginning of global necking. To analyze those situations several authors have provided post-necking corrections based on FE analyses [START_REF] Rosa | Postnecking elastoplastic characterization: Degree of approximation in the bridgman method and properties of the flow-stress/true-stress ratio[END_REF][START_REF] Zhang | Determining material true stress-strain curve from tensile specimens with rectangular cross-section[END_REF][START_REF] Cabezas | Experimental and numerical analysis of the tensile test using sheet specimens[END_REF][START_REF] Mirone | A new model for the elastoplastic characterization and the stress-strain determination on the necking section of a tensile specimen[END_REF][START_REF] Joun | Finite element analysis of tensile testing with emphasis on necking[END_REF][START_REF] Joun | A new method for acquiring true stress-strain curves over a large range of strains using a tensile test and finite element method[END_REF]. Other approaches exist based on the measurement of the critical section [START_REF] Fourmeau | Anisotropic failure modes of high-strength aluminium alloy under various stress states[END_REF]. Direct exploitations of DIC displacement fields were also reported [START_REF] Hoffmann | Determination of True Stress-Strain-Curves and Normal Anisotropy in Tensile Tests with Optical Strain Measurement[END_REF][START_REF] Han | The application of digital image techniques to determine the large stress-strain behaviors of soft materials[END_REF][START_REF] Hopmann | Determination of strain rate dependent material data for FEA crash simulation of polymers using digital image correlation[END_REF][START_REF] Tuninetti | Anisotropy and tension-compression asymmetry modeling of the room temperature plastic response of Ti-6Al-4V[END_REF]. This shows the need to obtain local information on the state of strains and the interest in the use of field measurement.

This aspect is even more crucial when trying to identify failure models as highlighted in [START_REF] Roux | On the interest of using full field measurements in ductile damage model calibration[END_REF]. Local information as those associated with the use of X-ray images allow obtaining more precise indication on the nature of the failure process (Cao et al., 2014b).

Digital image correlation techniques

2.9.1 The principle of image-based field measurements Digital image correlation Different optical techniques are available to measure displacement fields (Rastogi et al., 2012). Digital Image Correlation (DIC) has proven to be one of the most versatile methods. Another major advantage is the applicability to "any" scale of observation [START_REF] Grediac | Full-Field Measurements and Identification in Solid Mechanics[END_REF]Hild et al., 2012b;[START_REF] Sutton | Computer Vision-Based, Noncontacting Deformation Measurements in Mechanics: A Generational Transformation[END_REF]. DIC makes use of a distinctive pattern on the observed surface to measure the displacement fields between two or more recorded images.

The pattern can be the natural texture or consists of the application of paint speckles. The measurement of the displacement field relies on the premise of the conversation of gray level in the examined domain [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF]. Starting from its origins as an academic technique, there is now widespread use in academia and industry [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF][START_REF] Chu | Applications of digital-imagecorrelation techniques to experimental mechanics[END_REF][START_REF] Hild | Digital image correlation: From measurement to identification of elastic properties -A review[END_REF][START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF].

Digital image correlation consists in the measurement of displacement fields between a picture 𝑓 in the reference configuration and a series of pictures 𝑔 in the deformed configuration [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF]Hild et al., 2012a;[START_REF] Sutton | Computer Vision-Based, Noncontacting Deformation Measurements in Mechanics: A Generational Transformation[END_REF], where 𝑓 (𝑥) and 𝑔(𝑥,𝑡) describe measured fields of gray levels. 𝑥 is any pixels in the images and 𝑡 the considered time. The transformation of the deformed series 𝑔(𝑥,𝑡) with the spatio-temporal displacement field 𝑢(𝑥,𝑡) leads to the definition of the residual image series,

𝑟(𝑥,𝑡) = |𝑔(𝑥 + 𝑢(𝑥,𝑡),𝑡) -𝑓 (𝑥)|.
(2.36)

The problem of the recovery of the displacement field from the gray level is an ill-posed problem, as such it has an infinite number of solutions. To obtain a displacement field 𝑢(𝑥,𝑡) one then has to formulate the problem as an optimization one. One possible cost function that makes use of the residual image series 𝑟(𝑥,𝑡) reads

𝜂 2 (𝑡) = 1 |𝛺| ∑ 𝛺 (𝑔(𝑥 + 𝑢(𝑥,𝑡),𝑡) -𝑓 (𝑥)) 2 , (2.37)
where 𝛺 is the Region Of Interest (ROI).

The differences between the methods based on such cost function are related to the space within which the displacement field is searched and the regularization technique used within the method.

Acquisition noise affects both images, 𝑓 and 𝑔. A normalization with the sensor noise level (i.e. standard deviation 𝛾 𝑓 ) was proposed [START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF],

𝜒 2 (𝑡) = 𝜂 2 (𝑡) 2𝛾 2 𝑓 .
(2.38)

The expexted minimum value is 𝜒 = 1. Besides the displacement field 𝑢(𝑥,𝑡), the analysis of the residual images allows the gain of further information. For example, cracks and their evolution are clearly detectable through a local increase in residuals [START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF]. The speckle pattern becomes discontinuous and no optimum displacement can be found if one searches for continuous displacement fields. To illustrate the technique, a simple example of an image with a prescribed deformation is used, cf. figure 2.12(a, b). Figures 2.12(c, d) present the associated field of residuals for a displacement field equal to zero and for an optimal displacement field leading to zero residuals (the green line corresponds to the marking of the speckles). For convenience, the residual image is expressed with respect to the dynamic

Finite element updating methods

In order to obtain the parameters of a material model from the displacement field or the associated strain field, one has to define an identification procedure [START_REF] Avril | Overview of Identification Methods of Mechanical Parameters Based on Full-field Measurements[END_REF]. A popular one is based on the coupling of FE and DIC, known as Finite Element Model Updating (FEMU). FEMU compares the field obtained via DIC with its counterpart from FE analysis. A cost function between these two fields can be established. The minimization of this function leads to optimum values of the model parameters that describe the experiment. Here what is called an optimum is related to the cost function and again several proposals have been made in the literature. Various identification techniques have been applied in elasticity [START_REF] Grédiac | Applying the Virtual Fields Method to the identification of elasto-plastic constitutive parameters[END_REF][START_REF] Gras | Identification of the Out-of-Plane Shear Modulus of a 3D Woven Composite[END_REF], elastoplasticity [START_REF] Grédiac | Applying the Virtual Fields Method to the identification of elasto-plastic constitutive parameters[END_REF][START_REF] Cooreman | Elasto-plastic material parameter identification by inverse methods: Calculation of the sensitivity matrix[END_REF][START_REF] Robert | Identification of hardening parameters using finite element models and full-field measurements: some case studies[END_REF][START_REF] Réthoré | Robust identification of elasto-plastic constitutive law parameters from digital images using 3D kinematics[END_REF]Kim et al., 2014a;[START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF], as well as for damage [START_REF] Claire | Identification of damage fields using kinematic measurements[END_REF][START_REF] Chalal | Direct identification of the damage behaviour of composite materials using the virtual fields method[END_REF][START_REF] Réthoré | A fully integrated noise robust strategy for the identification of constitutive laws from digital images[END_REF][START_REF] Wu | Strain Localisation and Damage Measurement by Full 3D Digital Image Correlation: Application to 15-5PH Stainless Steel: Measuring Damage in Ductile Materials[END_REF][START_REF] Bouterf | Damage law identification from full field displacement measurement: Application to four-point bending test for plasterboard[END_REF][START_REF] Roux | On the interest of using full field measurements in ductile damage model calibration[END_REF]. However, the use of DIC displacements in the cost function implies several difficulties. In fact, the resolution of the DIC mesh is limited by the quality of the speckle pattern. Gray level distributions in the elements must be distinguishable. The element size restriction of the DIC method means that the displacement information is only available in a coarse mesh [START_REF] Besnard | Finite-Element" Displacement Fields Analysis from Digital Images: Application to Portevin-Le Châtelier Bands[END_REF]. If strain localization is very intense, the use of a coarse mesh results in an averaging of the localized displacements. Furthermore, if the same mesh is used in DIC and FE, the simulation may spatially not be converged. The use of a mesh, which ensures numerical convergence, requires the interpolation between the DIC mesh data and the FE mesh. The Integrated Digital Image Correlation (I-DIC) approach has been proposed to overcome those difficulties. Under a certain number of conditions it was shown that FEMU and I-DIC are equivalent [START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF].

Integrated digital image correlation

In order to simplify and automate the procedure, I-DIC identifies the sought parameters not from displacement or strain fields, but directly from the gray level images. The size of the elements is in theory without limits [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF] l. I-DIC transforms the recorded images with FE-based displacement field 𝑢 𝑐 .

Let 𝑝 denote the vector gathering all unknown material parameters. Then, the I-DIC approach minimizes, globally over space and time, the functional with respect to the sought material parameters,

𝑝 opt = argmin 𝑝 𝜒 2 tot (𝑝).
(2.40)

The total residual function consists of the global equilibrium gap 𝜒 2 𝐹 and global correlation residuals 𝜒 2 𝑐 ,

𝜒 2 tot (𝑝) = 1 2 𝜒 2 𝑐 (𝑝) + 1 2 𝜒 2 𝐹 (𝑝).
(2.41)

The correlation residuals, equation (2.37), depend on the computed displacement and hence on 𝑝,

𝜒 2 𝑐 (𝑝) = 1 𝑛 ∑ 𝑛 𝜒 2 𝑐 (𝑡,𝑝), 𝜒 2 𝑐 (𝑡,𝑝) = 𝜂 2 𝑐 (𝑡,𝑝) 2𝛾 2 𝑓 , 𝜂 2 𝑐 (𝑡,𝑝) = 1 |𝛺| ∑ 𝛺 ( 𝑔(𝑥 + 𝑢 𝑐 (𝑥,𝑡,𝑝),𝑡) -𝑓 (𝑥) ) 2 .
(2.42)

In addition, it is essential to add to the previous functional the equilibrium gap in order to take into account the information about the load level [START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF],

𝜒 2 𝐹 (𝑝) = (𝐹 𝑚 -𝐹 𝑐 (𝑝)) 𝑡 (𝐹 𝑚 -𝐹 𝑐 (𝑝)) 𝛾 2 𝐹 𝑛 , (2.43)
where 𝐹 𝑚 gathers all measured load levels for all considered images 𝑛. 𝐹 𝑐 is the vector gathering all computed resultant forces. 𝛾 𝐹 is the standard resolution of the load measurement device. The minimization of the correlation residuals allows the determination of the constitutive parameters of the chosen constitutive law.

Procedures of the digital image correlation methods General considerations

As already mentioned the correlation quality depends on many parameters. The latter ones comprise first experimental aspects:

• How is the alignment of the cameras to the sample?

• How is the lightning of the observed zone?

• Which are the optical equipment properties?

• Which are the properties of the speckle pattern? Which resolution does it allow?

• Is the correlation based on one camera or a stereo-correlation system?

For 2D correlation techniques, the tilt angle between the sample surface and the camera plane was studied by [START_REF] Lava | Error estimation in measuring strain fields with DIC on planar sheet metal specimens with a non-perpendicular camera alignment[END_REF]. The influence of out of plane movement for 2D DIC was studied by a comparison of 2D DIC and 3D stereo DIC [START_REF] Sutton | The effect of out-of-plane motion on 2D and 3D digital image correlation measurements[END_REF]. However, recent developments exploit the change in the image texture to extract 3D displacements from one single camera [START_REF] Réthoré | 3D displacement measurements using a single camera[END_REF]. The application of this method decreases the measurement error.

The literature is rich in studies concerning speckle pattern size and distribution. There are experimental studies [START_REF] Lecompte | Quality assessment of speckle patterns for digital image correlation[END_REF][START_REF] Triconnet | Parameter choice for optimized digital image correlation[END_REF] or studies based on computer generated images [START_REF] Bornert | Assessment of Digital Image Correlation Measurement Errors: Methodology and Results[END_REF]. The quality of the speckle pattern determines the minimum element size. Criteria for the assessment exist (Hild et al., 2012b). However, these studies concern the initial pattern quality. Little findings exist that report the pattern evolution and the eventual degradation under very large strains.

The second aspect is the correlation algorithm itself [START_REF] Schreier | Systematic errors in digital image correlation caused by intensity interpolation[END_REF][START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF]Hild et al., 2012a). The main issues can be summarized as follows:

• Are the images analyzed with locally independent windows or are global characteristics considered?

• What is the form and size of the windows?

• What type of interpolation is used for displacements and how are strains calculated?

Global digital image correlation

The displacement basis for 𝑢(𝑥,𝑡) is selected to minimize 𝜂 2 𝑐 with respect to the unknown amplitudes. This displacement basis can be local, as in most commercial packages [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF][START_REF] Chu | Applications of digital-imagecorrelation techniques to experimental mechanics[END_REF][START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF]. These local approaches evaluate every window independently. Global continuous displacement fields found in finite element discretizations are another choice [START_REF] Broggiato | Adaptive image correlation technique for full-field strain measurement[END_REF][START_REF] Sun | Finite-element formulation for a digital image correlation method[END_REF][START_REF] Besnard | Finite-Element" Displacement Fields Analysis from Digital Images: Application to Portevin-Le Châtelier Bands[END_REF]. Their application leads to lower uncertainties in the measured displacement fields when compared with local DIC approaches with the same local interpolations, that is, the same mesh size (Hild et al., 2012b). In addition, these global methods facilitate the implementation of a multi-scale approach. The displacement is first identified on a coarse scale, using superelements. Once this scale is converged, the superelement displacements serve as initial values for the lower scales. Several scales may be used, the smallest scale is equal to the local elements. Again, better and faster convergence is the result [START_REF] Hild | Multiscale displacement field measurements of compressed mineral-wool samples by digital image correlation[END_REF][START_REF] Besnard | Finite-Element" Displacement Fields Analysis from Digital Images: Application to Portevin-Le Châtelier Bands[END_REF]. In addition, the interest of using information at different scales and therefore different DIC meshes is the possibility to use information of cameras with different resolutions as done in [START_REF] Passieux | Multiscale Displacement Field Measurement Using Digital Image Correlation: Application to the Identification of Elastic Properties[END_REF].

Figure 2.13 provides an illustration of the principle of local and global approaches. In practice, the correlation of the image series is performed in sequential order. For a better readability, the time dependence of the terms is omitted. For every element, the displacement of any pixels is interpolated by the nodal displacements 𝑢 𝑖 and the shape functions 𝜓 𝑖 (chosen according to T3, Q4 interpolation),

𝑢(𝑥) = ∑

𝑖 𝑢 𝑖 𝜓 𝑖 (𝑥).

(2.44)

The interest of using improved quadrature rules for FE DIC to take into account pixels belonging to more than one element has been shown [START_REF] Pierré | Unstructured finite element-based digital image correlation with enhanced management of quadrature and lens distortions[END_REF]. Together with equation (2.37) the matrix form of the linearized minimization problem is

𝑀 𝛿𝑢 = 𝑏, (2.45)
where 𝛿𝑢 is the vector of nodal displacement increments. Using the global shape function vector 𝜓, the matrix 𝑀 is assembled from the image gradients and the vector 𝑏 from image gradient and residual image 𝑟,

𝑀 𝑖𝑗 = 1 |𝛺| ∑ 𝛺 ( 𝜓 𝑖 (𝑥) • ∇𝑓 (𝑥) ) ( ∇𝑓 (𝑥) 𝑇 • 𝜓 𝑗 (𝑥) ) d𝑥, 𝑏 𝑖 = 1 |𝛺| ∑ 𝛺 ( 𝜓 𝑖 (𝑥) 𝑇 • ∇𝑓 (𝑥) ) 𝑟(𝑥) d𝑥.
(2.46)

The Newton-Raphson scheme allows an efficient solution of this system to determine the displacement field [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Bruck | Digital image correlation using Newton-Raphson method of partial differential correction[END_REF]Hild et al., 2012a).

Integrated digital image correlation

The residual functional 𝜒 2 tot is highly non-linear with respect to 𝑝. Similar to DIC, a modified Newton-Raphson method is used for the minimization scheme. To update the current estimate of the parameters, the sensitivities of the computed displacement field 𝑢 𝑐 (𝑝) and load 𝐹 𝑐 (𝑝) have to be assessed iteratively,

𝑢 𝑐 (𝑥,𝑝 (𝑖) ) =𝑢 𝑐 (𝑥,𝑝 (𝑖-1) ) + 𝜕𝑢 𝑐 𝜕𝑝 (𝑥,𝑝 (𝑖-1) )𝛿𝑝, 𝐹 𝑐 (𝑝 (𝑖) ) =𝐹 𝑐 (𝑝 (𝑖-1) ) + 𝜕𝐹 𝑐 𝜕𝑝 (𝑝 (𝑖-1) ),
(2.47)

where 𝑝 (𝑖-1) is the set of parameters at iteration 𝑖 -1, and 𝛿𝑝 the parameter increment.

The computation of the sensitivity fields,

𝑆 (𝑖) 𝑢 = 𝜕𝑢 𝑐 /𝜕𝑝, (2.48)
is either analytical [START_REF] Cooreman | Elasto-plastic material parameter identification by inverse methods: Calculation of the sensitivity matrix[END_REF][START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF] or numerical. The numerical sensitivities are achieved by computing the displacement fields for small variations of each parameter [START_REF] Mathieu | Analyse de la tenue mécanique d'un liner en titane : apport des mesures de champs cinématiques[END_REF]. The sensitivity matrix allows the projection of the DIC quantities 𝑀 and 𝑏 in the space of the parameters, 𝑁 (𝑖) = 𝑆 (𝑖) 𝑢 .𝑀 .𝑆 (𝑖) 𝑢 , 𝐵 (𝑖) = 𝑆 (𝑖) 𝑢 𝑏.

(2.49)

The DIC problem, equation (2.45), with its many degrees of freedoms is reduced to a much more manageable problem, 𝑁 (𝑖) 𝛿𝑝 = 𝐵 (𝑖) .

(2.50)

The force sensitivities 𝑆 𝐹 = 𝜕𝐹 𝑐 /𝜕𝑝 are obtained analogously and the combined problem is solved [START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF].

The strong integration of the images in the I-DIC identification is a direct extension of DIC. DIC is frequently perceived as a measurement technique. In the same manner, I-DIC may be considered as a measurement and identification technique. Further, I-DIC may be considered as a measurement method that allows the evaluation of the local stress state. The residual image indicates the local errors.

The boundary conditions can be included in the parameters or taken from a DIC correlation. The application of DIC and I-DIC in the same procedure results in a multiscale like approach. The boundary displacement can be identified with a very coarse mesh. The refined scale is entirely performed using I-DIC.

Due to the complexity, the I-DIC method is prone to suffer from several error sources. One may separate the errors in two groups, namely, input and method errors. Input errors relate to the prescribed model parameters. These are the uncertainties of the geometry, the boundary conditions or the material model. Furthermore, if more than one camera is used, the relative positioning of the cameras must be accounted for. A CAD-based positioning procedure was developed [START_REF] Beaubier | CAD-based calibration and shape measurement with stereoDIC: Principle and application on test and industrial parts[END_REF][START_REF] Dufour | CAD-based Displacement Measurements with Stereo-DIC: Principle and First Validations[END_REF].

Method errors are due to the finite element discretization and errors in the solution of equation (2.50). The finite element errors can be assessed with error criteria and convergence criteria [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineerng analysis[END_REF][START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF][START_REF] Pled | Vers une stratégie robuste et efficace pour le contrôle des calculs par éléments finis en ingénierie mécanique[END_REF]. An assessment of I-DIC can be done with FE-based virtual images [START_REF] Mathieu | Analyse de la tenue mécanique d'un liner en titane : apport des mesures de champs cinématiques[END_REF][START_REF] Neggers | Time-resolved integrated digital image correlation[END_REF]. These procedures allow every error source to be addressed separately.

Conclusion of the state of the art

The two materials of interest, TA6V and Udimet 500 are reported to both exhibit plastic strains greater than 10 %. Literature findings indicate that fracture can depend severely on the stress state and in particular the triaxiality. In addition, the application of local stress-state based criterion to disk calculations showed a significant influence on the burst prediction. To improve such predictions it is necessary to study the material at the stress state that occurs in the disks.

The burst prediction conducted so far were using constitutive models identified on smooth tensile samples. Two possible limitations about this type of experiments have to be considered. The first one is the localization of stress and strain from the onset of necking. As a consequence, an identification based on global measurements can be conducted only up to a limited level of strain. The second limitation concerns the state of triaxiality in these experiments that is not representative of those in the disks.

To overcome these limitations we decided to conduct experiments on non-smooth tensile sample geometries allowing much more representative levels of triaxiality to be obtained. The use of DIC and I-DIC techniques will, presumably, enable for the acquisition of as much as possible relevant local information with regard to the identification of the constitutive behavior and of the local failure criteria.

CHAPTER 3

Identification of plasticity at high strain levels and for heterogeneous samples using I-DIC One axis of the thesis concerns the identification of the material behavior up to fracture under multiaxial loading for situations representative of those in service. The triaxiality was identified as a possible important factor. Therefore, special attention was made to the state of stress in the tested geometries. An analysis of disk geometries clarified the range of interest for the triaxiality. This led to the definition of two sample geometries (i.e. thin and thick, cf. section 3.1).

The experimental procedures are detailed in section 3.2. The experimental results and their exploitation are presented separately for each material in sections 3.3 and 3.4.

Design of uniaxial tests at high triaxiality values

Triaxiality in turbomachine disks

Qualitative analysis

The range of triaxiality, which occurs in the test, should correspond to the predominant values in turbomachine disks. To this end, two different disk geometries were analyzed at critical rotational speeds. The two examined disks present characteristically different geometries. A "non-drilled" disk 𝐴 with one central hole and a "drilled" disk 𝐵 with five decentralized holes, cf. figure 3.1(a, b). Disk 𝐴 was simulated under the assumption of axisymmetry with a special formulation for the non-axisymmetric outer region. The drilled disk was simulated with a three-dimensional analysis. Only an angular sector of the disk was modeled with cyclic symmetry.

Both simulations were performed using the ANSYS software with the plastic behavior of a nickel-base alloy. Non-linear geometrical effects were considered. The axisymmetric mesh is made of quadratic quadrilateral and triangular elements. For the three-dimensional mesh, quadratic hexahedral and tetrahedral elements were used.

Figure 3.1(c-f) shows the fields of normalized equivalent plastic strain and triaxiality obtained from these simulations. Normalized equivalent plastic strain is defined with the maximum value as 𝑝 ′ = 𝑝/𝑝 𝑚𝑎𝑥 .

Both disks possess areas of strain concentrations especially around holes. For the nondrilled disk, the maximum of the equivalent plastic strain is located around the central hole corresponding to the area marked with (1) in figure 3.1(c). For the drilled disk, the equivalent plastic strain is high along the radial plane passing through the hole at points (3, 4). Additional critical areas exist. The neck of disk 𝐴 shows a small area of high strain at point (2). For disk 𝐵, the areas (5, 6) are related to sharp changes in the geometry. These lead to stress concentrations and high strains. High plastic strains are a problem as such. Moreover, it follows from the literature review that the combination with high triaxiality adversely affects the failure. high triaxiality is found at point (2) and at point (7) for the non-drilled disk. For disk 𝐵, the values are high at the locations (5, 8 and 9).

In conclusion, the areas ( 2) and ( 5) present a critical combination of equivalent plastic strain and triaxiality, both roughly around 𝜂 ≈ 0.75.

Quantitative analysis

Two key aspects are to be taken into account for the design of the experiments. The first one concerns the level of strain. The second aspect concerns the level of triaxiality. To compare different specimen an indicator of the triaxiality level is used,

𝜂 𝑎𝑣 = 1 𝑝 𝑚𝑎𝑥 𝑉 ˆ𝑅𝑂𝐼 ˆ𝑝𝑚𝑎𝑥 0 𝜂 d𝑝 d𝑉. (3.1)
This integral is evaluated on FE calculations but two issues have to be considered. First, while the selection of the Region of Interest (ROI) in a tensile geometry is rather straightforward, this is not the case in the disk. In the latter, the values are very heterogeneous and the influence of the extrema is diminished if one considers an extended ROI. For example if one consider a ROI corresponding to the entire disk, one obtains 𝜂 𝑎𝑣 1 = 0.57 for disk 𝐴 and 𝜂 𝑎𝑣 1 = 0.54 for disk 𝐵, which is far from the maximum values, cf. appendix C.

In order to have a more precise way of analysis, the results of all elements were plotted in a diagram of triaxiality and normalized equivalent plastic strain, cf. figure 3.2. This representation allows for the precise identification of the maximum points for each disk. For disk 𝐴, the total maximum strain is identifiable at a stress triaxiality close to 𝜂 = 0.4. Another local maximum exists at 𝜂 = 0.7. Based on the assumption that strains become more critical at a higher stress triaxiality, one should also consider the presence of intermediate strain levels with a high level of triaxiality close to 𝜂 = 0.8. For disk 𝐵, maxima are identified for 𝜂 = 0.4 and 𝜂 = 0.75. From these analyses, it is clear that there is no unique critical value. Instead, it seems more relevant to define a target triaxiality range. From the above considerations this target range should extend from 𝜂 = 0.4 to 𝜂 = 0.75. For an even more detailed analysis cf. appendix D

Determination of representative sample geometries

A first concern is to be able to measure the area where fracture occurs, this is a priori achieved by introducing important geometrical strain concentrations. The second one is to obtain representative levels of triaxiality in this area.

Definition of sample geometries

Several geometries allow testing to be in the range of 0.4 ≤ 𝜂 ≤ 0.75:, namely, round notched samples, flat samples with notch or hole and biaxial samples, cf. section 2.6. In order to master as well as possible the loading conditions, only uniaxial tension tests are considered. The use of 2D-DIC requires flat surfaces or a very small curvature. Therefore, only flat geometries were considered. Within these, notches allow higher triaxialities than holes. With the employed equipment, the sample maximal cross-section was set to 8 x 7 mm, cf. section 3.2.

Taking these considerations into account led to the design of a sample of 7 mm in thickness. In order to test the procedure on a sample with low three-dimensional effects and possibly higher strains it was also decided to test thin samples (0.7 mm). Figure 3.3 depicts the final designs. in the thickness direction. However, this plane stress state leads to an overall decrease in triaxiality. Maximum equivalent plastic strain shifts to 𝜂 = 0.38, cf. figure 3.4(d).

Pre-test analysis of sample geometries made of U500 For U500, the same geometry leads to completely different fields, cf. figure 3.5. For the thick sample the maximum equivalent plastic strain corresponds to 𝜂 = 0.42, cf. figure 3.5(c).

Estimates for the thin sample show more strain concentration. The total maximum strain at 𝜂 = 0.38 is far greater than most other strains, cf. figure 3.5(d).

Differences between the local strain and stress field of TA6V and U500 The differences between the stress and strain fields are due to the different constitutive behavior. The Young's modulus of U500 is two times that of TA6V while the yield stress of U500 is lower than that of TA6V. Maybe more importantly the two materials exhibit significant yield behavior. Figure 3.6 shows the normalized yield curves of both materials. The hardening 𝑅 is expressed relative to the maximum equivalent stress 𝑅 𝑚𝑎𝑥 . The TA6V alloy is characterized by a very small strain hardening slope, contrary to U500, which demonstrates very large capabilities of strain hardening. The implication of those differences are further discussed with the fields identified from the experiments on TA6V in section 3.3.2 and on U500 in section 3.4.3. 3.2 Procedures for the testing with digital image correlation methods

Experimental setup Tensile test equipment and sample preparation

The tests were conducted under uniaxial tension with an MTS 100KN machine. To guarantee constant external load conditions up to fracture, the experiments are displacement controlled. The inter-clamp strain rate was set to ε = 0.001 1 𝑠 . In addition, some samples were subjected to unloading-reloading cycles to study the potential evidence of damage.

The samples were fabricated from an extruded bar via electronic discharge machining. In a last manufacturing step, the notch section of the samples was polished to avoid fracture due to surface defects. The two principal directions (i.e. longitudinal (L) and transverse (T)) of the bar were tested to account for possible anisotropy. Samples of every possible combination of material (TA6V, U500), thickness (0.7, 7 mm) and direction (L, T) were tested three times.

The polishing of the surfaces induced a non-optimal adhesion of the speckle pattern paint. The modification of the standard painting method allowed circumventing this problem. The iterative application of black and white droplets with intermediate pauses promoted a better adhesion of the paint. A resulting disadvantage was a heterogeneous pattern quality.

Loading-unloading tests were limited, because this kind of loading led to increased premature paint cracking. Further studies are necessary to clarify the best choice between surface condition, strain gradients and paint particle sizes. is 8.8 x 6.6 mm. The used telecentric lenses have superior optical properties in order to limit the effect of out-of-plane motion. With this optical setup, the physical size of one pixels is 6.4 𝜇𝑚 in single magnification and the standard displacement resolution is .2𝜇𝑚 per pixel. The corresponding strain resolution is .05 %. It was obtained by correlating two images at the very beginning of the test.

The observation from all four sides requires correct lighting of each surface with a separate spotlight. Two cameras were equipped with adjustable telecentric lens to allow more freedom in the positioning of the light sources. The observation setup was complemented by a fifth camera, which records the notch zone with a double magnification lens, cf. figure 3.7. The exploitation of this camera's data is not part of this study.

The cameras captured the scene at a frequency of one image per second. This restricts the amount of data stored, but ensures a sufficient amount of detail during the plasticity phase. Unfortunately, a DIC analysis of rupture requires an acquisition frequency only possible with high-speed cameras.

A priori assessment of 2D DIC An a priori assessment of the influence of out-of-plane motion was performed to check that the use of 2D DIC delivers accurate results. Supposing that the local maximum longitudinal plastic strain is less than 0.5, the incompressible out of plane strain is of the order of -0.25. For thin samples, the out-of plane displacement results to 0.7 𝑚𝑚/2 * 0.25 = 0.0875 𝑚𝑚. This displacement can lead to a magnification change in the image. The frontal distance of the lens is of the order of 100 mm. From Newton's lens formula, it follows that the magnification change leads to a maximum spurious dilatational strain of less than 0.1 % [START_REF] Fayolle | Contrôler une machine d'essai avec une caméra[END_REF]. For the thick sample one obtains analogously 1 %. For telecentric lenses, as used herein, these effects are significantly reduced in this displacement range. It is therefore expected that the effect of out-of-plane motion is negligible in the results reported hereafter.

The sample and the load cell of the testing machine may not be 100 % parallel to each other. In addition, the optical axis and the surface of the sample may not be perpendicular. This misalignment can be static or vary during the test.

Only the alignment about the out-of-plane axis, z-axis, is controllable in 2D DIC. Rotations around the x-axis, horizontal axis in the image plane, are visible from the perpendicular cameras. The alignment angle about the y-axis, vertical axis in the image plane, requires the analysis of several cameras together.

Definition of the region of interest and image selection

The region of interest (ROI) was defined to comprise the entire notch section for the thin samples. For the thick samples, the selected zone was reduced to avoid sharp boundary angles, which prevented the identification of DIC displacement boundary conditions. For the cameras facing the flat side of the samples, examples of both types are presented in figure 3.8. For the side cameras, the vertical extent of the ROI was chosen to closely correspond to the face camera ROI. The width of the ROI corresponds to the sample thickness.

One of the main factors that influence the time necessary to perform an I-DIC identification of plasticity is the number of non-linear FE steps. To limit the amount of necessary steps, the images were pre-selected based on force and displacement values as measured by the testing machine, cf. figure 3.9. It was considered that ten images selected in the linear range guarantee sufficient information in a state of pure elasticity. Another ten images from the end of the test were included, to have access to enough data with high strains. A uniform interval amongst the remaining images complemented the series, resulting in approximately 50 images per analysis. It was verified that the selected sequence and the entire image series led to the same residuals if DIC is used to identify the displacement fields. I-DIC and DIC are always performed on this sub-selection. A discussion of variation of the identification, if less than the maximum amount of data is used is presented in appendix E.

DIC parameters and mesh resolution

The correlation parameters of the DIC multi-scale algorithm were set to use three scales and a maximum of 500 iterations per scale. The correlation was performed incremental ny, i.e. the results of the previous image were used to find the displacement increments 𝛿𝑢.

The convergence limit was set to ||𝛿𝑢|| = 10 -4 pixel. With these restrictions, the minimum element size, which allows a complete measurement of the displacement fields on the entire test, was identified from preliminary tests, on several of the recorded speckle patterns.

The characteristic element size on the boundary of the ROI and for the entire DIC mesh was then uniformly set to 60 pixels. Within the geometry boundaries, meshes were automatically created using the FE-code ABAQUS. Meshes were constructed with quadrilateral elements, which were then divided to produce triangular elements. Figure 3.10 shows the meshes for the two presented ROI. For the thin sample, the number of degrees of freedom (DOF) is 376 for 314 elements. For the thick geometry one counts 216 DOF and 174 elements. For the two refined meshes, the central section is specified to have a characteristic element size of 25 and 10 pixels, respectively.

For the tested data, triangular (T3) DIC showed superior convergence capabilities than the quadrilateral (Q4) counterpart. Previous I-DIC applications were limited to low plastic strains. Mesh distortion was limited and triangular meshes were used throughout these studies. However, linear triangular elements lead to a pressure incompatibility at large plastic strain; quadrilateral elements are necessary. To combine the advantages of both Q4 and T3 elements, the FE analysis was performed with linear Q4 elements. These were divided into two triangles for the correlation of the images. Figure 3.11 shows examples of the quadrilateral meshes. The number of elements are 157, 549 and 2660 in a planar 2D mesh.

Furthermore, to prevent volumetric locking, reduced integration elements with hourglass stabilization were used as recommended in the ABAQUS user manual.

Dimensionality of I-DIC simulations

The I-DIC simulations were carried out with 2D plane stress conditions and in 3D using ABAQUS/STANDARD. For 3D FEM models, the meshes were extruded to 1, 4 and 6 layers of brick elements, respectively, leading to 157, 2196 and 15960 elements. The number of layers was chosen to ensure a good aspect ratio of the elements.

The thickness of the sample was measured mechanically before the test. In addition, the thickness was controlled with the recorded images.

Global parameters for I-DIC

The I-DIC sensitivity fields were obtained numerically. For their calculation, the variation of the model parameters was 1 %. To speed up the identification process, the initial parameter values corresponded to the best know fit, e.g. literature data, of similar sample or data from a coarser mesh. To limit possible divergence in the identification procedure the parameter increment was limited to 20 % from one iteration to the next. The algorithm terminated if 50 iterations were exceeded or if the norm of the total parameter variation was ||𝛿𝑝|| < 0.001.

The measurement noise, which is used for the dimensionless residuals, depends on the quality of the speckle pattern. A typical value is 𝛾 𝑓 = 230. The load cell noise was measured as 𝛾 𝐹 = 2.2𝑁 .

Definition of the boundary conditions

In the used I-DIC analyses, boundary conditions were prescribed. The Dirichlet boundary conditions of I-DIC analysis are identical to DIC displacements on the ROI boundaries. An error in DIC measured displacement of these boundary nodes is therefore inherent to all analyses reported hereafter.

The combination of the measurements of several cameras requires the relative positioning of the cameras in a global coordinate system. This procedure requires further studies for the presented 2D camera approach. Therefore, the presented analyses concern only the case of boundary conditions homogeneous in thickness.

For 3D models, the extrusion of the DIC boundary displacements ensured this condition. The use of either camera facing the two flat surfaces is possible. This procedure resulted in two different analyses for every sample.

Given that the speckle-pattern on both sides should have an equal quality, it is possible to comment on possible misalignment of the sample by comparing these two analyses. It is therefore assumed that a problem of alignment exists if the identified models vary too much from one side to the other.

Parameters of the selected material models

The combined inverse identification of elasticity and plasticity proved being difficult in previous studies [START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF]. In addition, Poisson's ratio, which determines the elastic kinematics, does not vary much for metals. Therefore, the reference values were chosen and kept constant in the entire study: 𝜈 = 0.32 for Ta6V and 𝜈 = 0.3 for U500. Elastic modulus and plastic law were identified together.

The selection of the material models is made mainly among those available in the ABAQUS/STANDARD software, cf. table 3.1. The three hardening laws Ludwik, Voce and extended Voce are studied using ABAQUS parameterization option for the yield curves.

To analyze the effect of the stress state, the criteria of Drucker-Prager, Hill and Hosford were analyzed for the extended Voce hardening. The studies were performed on the intermediate 3D mesh (I-DIC-i). The initial parameters were set to the ones obtained from the analysis with the von Mises criterion. For Hill's anisotropic yield criterion, yield ratios are defined; 𝑅 𝑖𝑗 , is defined as the ratio between yield stress in the direction of the stress component 𝜎 𝑖𝑗 and for a reference yield stress 𝜎 0 (𝑝),

𝑅 𝑖𝑗 = 𝜎 𝑖𝑗 𝜎 0 (𝑝) . (3.2)
The relation of yield ratios and the criterion coefficients is detailed in appendix H. The longitudinal sample x-axis coincides with the principal material direction 𝑅 11 , cf. figure 3.12; this direction is set to be constant, 𝑅 11 = 1. The axis of primary loading is the y-axis which is associated with 𝑅 22 . An ABAQUS UMAT routine allowed for the application of the Hosford criterion. However, to mitigate convergence issues, the optimization was not performed on the model parameter 𝑛, cf. equation (2.19). Instead, the exponent of the criterion was prescribed and the best fit for the hardening model was obtained. Nine different exponents were tested, 𝑛 = [4, 8, 16, 24, 34, 44, 64, 74, 84]. The comparison of the total residuals indicated the best fitting exponent 𝑛. 3.3 Results of samples made of TA6V

R1

Digital image correlation

In the following sections, a thin sample serves for the explanation of various characteristics of the measurement and identification techniques used. Many of the reported results are also valid for the other tested samples.

Definition of measurement scales

The DIC method has been frequently used for direct identification of the stress-strain law. The following demonstration highlights the drawback of such an identification. The force-displacement curve, as measured by the testing machine, is displayed in figure 3.13. The influence of the testing machine leads to a very long linear phase of almost half the total displacement. A descent of the curve can be observed. where, the reference length 𝑙 0 is the vertical extent of the extensometer. A macroscopic length was defined as the boundaries of the ROI, 𝑙 0𝑚𝑎𝑐𝑟𝑜 = 3.97 𝑚𝑚, (A) in figure 3.14(a). The mean longitudinal boundary displacement ū𝑦 is obtained from the DIC displacements.

A second scale is defined as a small ligament, approximately equal to the size of one DIC element. It is called mesoscopic scale, 𝑙 0𝑚𝑒𝑠𝑜 = 0.23 𝑚𝑚, (B). Due to the non-structured character of the mesh, the average axial strain of the elements in this band is considered.

The main challenge lies in the definition of the stress. As explained, the transformation of the engineering stress to true stress assumes a uniform field. Due to the strong variation of strains, this hypothesis is not valid. In the following a global measure of the stress will be defined as the ratio of the applied force to the initial section of the ligament An indication of the strain concentration in the specimen is clearly visible from the differences between the macroscopic strain measurement, 𝜀 𝑚𝑎𝑐𝑟𝑜𝑚𝑎𝑥 = 0.07, and the mesoscopic one, 𝜀 𝑚𝑒𝑠𝑜𝑚𝑎𝑥 = 0.28. In what follows the macroscopic strain 𝜀 𝑚𝑎𝑐𝑟𝑜 will be used as loading parameter.

𝑆 0 = 𝑡𝑤 = 0.72 * 3.11𝑚𝑚 2 = 2.24𝑚𝑚 2 , cf.
Finally, the actual local maximum strain, as measured by DIC, is above 𝜀 𝑦𝑚𝑎𝑥 = 0.35, which indicates that a measurement at the meso scale does not lead to a precise estimation of the strain level within the ligament.

Displacement and strain fields via DIC Figure 3.15 displays the longitudinal displacement and strain fields at three different instants, namely, the instant of maximum force in image 28, 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.034, an instant of high strains in image 37, 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.061, and the final instant, just before fracture, in image 47, 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.072. As can be seen from figure 3.15 the concentration of the strain increases during the test due to plasticity.

The DIC data in the deformed mesh are projected on the associated image. The fields fit quite accurately the recorded shape. The displacement fields are smooth and show a strong variation in the notch region. Due to the mesh coarseness, the strain field interpretation is not obvious. For the instant of the maximum force, the strain is still relatively homogeneous. The deformation at the notch tip is slightly higher than in the middle of the notch section. As deformation progresses, the strain increases and the final measurement exceeds 0.35 %. The higher strains are located on the right side and in the middle of the notch. However, the very strong gradients between the elements indicate the occurrence of speckle pattern defects. Figure 3.16 shows additional displacement fields from this test. The asymmetry of the transverse displacements indicates the unperfect boundary conditions. Even though the geometry is not very thick, the sides of the sample were observed and an attempt of the measurement of the displacement fields was made. Due to the low number of elements and the curvature of the notch the accuracy is limited. These fields indicate that necking in thickness took place at the final instant at the smallest cross section. From the recorded transverse displacement, it is possible to assess the influence of out-of-plane displacement.

For the presented example, the measured transverse displacement remains less than 50 pixels or 320 𝜇𝑚. The frontal of the face cameras objective is 100 mm. These values introduce a maximum spurious strain of 0.32 % for a conventional lens.

The variation of the longitudinal displacement field remains very low in the thickness direction and constant on the boundary. Therefore, the extrusion of the DIC boundary displacements for the I-DIC method seems to be acceptable.

Measuring with residuals

The actual gray level balance, which is obtained in the displacement identification, can be controlled through the examination of the residual image. For the three chosen instants, the normalized residuals 𝜙 were calculated, cf. figure 3.17. The dynamic range of the initial image was used for the normalization, 𝛥𝑓 = 41423 gray levels. High localized residuals are observed in the ligament region (dark spots) at the instant of the maximum force, image 28, 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.034. These residuals increase further for 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.061. However, an additional increase for 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.072 remains limited.

The residuals in the top half of the sample are low, despite the high levels of displacements. In contrast, the overall increase in the residuals in the notch region is also an indication that the chosen kinematics, associated with the coarse FE mesh, is not capable to fully capture the actual displacement fluctuations. A more comprehensive overview on the quality can be gained from the time-dependent dimensionless residuals, 𝜙 𝑐 (𝑡), shown in figure 3.18. The three chosen instants are indicated with vertical dash-dotted lines. After being almost constant, the residuals start to deteriorate after 20 images. However, this deterioration process reduces towards the end of the test. A fact also reflected by the residual images that is similar for image 37 and 47.

For camera 1, the total mean residual is φ𝑐 = 0.028. This indicates that the measured displacements are overall trustworthy; the displacement field fits reasonably to the image.

General assessment of integrated digital image correlation

The test data were analyzed with integrated digital image correlation. In the following paragraphs, some characteristic results are presented.

Decreasing the measurement scale Figure 3.19(a) presents the measured force and the resulting forces corresponding to the identified model (Ludwik hardening model, I-DIC-3D-f mesh). The loading parameter 𝜀 𝑚𝑎𝑐𝑟𝑜 is used for the abscissa. The three instants are marked again with vertical lines. The overall estimation of the forces is globally satisfactory. This means that the macroscopic softening can be described accurately by making use of a plastic model.

The reported literature values for the maximum equivalent plastic strain remain below 16 %, cf. section 2.1. This I-DIC analysis identifies a maximum value of 𝑝 𝑚𝑎𝑥 = 0.37, cf. figure 3.19(b). At the maximum force, the equivalent plastic strain reaches values up to 14 %. For the second depicted instant it rises up to 27 %.

Assessment of mechanical fields

Three dimensional simulations allow the visualization of various fields on the monitored surface and within the bulk. In the following discussions, the surface and the midsection plane (parallel to the observed surface at half the sample thickness) are considered. Similarly, to the longitudinal strain, equivalent plastic strain is localized within the notch vicinity. Even though the overall levels are close between surface and midsection planes, differences exist. In particular, the field in the midsection plane is more localized than on the external surface. Moreover, the maximum level of equivalent plastic strain is higher (𝑝 = 0.37 vs. 𝑝 = 0.29).

The associated von Mises equivalent stress field 𝜎 𝑒𝑞 (𝑥,𝑡) is displayed in figure 3.22 in the midsection plane for two instants in the loading history. The first one corresponds to an elastic step, 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.003. The second corresponds to the final instant just before the final fracture, 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.072. In elasticity, stress concentrations occur at the notch root.

In contrast, at the end of the test, the entire notch section has an almost uniform stress level. The triaxiality 𝜂(𝑥,𝑡) is another field that can be extracted from the FE simulations. In the elastic regime, the field is very similar in terms of its distribution in both observed planes, cf. figure 3.23. The maximum levels are around 𝜂 = 0.5. In contrast, in the final step, the triaxiality field is more localized in the midsection plane than on the external surface, cf. figure 3.24. Further, the maximum levels are no longer similar (i.e. 𝜂 = 0.7 vs. 𝜂 = 0.59). To study this evolution, the triaxiality of four characteristic points was plotted as a function of the corresponding cumulated plastic strain, cf. figure 3.25. Two of the considered points are located on the notch tip and the other two in the center of the two analyzed sections. The two points close to the notch tip remain virtually at a constant stress triaxiality, around 𝜂 ≈ 0.34. The free-edge condition is the primary reason, as it limits the maximum achievable triaxiality. Conversely, the two central points experience significant fluctuations of stress triaxiality during the test. The point located in the midsection plane reaches the highest levels (i.e. 𝜂 = 0.63 close to the total maximum in the sample of 𝜂 = 0.7) at the end of the experiment.

Temporal evaluation of residuals

As for DIC, the registration quality is assessed with the residuals. The residual images at the three marked instants are depicted in figure 3. 26(a-c). To get a more quantitative comparison between DIC and I-DIC, figure 3. 26(d-f) reports the absolute difference of correlation residuals between the two analyses, 𝛥𝜙 = 𝜙 𝐼-𝐷𝐼𝐶 -𝜙 𝐷𝐼𝐶 . Except in the central part of the sample, the levels are very small and close to noise contributions. While the differences in the residual images between 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.061 and 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.072 are not distinguishable, the total difference 𝛥𝜙 for 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.072 in figure 3.26(f) shows higher values.

The total mean residual for I-DIC is φ𝑐 = 0.033, which also indicates trustworthy measured displacements. This level is slightly higher than the one observed with DIC ( φ𝑐 = 0.028). Such a result is to be expected in case a model error occurs. Here, model error refers to the choice of the constitutive law but also to the approximation of the geometry, to the boundary conditions as well as to the discretization errors. Nevertheless, the residual with I-DIC is close to that with DIC, this indicates that those errors are limited.

Figure 3.27(a) shows the evolution of the mean dimensionless correlation residuals 𝜙 𝑐 (𝑡) for the entire image sequence. The I-DIC analysis follows the trend of the DIC residuals. For 𝜀 𝑚𝑎𝑐𝑟𝑜 < 0.02, the dimensionless residuals are close to the acquisition noise. Both approaches are almost identical in terms of residuals. In this part of the experiment, the chosen model is fully validated according to the correlation residuals. In contrast, for 𝜀 𝑚𝑎𝑐𝑟𝑜 ≥ 0.2 there is a clear degradation of both residuals which corresponds to a more or less progressive degradation of the speckle pattern.

The total residual consists of the contribution of the equilibrium gap and the correlation residuals. All expressions are relative to the acquisition noises, which for the tested sample are 𝛾 𝑓 = 224 gray levels, 𝛾 𝐹 = 2,2𝑁 , cf. figure 3.27(b). For the dimensionless force residual 𝜒 𝐹 (𝑡), most of the difference occurs in the elastic regime. The total value is 𝜒 𝐹 = 11. The force residuals rise towards fracture at 𝜀 𝑚𝑎𝑐𝑟𝑜 > 0.069. The total identification residual is 𝜒 tot = 9.26. This level is one order of magnitude greater than a level expected if only noise contributed to the residuals. This confirms the hypothesis of a model error. In a simpler situation (i.e. strain localization did not occur in a dog-bone sample), levels of the order of 5 have been reported [START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF]. It is believed that the model even though not perfect is a good approximation of the experimental observation.

Comparison of 2D and 3D analyses and influence of the element size Mesh size effects as well as dimensionality (2D vs. 3D simulations) were analyzed. Table 3.2 presents an overview of possible gains. There is a clear gain switching from 2D plane stress simulations to 3D simulations, for any type of discretization. Consequently, even though the thickness of the sample is very small (i.e. 0.72 mm), a plane stress computation induces a model error that is detected.

The observed variation of the mechanical fields in the thickness emphasizes the interest of 3D simulations. Of course mesh sizes have to be adapted differently for 2D and 3D analyses.

In the plane stress case, the intermediate mesh provided a lower identification residual, but the maximum strain almost doubled from 𝑝 𝑚𝑎𝑥 = 0.39 to 𝑝 𝑚𝑎𝑥 = 0.67. The finest

Approach

𝜒 tot 𝑝 max I-DIC-2D-c 10.68 0.39 I-DIC-2D-i 10.51 0.67 I-DIC-2D-f 11.43 0.92 I-DIC-3D-c 9.58 0.29 I-DIC-3D-i 9.26 0.36 I-DIC-3D-f 9.26 0.37 mesh led to significantly higher residuals and higher strain localizations. The maximum equivalent plastic strain was 𝑝 𝑚𝑎𝑥 = 0.92. In contrast, for 3D simulations, the two finer discretization led to an almost constant maximum equivalent plastic strain and equal identification residuals. This result indicates that numerical convergence of the FE simulation is reached. In addition, there are almost no further gains from I-DIC-3D-i to I-DIC-3D-f. Figure 3.28 presents the hardening curves, which are associated with these models. The yield stress is relatively constant for all identifications. 2D identifications led to significantly higher stresses.

For 2D simulations, the refinement resulted in a rise and "linearization" of the curve. The mesh refinement leads to higher equivalent plastic strains in the analysis. This higher plasticity then induces a variation in the hardening curve, which in turn influences the equivalent plastic strain. In contrast, for 3D calculations, the shape of the curve remains almost constant for fine and intermediate meshes, only the stress level for the coarse mesh is lower.

Control of the sample alignment

Assuming an equal quality of the speckle pattern on either sample side, optimal alignment of the sample would lead to the same identification result independently of the sample surface that is used for I-DIC analyes. Figure 3.29 presents two different samples. Good alignment resulted in nearly identical displacement fields and thus hardening curves. In contrast, a test with misaligned components led to more important differences in the identified hardening curves. Allowing more freedom in the hardening law does not induce a significant improvement of the identification quality. The variation is of second order in the global identification residual 𝜒 tot , cf. table 3.3. One may state that the overall model error is not due to a too poor hardening model. When an elastic law is assumed during the whole experiment, a very significant degradation of the global identification residual 𝜒 tot is observed. Figure 3.30 shows the equivalent stress/equivalent plastic strain responses for these hardening models. The predicted equivalent stress levels are very similar below 𝑝 = 0.2. In contrast, at 𝑝 = 0.37 the curves predict the values 𝑅(𝑝 = 0.37) 𝐿𝑢𝑑𝑤𝑖𝑘 = 1253 𝑀 𝑃 𝑎, 𝑅(𝑝 = 0.37) 𝑉 𝑜𝑐𝑒 = 1220 𝑀 𝑃 𝑎, 𝑅(𝑝 = 0.37) 𝑒𝑥𝑡.𝑉 𝑜𝑐𝑒 = 1249 𝑀 𝑃 𝑎. The hardening law with the higher number of parameters lies between the other two postulates. These differences are enough to invoke a different equivalent plastic strain at fracture, cf. table 3.3.

Construction of average hardening curves

The extended Voce law was identified on all thin samples. As the identified hardening curves are not identical, an average curve was calculated from the test data. This is illustrated for the longitudinal direction in figure 3.31(a).

For each sample, both Camera 1 (C1) and Camera 4 (C4) were analyzed. Due to the non-linearity of the laws, the mean fit curve is calculated from the six curves and not from the coefficients of the curves. The maximum equivalent plastic strain is set to the average of all six analyses.

The average material behavior is identical for both directions, cf. figure 3.31(b). At 𝑝 = 0.1, the predicted hardening stresses are R(𝑝 = 0.1) 𝐿 = 1112 𝑀 𝑃 𝑎, R(𝑝 = 0.1) 𝑇 = 1132 𝑀 𝑃 𝑎. The relative differences between the directions were of the same order as the variation for samples of the same direction. Similar observations were made for the maximum equivalent plastic strain, p𝑚𝑎𝑥 𝐿 = 0.39, p𝑚𝑎𝑥 𝑇 = 0.38.

Identification of yield criteria

To analyze the effect of the stress state, the criteria of Drucker-Prager, Hill and Hosford were analyzed for extended Voce hardening on a sample of longitudinal orientation. The intermediate mesh was used as it guarantees already a decrease in the residuals and allows a faster computation. The minimization algorithm was not able to identify the Drucker-Prager criterion. The residuals were always higher than the von Mises model and the model parameter 𝛽 was close to zero, which suppresses any influence of the triaxiality in the criterion.

The application of the Hosford and Hill criterion led to a reduction of the total residuals compared to the identification with the von Mises criterion, cf. table 3.4. As expectable from the macroscopic results, the application of the Hill criterion to TA6V only led to a gain of 0.8 %.

R11 R22 R33 R12 R13 R23

1.00 0.87 0.91 0.82 0.50 0.49

The hardening curves for both materials are presented in figure 3.32(a) and the coefficients of the Hill criterion in table 3.5. The shape of the hardening curve remains similar; only the stress level changes significantly: 𝑅(𝑝 = 0.1) 𝑀 𝑖𝑠𝑒𝑠 = 1146 𝑀 𝑃 𝑎, 𝑅(𝑝 = 0.1) 𝐻𝑖𝑙𝑙 = 1346 𝑀 𝑃 𝑎.

For the Hill criterion, the transverse direction was set constant, 𝑅 11 = 1. The longitudinal coefficient is 𝑅 22 = 0.87. However, the macroscopic analysis of plasticity does not show any sign of anisotropy of this order. One reason could be that the displacement field is not heterogeneous enough to analyze precisely a complete anisotropic law by one uniaxial test. While it remains difficult to compare the curve of the anisotropic model with the isotropic ones, one can observe a shift of the associated hardening curve.

For the Hosford criterion, the exponent 𝑛 was prescribed. The total residuals were normalized with respect to the level when at 𝑛 = 4, cf. figure 3.32(b). After an initial decrease the change remains almost independent of 𝑛. As 𝑛 approaches 80 the residuals increase again. The minimum value was obtained for 𝑛 = 64. The identification of the hardening law led to an increase in stresses, 𝑅(𝑝 = 0.1) 𝐻𝑜𝑠𝑓 𝑜𝑟𝑑 = 1236 𝑀 𝑃 𝑎. These results are similar to the findings reported in literature [START_REF] Mazière | Eclatement des disques de Turbomachines -Burst of turboengine disks[END_REF]. However, compared to the von Mises identification, the obtainable residual gain is only 5.6 %. In addition, similar gains are already obtained for any 𝑛 > 16. The high residuals have most likely another cause.

The use of the Hosford or the Hill criterion raised the stress level. It is important to note that this led to a decrease of the fracture strain, 𝑝 𝑚𝑎𝑥 𝑀 𝑖𝑠𝑒𝑠 = 0.36, 𝑝 𝑚𝑎𝑥 𝐻𝑜𝑠𝑓 𝑜𝑟𝑑 = 0.35, 𝑝 𝑚𝑎𝑥 𝐻𝑖𝑙𝑙 = 0.29.

Influence of the thickness on the identification DIC displacement fields for thick samples

The characteristic displacement fields acquired for a thick sample are presented in figure 3.33. The DIC fields on the flat surface are similar to those of the thin sample. Of primary interest in this case are the fields of the side camera. Due to the larger thickness, DIC of the side camera images results in a higher quality of the displacement fields.

The side cameras are not used for I-DIC, instead they serve two other purposes: the verification of the negligible out-of-plane displacement and of the boundary condition extrusion for I-DIC.

The preliminary assessment of out-of-plane strain of the thick sample resulted in a 1 % spurious strain. For the presented example, the measured transverse displacement remains less than 25 pixels or 160 𝜇𝑚. The frontal of the face camera objective lens is 100 mm. These values introduce an approximate spurious strain of 0.16 %, which is very low.

The second verification concerns the extrusion of the boundary condition for 3D I-DIC models. In fact, the horizontal displacement is slightly inclined through the thickness. It is higher on the left side than on the right side. Future studies on thick samples should consider this difference to evaluate the resulting residual error. No strong variation can be found on the boundary in figure 3.33(b).

Identification of plastic hardening law for thick samples

The analysis of the thick samples was restricted to the identification of the plastic hardening behavior with the von Mises criterion. The most important point, which will be discussed in the next chapter, is the difference in the failure strain. The mean value was p𝑚𝑎𝑥 𝐿 = 0.17, cf. appendix F. The residuals of the thick samples were on average a lot higher. For the example presented in the following, the value was 𝜒 tot = 24.84. This means that the results as such are less trustworthy than for the thin samples. However, the calculated forces from the identification are still close to the experimental forces, cf. figure 3.34. This emphasizes the high sensitivity of the residuals.

To assess the difference between thin and thick samples, the average hardening law of thin samples was verified on the data of a thick sample. For the chosen example the verification residual rose to 𝜒 tot = 159.97. There is a clear degradation of the overall quality. In addition, the forces are significantly lower. We will discuss the assumed origin of this difference after the presentation of the results for U500.

I-DIC results of samples made of U500

Identification of plastic hardening

A thin sample of Udimet 500 was tested with the Ludwik and the extended Voce law. Table 3.6 shows that the differences in the total residuals and in the maximum equivalent plastic strain are of similar order as for TA6V and that the extended Voce law produces the lower residuals. Therefore, all thin samples were examined with the extended Voce law and an I-DIC mesh of the finest resolution. longitudinal direction, the average maximum equivalent plastic strain is p𝑚𝑎𝑥 𝐿 = 0.31, whereas for the transverse direction only an average of p𝑚𝑎𝑥 𝑇 = 0.12 was obtained. This is about 39 % of the longitudinal value. In addition, the levels of stress in the transverse direction are significantly lower, R(𝑝 = 0.1) 𝐿 = 1243 𝑀 𝑃 𝑎, R(𝑝 = 0.1) 𝑇 = 1019 𝑀 𝑃 𝑎. This indicates that the transverse direction could be as much as 18 % weaker.

Influence of the yield criterion

The identification was carried out with the intermediate mesh. In contrast to TA6V, the residuals were minimally different from the identification with the finest mesh for the von Mises criterion. The maximum equivalent plastic strain reduced to 𝑝 𝑚𝑎𝑥 = 0.28. As the test samples showed a larger macroscopic anisotropy, the identification of the Hill criterion lead to a larger reduction of the total residuals compared to the identification with the von Mises criterion, cf. table 3.7. The gain is 7.6 %, while for the Hosford criterion, it is 3.1 %, which is less than for TA6V.

The hardening curves are presented in figure 3.36 and the coefficients of the Hill criterion in table 3.8. As for TA6V, the direct comparison of anisotropic and isotropic model parameters is delicate. For the Hill criterion, the curve is lower than the isotropic model curves and the maximum equivalent plastic strain increases. This behavior corresponds to the observations made for TA6V.

For the Hosford criterion, the total residual remains almost independent of 𝑛, after an initial decrease. The overall minimum is identified at 𝑛 = 74. The hardening curve is raised upwards 𝑅(𝑝 = 0.1) 𝐻𝑜𝑠𝑓 𝑜𝑟𝑑 = 1315 𝑀 𝑃 𝑎. In addition, the maximum equivalent plastic strain changes from 𝑝 = 0.28 to 𝑝 = 0.29, which contradicts the previous observations.

Assessment of mechanical fields

For Udimet 500, the identified hardening curves led to higher stresses and have a higher hardening modulus than for TA6V. In the following, it is examined how this behavior influences the mechanical fields in the sample. Figure 3.37 shows selected fields from the midsection of the longitudinal sample in figure 3.35(a) at two different instants.

The first instant is when 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.029. Plasticity is established and high values are concentrated at the notches. However, the plastic field extends almost through the entire

R11 R22 R33 R12 R13 R23

1.00 1.13 1.02 1.25 0.13 0.96 As for TA6V, four characteristic points were analyzed, cf. figure 3.38. The notch tip points have a higher triaxiality in the intermediate strain state than for TA6V. The points in the center are initially identical to the TA6V points. In addition, little difference exists for the point on the surface and the point in the midsection plane. However, all points undergo an evolution to lower triaxiality, 𝜂 ≈ 0.38 for the notch tip, 𝜂 = 0.47 for the central points.

From these findings one can draw the following conclusion: steep hardening curves lead to further strain concentration and, in the used samples, to an evolution towards lower triaxiality values. In contrast, a material with a small hardening slope leads to a de-localization of plasticity in the notch section. In the used samples, the triaxiality evolves to higher values.

Influence of the thickness on the identification

The same procedure as for TA6V was applied to identify the difference between thin and thick samples. As for TA6V, the increase in thickness led to a lower maximum equivalent plastic strain, p𝑚𝑎𝑥 𝐿 = 0.18. For the thick sample used in the following, the total identification residual was 𝜒 tot = 26.77. The verification residual rose to 𝜒 tot = 451.23, a value significantly higher. There is a significant effect between the two sets of material parameters. A circumstance also visible on the force -elongation curve, cf. figure 3.39.

Discussion of the results of the thick samples

The difference between identification and verification forces is higher for U500 than for TA6V. To examine the possible origin of such a difference, the hypothesis of the extraction of the boundary conditions was verified on a FE model of the entire sample geometry. The sample was modeled once completely and once only as section of 4 𝑚𝑚 length, similar to the extent of the used ROI in I-DIC, cf. figure 3.40(a). The nominal geometry and the identified mean constitutive model parameters of the thin samples were used. Mesh constraints were applied to ensure the mesh of the ROI corresponds to the mesh in the same region of the complete model. In the case of the complete model, the sample was loaded uniformly under tension at its upper end, marked with (1), and fixed at the lower end, point (2). This loading leads to similar load conditions on the ROI boundary as in the experiment and is therefore denominated as "true" boundary condition (BC). After the calculation of the model the longitudinal displacement at the limits of the ROI were extracted on the surface, point (3), and a 4th order polynomial was fitted to these data. The obtained function was then applied as load condition to the boundary surface nodes of the ROI model and extruded for all nodes in thickness direction. A procedure that is assumed to be almost equivalent to the processing of the experimental data in I-DIC.

The mean displacement of the ROI boundary and the associated macroscopic elongation was calculated for both models (true and extruded BC). Figure 3.40(b) shows the forceelongation curves. For the ROI model with the extruded boundary conditions the forces are lower. As in the verification, the differences between true and extruded boundary conditions are larger for U500 compared to TA6V.

This leads to the following explanation. In the I-DIC identification, the lower forces where compensated by the stress level of the material model, which was higher for the thick samples compared to the thin samples.

The difference in the resulting errors between TA6V and U500 can be explained as follows. U500 has a lower yield stress and a stronger hardening modulus. This results in a greater extent of the plasticity field as observed for the thin sample in figure 3.37. In addition, the plastic strain field is not homogenous in thickness. For a similar size of the ROI, the experimental displacement boundary surface is then more deformed in thickness compared to TA6V and the assumption of the extrusion is less correct. To improve accuracy, more elaborate procedures are necessary to deal with these kind of boundary conditions. One solution is the enlargement of the ROI. However, this was difficult as the strong notch curvature and the resulting element form made DIC impossible. An improvement could be achieved through the interpolation from separate coarse DIC solutions to fine boundary distributions for I-DIC. Another solution would be the introduction of the boundary conditions to the unknowns. In this case DIC would not be necessary and no coarse meshes need to be used.

Summary and conclusion from I-DIC analyses

This chapter presents the design of the experiments for tensile tests at high triaxiality values and the analyses of the tests with integrated digital image correlation. The triaxiality state of the disk was analyzed. The range of interest is between 𝜂 = 0.4 and 𝜂 = 0.75. Two different notched tensile samples were designed for the two materials TA6V and Udimet 500: a thin sample with a design range of triaxiality of 0.33 < 𝜂 < 0.5 and a thick sample with a design range of 0.33 < 𝜂 < 0.8. Gray level images were recorded during the test.

The application of the digital image correlation methods to the tests has notably revealed the difference of the maximum strain for thin and thick samples, which is further discussed in the next chapter. The importance of using 3D models even for thin samples must be underlined. These models allow for a better identification as can be seen from the associated residuals. The advantage of I-DIC over other FEMU identifications was shown with the decrease in the residuals with I-DIC allowing for the use of finer meshes.

Considering the two materials, the tested samples were first analyzed using a von Mises criterion. Then the Hill and Hosford criteria were identified on a thin sample. The study of TA6V revealed that the equivalent plastic strain reaches very high values (almost 40 %). For the thick samples, the maximum equivalent plastic strain reaches roughly only half of the values reached on thin samples. Thin and thick samples are characterized by different values of triaxiality. Therefore the influence of this parameter on the failure criteria will be studied in the next chapter. For Udimet 500 similar results were observed.

In addition, the study has highlighted where improvements should be researched to increase the accuracy and reliability of the method. This concerns first the application of speckle patterns to zones of high strain gradients. Second, for the thick samples the boundary conditions need to be improved through separate coarser DIC analysis of these regions or taking advantage of I-DIC and including the boundary conditions to the identification.

For the analysis of the material models, another improvement would be the creation of one large I-DIC analysis, which considers the entire recorded image data from all four cameras. This would allow improving the boundary conditions applied for the DIC and I-DIC analyses.

CHAPTER 4

Towards the identification of material failure criteria Up to now, the analysis focused on the determination of the constitutive law prior to failure. In this chapter, we seek possible failure criteria based on the identified constitutive law. First, a microscopic analysis of the failure surface and a macroscopic analysis are proposed. In addition, we tried an image-based analysis of the local state of stresses and strains in the vicinity of the failure surface at the instant just before failure. This analysis allowed the evaluation of different criteria to be performed with a large amount of data. Even though this approach is at its outset, it seems to us that it is a first step towards what could be a "big data" procedure for the identification of failure models.

Characterization of the fracture of TA6V

4.1.1 Microscopic analysis An important step in the analysis of the materials is the characterization of the fracture surface. The examination of the broken sample surface was performed with scanning electron microscopy. Figure 4.1 shows an overview on the fracture surface and a zoom.

The presented examples are representative of thin and thick samples. The fracture surface is not flat and resembles the cup cone structure, which occurs in ductile materials. On a closer inspection of the surface, two distinguishable features are revealed. Considerably larger voids exist on the surface, i.e. at point (1). The size of the holes reaches almost 25𝜇𝑚. A possible explanation is the decohesion of grains or inclusion particles. In addition, repetitive structures cover the entire surface at point (2). These features are dimples typically observed in ductile rupture. Therefore, the fracture may be described as predominantly transgranular.

These observations indicate that the final fracture mechanism of TA6V corresponds to the definition of ductile rupture.

Macroscopic analysis

Ductile fracture is the result of void nucleation and growth, which are associated with the damage state variable. If substantial and diffusive amount of damage occurs in the tested sample, the measure of the global elastic stiffness allows the damage evolution to be identified. If such a method is commonly used for composite materials, such a global decrease in the apparent modulus of elasticity is not documented in the literature for metallic materials. Nevertheless, the DIC and especially I-DIC procedures have proved to be sensitive to small changes unseen with other techniques by using residuals.

To study possible occurrence of damage with the used test equipment, several unloading and reloading cycles were included for a thin sample. In a first step, the identification of the optimal elasto-plastic model parameters was performed for an image selection of continuously growing elongation, excluding the cycles. In a second step, the model parameters were applied to the entire test data, including the cycles. The entire forceelongation curve is well approximated including the macroscopic softening as well as the cycles, cf. shows the maximum equivalent plastic strain that was reached at every macroscopic elongation in I-DIC. The instants of the last two unloadings are marked in both diagrams with 𝜀 𝑚𝑎𝑐𝑟𝑜 𝐼 = 0.032 and 𝜀 𝑚𝑎𝑐𝑟𝑜 𝐼𝐼 = 0.053. At these instants, the equivalent plastic strain is 𝑝 𝐼 = 0.163 and 𝑝 𝐼𝐼 = 0.303. The absence of significant deviation of I-DIC from the experiment means that the global effect of a possible damage is not noticeable at 𝑝 𝐼𝐼 . Unfortunately, due to the strong heterogeneity of the test, this does not imply that no local damage exists in the specimen for these values of equivalent plastic strain. Additional cycles after 𝜀 𝑚𝑎𝑐𝑟𝑜 𝐼𝐼 could deliver more precise information. However, this would require a higher cycle frequency. In the performed test campaign, this was avoided due to the increased likelihood of paint cracking.

Due to the lower level of plastic strain in the thick samples compared to the thin ones, no visible global effect of possible damage was expected prior to failure. Therefore, no cyclic tests were performed on the thick samples.

Instead of analyzing the possible damage behavior, we focused on the identification of a failure criterion. It seems that the damage nucleation and growth phases, which could precede the TA6V fracture, can be neglected in a first instance, as it occurs very late and for a short duration. The fact that thin and thick samples fail for very different values of the equivalent plastic strain led to the investigation of various criteria, notably those involving the effect of triaxiality.

The principle followed for the identification of a failure criterion was first to determine the state of strains and stresses "on" the failure surface at an instant as close as possible to the time of failure. Then from these data, one can seek a possible criterion and check its consistency with all the information available in the computation; that is to check that the proposed failure criterion was not cracked at an earlier time and/or in another location.

Image-based analysis of the cracked surface

The image after fracture cannot be analyzed by the current I-DIC procedure, as the sample is discontinuous. However, it is possible to correlate one part of the broken sample with the preceding image, that is, the last image of the I-DIC analysis. The final displacement field is mostly limited to elastic unloading and a large rigid body motion. The rigid body motion results from the feedback response of the testing machine as the force drops to zero. This post-fracture displacement field could be identified by the restriction of a DIC analysis to one half of the sample.

However, to demonstrate the principle of the selection process, this is not necessary and such a DIC procedure was discarded. The elastic unloading was neglected and only the rigid body motion was identified. The last image was shifted, until the upper part of the broken sample corresponded globally to the location of the unbroken sample in the previous image, cf. figure 4.3(a). The data of the last I-DIC step were superimposed onto this shifted image of fracture, cf. figure 4.3(b). As the upper part was used for the manual correlation, the mesh and the recorded sample limits in the lower part do not correspond to each other.

In the upper part, the superposition allows for the selection of the elements, which are closest to the fracture surface. The selected elements in this critical zone are shown in figure 4.3(c). These selected elements were extracted from the mesh, cf. figure 4.3(d). The color indicates the equivalent plastic strain within the elements; a darker color symbolizes a higher value.

For the thick sample the same procedure was performed, cf. figure 4.4. The selected zone depends on the actual test conditions. The fracture surface can vary in the out-of-plane direction. However, for 3D simulations the boundary conditions were a simple extrusion of the face boundary conditions. Therefore, the same procedure was applied for the 3D fracture surface. All elements that are an extrusion of the surface selection were extracted from the 3D model.

Construction of the fracture locus Choice of quantities for the fracture locus

The quantities that are used for the prediction of fracture vary for existing fracture criteria.

The quantities are either stress, strain or energy-based. Amongst the stress based quantities are the different equivalent stresses, as for example the von Mises stress or the Tresca stress.

The maximum principal stress is frequently used to describe brittle fracture. Finally the stress invariants are used in more complex fracture criteria.

The strain based quantities are associated to these stresses. The most commonly employed quantity is the equivalent plastic strain, which is used for ductile materials. Other quantities are the principal strains. A particular case is the use of damage criteria. Damage can be driven by the state of stress and strain. In this aspect damage is similar to an energy based criterion. For example, the strain energy density for small strain theory is defined as

𝑊 = 1 2 3 ∑ 𝑖=1 3 ∑ 𝑗=1 𝜎 𝑖𝑗 𝜀 𝑒 𝑖𝑗 = 1 2 (𝜎 𝑥 𝜀 𝑒 𝑥 + 𝜎 𝑦 𝜀 𝑒 𝑦 + 𝜎 𝑧 𝜀 𝑒 𝑧 + 𝜏 𝑥𝑦 𝛾 𝑒 𝑥𝑦 + 𝜏 𝑦𝑧 𝛾 𝑒 𝑦𝑧 + 𝜏 𝑥𝑧 𝛾 𝑒 𝑥𝑧 ). (4.1)
Fracture locus based on literature findings From these quantities it is possible to create any possible combination for a fracture locus.

The most common representation describes the maximum equivalent plastic strain over a range of triaxiality, cf. section 2.7. The values of the selected elements can be plotted in a fracture locus, cf. figure 4.5(a). The graph contains the data points of both sample types, that is, thin and thick. Despite the tremendous reduction of data, the states of stress and of strain are still very heterogeneous. For the thin sample, the location of maximum equivalent plastic strain and maximum triaxiality coincide and it seems very likely that this combination leads to the fracture initiation. For the thick sample the plastic strain is almost uniform in the ROI over the entire triaxiality range. The "maximum" equivalent plastic strain is located at 𝜂 = 0.48. However, this state cannot induce fracture, as it also exists in the thin sample. It seems therefore more likely that failure is associated with the strains at the highest triaxiality.

The definition of a criterion is difficult in this case. If the two points were to be connected by a line this would indicate that the material fails almost at a critical triaxiality.

This contradicts the findings in the literature, namely, ductile materials usually show a very distinctive dependence of the maximum plastic strain on triaxiality. This fracture locus does not seem to be appropriate for the current tests of TA6V.

Another proposition for a fracture locus from the literature is the use of the maximum principal stress and equivalent plastic strain, cf. figure 4.5(b). In contrast to literature results, no influence of the level of the equivalent plastic strain is found. A limit criterion was drawn at 𝜎 1 = 1630 𝑀 𝑃 𝑎. The fracture of TA6V could only be sensitive to the maximum principal stress, given that the computation should consider plasticity and large strains.

For ductile damage, it is possible that the triaxiality increases the amount of damage that exists for a certain equivalent plastic strain. For the Lemaitre damage law, the influence is described by the triaxiality function 𝑅 𝜈 . Due to the uncertainty on the role of damage we propose a simplified damage indicator, which weights the equivalent plastic strain with the triaxiality function,

𝐷 𝑖 = 𝑝𝑅 𝜈 .
(4.2)

Figure 4.6 presents two possible variations for the triaxiality function, the original Lemaitre function 𝑅 𝜈 1 and a proposition from the literature 𝑅 𝜈 4 , cf. section 2.5.3,

𝑅 𝜈 1 = 2 3 (1 + 𝜈) + 3(1 -2𝜈) (𝜂) 2 , 𝑅 𝜈 4 = (1 + 𝛽 𝐷 [𝜂 2 -(1/3) 2 ]), (4.3)
where the parameters are 𝜈 = 0.32 and 𝛽 𝐷 = 6. Compared to the equivalent plastic strain -triaxiality locus, the function 𝑅 𝜈 1 only marginally increases the importance at high strains. The proposition 𝑅 𝜈 4 imposes a much stronger modification. For both samples the maximum damage correlates with the maximum triaxiality.

Other propositions for the fracture locus

As the construction of the fracture locus is based on I-DIC simulations, the choice of any other mechanical property is possible.

1968), ˆ𝑝𝑚𝑎𝑥

0

𝜎 1 𝜎 𝑒𝑞 d𝑝 = 𝐶 𝐿𝐶 , (4.4)
where 𝐶 𝐿𝐶 is a material constant. Similar to the combination with plasticity, almost no influence of the von Mises stress exists. The same limit criterion as above was drawn in the figure .  A criterion, which uses the principal strain was also proposed [START_REF] Bourgeon | Etude et modélisation des mécanismes d'endommagement en forge à froid[END_REF]. The circumferential strain and the axial strain at fracture in a smooth round tensile sample are related by ε𝜃𝜃 + 1 2 ε𝑧𝑧 = 𝐶, (4.5)

where 𝐶 is a material constant. The use of such a criterion seems equivalent to a dependence of the maximum principal stress to the other principal stresses. If the magnitudes of the minimum and the intermediate principal stresses become closer to the maximum principal stress, the triaxiality becomes higher. We therefore used the maximum principal stress and the triaxiality, cf. figure 4.7(c). However, the critical points are too close together to identify any variation. Instead, the same criterion as before was drawn in the diagram.

The exact relationship between the diagrams is complex. It seems therefore necessary to investigate all diagrams together to acquire a comprehensive overview on the fracture cause.

Validation of the fracture locus

The final step in the proposed procedure is the evaluation of the fracture locus for the entire data. Figure 4.8 shows this step for the last step of the analysis. It was checked that the final step leads to higher values than the preceding ones, cf. appendix I.

For the maximum principal stress -equivalent plastic strain locus, all points are below the limit lines. This means the selected elements contained the fracture surface. For the locus equivalent plastic strain -triaxiality, no obvious criterion can be deduced. However, 4.2 Characterization of the fracture of U500

Microscopic and macroscopic analyses

The previous methodology was applied to the U500 test results. The examination of the fracture surface clearly reveals the structure of the grains. In addition, the macroscopic force elongation curve for any sample showed no significant softening. Both results are as expected for brittle fracture. Therefore no special cyclic tests were performed. Nevertheless, the fracture surface is not straight. Larger holes mark the surface, i.e. close to the notch at point (1). The occurrence of bands with a distorted microstructure is the microscopic cause of the mentioned macroscopic anisotropy (2). On closer inspection of the grains, one finds a coarse texture, which covers the entire grain surface (3). The most likely explanation is the existence of micro-voids on the grain boundary.

Construction of the fracture locus

The selection of the elements in the vicinity of the fracture surface was performed analogously to the selection for TA6V, cf. figure 4.10. The same fracture loci as for TA6V were evaluated. Even though the material fails in a brittle manner, the fracture locus usually employed for ductile rupture, namely equivalent plastic strain -triaxiality, is quite reasonable. For the thin sample, the maximum equivalent plastic strain is not located at the maximum triaxiality. The points seem to follow to the tendency described in [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF][START_REF] Papasidero | Ductile fracture of aluminum 2024-T351 under proportional and non-proportional multi-axial loading: Bao-Wierzbicki results revisited[END_REF] and shown in figure 2.11(a, b). Due to the lack of data at a triaxiality 𝜂 > 0.6 no criterion was defined. Further tests at even higher triaxialities are necessary.

Similar to TA6V, the locus of maximum principal stress -equivalent plastic strain allowed for the construction of a limit line at 𝜎 1 = 1830 𝑀 𝑃 𝑎. Such a criterion seems therefore also relevant for U500. Let us emphasize again that the results implicitly depend on the plastic response and that the calculation should be performed with large strains. For the evaluation of the critical damage indicator fracture locus the proposition 𝑅 𝜈 4 leads to a critical value of the damage indicator. We chose the maximum value of the thin sample due to the higher precision of its identification. The coefficients were chosen as 𝜈 = 0.3 and 𝛽 𝐷 = 6. For the first proposition the modification is too weak to change substantially from the locus with the equivalent plastic strain. The three additional propositions for the fracture loci are presented in figure 4.13(a-c). The results are globally similar to those of TA6V. However, the points of the two sample types are farther apart. The limit line was drawn for the loci that depend on the maximum principal stress. In figure 4.13(d) the criterion maximum principal stress -equivalent plastic strain was satisfied for all data points. Let us note that the failure surface was idealized and that as a result the approximation is more pronounced in the thick samples than in the thin samples. In fact, there are some points for the thick sample, which are slightly above the line.

Summary

This chapter concerns the failure behavior of TA6V and U500. In a first step the fractographies of TA6V and Udimet 500 are presented. TA6V ruptures in a transgranular manner; no grains are visible. The fracture surface indicates that void nucleation and growth are the primary reason of the found fracture. No indication of macroscopic damage was found before the occurrence of failure. For U500, the grains are visible and the fracture is intergranular. Therefore in a first instance it is assumed that for both materials the material failure can be appropriately described by a criterion. The eventual occurrence of microscopic damage is neglected.

The chapter then focuses on the use of I-DIC approaches with regard to the identification of possible failure criteria. A procedure was developed to analyze the mechanical fields in combination with the image of the fractured sample. The elements from the fracture surface are extracted, i.e. in the proximity of an idealized surface. The data of the extracted elements can then be analyzed to construct a fracture locus.

Several propositions for failure criteria from the literature were tested, namely, equivalent plastic strain -triaxiality and maximum principal stress -equivalent plastic strain. A damage indicator was proposed, which weights the equivalent plastic strain with a triaxiality function.

Other alternatives were studied and it appears that the maximum principal stress criterion fits very well the extruded data keeping in mind that any such analysis should consider plasticity calculations. In addition, it has to be noted that the fact that only one point for thick and thin samples reaches the limit value is attributed to the fact that once the failure initiates it propagates in a catastrophic manner. The analysis of all the data in time and space available from the I-DIC analysis was conducted to confirm that the criteria were not reached before and/or in other locations of the tested samples.

Perspectives for future applications

The previous conclusion, especially the one concerning the maximum principal stress criteria should be checked on other test cases, on sample geometries involving other and various states of stresses.

The simultaneous exploitation of all four cameras could improve the overall I-DIC residuals, which in turn could lead to a better identification. Of course a step forward would be to analyze the result on the experimental fracture surface and not on an idealized one.

Another improvement of the method would be the introduction of these fracture criteria in FE simulations. The erosion of elements from the simulation could lead to a virtual fracture surface. With additional procedures these fracture surfaces could be compared to the recorded images of the fracture. One could also check whether once the failure initiates, it propagates in a catastrophic manner Further, with an adapted I-DIC procedure, the crack can be recreated in the computergenerated images. Such a procedure would enable the correlation of the fracture image with the I-DIC approach. The identification of plasticity and fracture could be performed simultaneously. This is similar to the ideas issued by [START_REF] Roux | Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2D and 3D cracks[END_REF]Cao et al., 2014b) which introduced more local information to the identification functional.

Another important point considers the identified constitutive model. For the thin and thick samples the average curves were not identical. First, the reason for such a difference needs to be investigated. In this context, the sensitivity of the triaxiality fields to minimal changes in the identified law should be studied. Second, it should be possible to identify a unique constitutive model for all samples in a simultaneous identification.

Together with the described improvements, the model identification for plasticity and fracture would then be one "big data" identification.

CHAPTER 5

Perspectives for the simulation of material failure Compared to the use of global criteria, local approaches may give insight into the initiation of failure and the possible subsequent crack pattern in the disk if subjected to overspeed burst. This could help the engineer to improve the design. From such approaches a potential gain in the precision of the limit speed prediction could also be expected, which is the main concern of this thesis.

Depending whether the failure is brittle or ductile, two types of situations have to be considered. Let us recall however, that, from the analysis of the tests conducted in this thesis on TA6V and Udimet 500, no clear indication of ductile failure was obtained.

Plastic strain leads to a decrease in the local stiffness. This effect is increased by the spin softening, in case of the disks, or by the existence of damage, in case of ductile failure. If too much softening occurs, the non-convergence of implicit static simulations takes place for most algorithms. A simulation, which does not terminate in a stable predefined manner, presents a major inconvenience for the numerical analyst. If a local damage model is used, an even more problematic phenomenon occurs, the one of spurious mesh dependence.

Regularization of the mesh dependence

Background

In today's industrial environment, due to the lack of robustness of material failure models with respect to structural simulation (called spurious mesh dependence in the literature), important numerical parameters such as the mesh density are fixed in order to calibrate the model with respect to some reference test. Due to the importance of the topic, a large body of literature has been devoted to the problem of objective simulation of failure since the beginning of the 80's.

The following paragraph proposes a very brief overview of non-local approaches, the most widespread approach to overcome the lack of consistency of the material model with respect to the simulation of failure. The focus is then shifted to rate dependent models and more specifically to bounded rate models.

Non-local model

The fundamental problem of spurious mesh dependence has been studied for more than thirty years now and is well understood and analyzed. A huge literature has been devoted 81 to non-local models as non-local integral approaches, explicit and implicit gradient model or Cosserat models [START_REF] Bazant | Instability, ductility, and size effect in strain-softening concrete[END_REF][START_REF] Borst | Simulation of srain localization: a reappraisal of the Cosserat continuum[END_REF][START_REF] Lasry | Localization limiters in transient problems[END_REF][START_REF] Pijaudier-Cabot | Nonlocal Damage Theory[END_REF]. Reviews are available for integral approaches, [START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: survey of progress[END_REF], and for various aspect of damage formulation in this context in [START_REF] Besson | Continuum models of ductile fracture: a review[END_REF].

Many improvements have been made since the middle of the eighties and non-local approaches are now well-mastered from the mathematical and numerical point of view. The development of non-local approaches in an industrial context is nevertheless seldom. The main reason is probably the fact that non-locality implies many and non-obvious code developments. The identification practice is also an issue. Therefore, a large body of theories is under development to relate internal length scales and associated boundary conditions with the underlying physics at small scale see e. g. [START_REF] Al-Rub | A Finite Strain Plastic-damage Model for High Velocity Impact using Combined Viscosity and Gradient Localization Limiters: Part I -Theoretical Formulation[END_REF][START_REF] Ricci | Numerical Analysis of Nonlocal Anisotropic Continuum Damage[END_REF]. Let us nevertheless note that an example for an implicit implementation in ABAQUS/STANDARD is proposed in [START_REF] Sornin | About Elastoplastic Nonlocal Formulations with Damage Gradients[END_REF]. The authors reported convergence difficulties. An implementation of non-local models for ABAQUS/EXPLICIT was proposed in [START_REF] Saanouni | Micromorphic approach for finite gradientelastoplasticity fully coupled with ductile damage: Formulation and computational aspects[END_REF].

Rate dependent model

Needleman [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF] was possibly the first to discuss how, in statics, the use of viscosity can help to conserve the elliptic property of the incremental equilibrium equations and thus eliminate pathological mesh-sensitivity. Several models have been proposed in order to control localization using viscosity. In particular, this approach was applied to ductile materials with negative hardening [START_REF] Sluys | Wave propagation and localization in a rate-dependent cracked medium-model formulation and one-dimensional examples[END_REF]). Nevertheless, it has been observed that the crack growth behavior predicted by simulations based on a visco-plastic version of the Gurson-Tvergaard-Needleman (GTN) model, is, in general, mesh sensitive [START_REF] Needleman | Mesh effects in the analysis of dynamic ductile crack growth[END_REF].

Other experimentations have led to deceptive results, as shown for example in the work of [START_REF] Comi | On gradient regularization for numerical analyses in the presence of damage[END_REF]. More recently, an extension of the Johnson-Cook model with damage was carefully studied [START_REF] Flatten | Non-local modeling of thermomechanical localization in metals[END_REF][START_REF] Flatten | Lokale und nicht-lokale Modellierung und Simulation thermomechanischer Lokalisierung mit Schädigung für metallische Werkstoffe unter Hochgeschwindigkeitsbeanspruchungen[END_REF]. Here again, it was shown and explained that pathological mesh dependence is not prevented, even though the model is a rate-dependent one. This seems in conflict with theoretical studies, showing that the use of viscosity allows the problem to remain hyperbolic in dynamics. In some cases, it was shown that for the time-discretized problem, the difficulty could be a pure numerical one [START_REF] Benallal | A note on ill-posedness for rate-dependent problems and its relation to the rate-independent case[END_REF]. A critical time step, much lower than the ones generally used in fracture simulations, is needed to achieve sensible results.

The bounded rate approach

In the case of infinitesimal strain and for an elasto-damage model, [START_REF] Allix | Delay damage modelling for fracture prediction of laminated composites under dynamic loading[END_REF] have tried to overcome the limitation of viscous-regularization by introducing the concept of a bounded damage rate model. A physical interpretation of the model is that a continuous damage variable results from the average effects of micro flaws. As each flaw has a finite propagation velocity, "any" averaging process should give rise to a bounded rate of damage. This idea is to be related with the concept of incubation time introduced in failure criteria for dynamic loading [START_REF] Curran | Dynamic failure of solids[END_REF]. The introduction of the maximal damage rate intervenes only during the localization phase, during which the damage rate may approach its maximal value. If the material exhibits more classical viscous effects, these should also be incorporated. In the numerical experimentation performed with the bounded rate damage model no spurious mesh dependence was observed.

In [START_REF] Guimard | Characterization and modeling of rate effects in the dynamic propagation of mode-II delamination in composite laminates[END_REF] experiments on mode II dynamic delamination allowed the identification of the model for interlaminar interfacial damages. An explanation of the regularization effect of bounded rate damage model was proposed in [START_REF] Allix | The bounded rate concept: A framework to deal with objective failure predictions in dynamic within a local constitutive model[END_REF]. A comparison of the capabilities to deal with localization in concrete was performed for non-local and bounded rate damage model [START_REF] Desmorat | Delay-active damage versus non-local enhancement for anisotropic damage dynamics computations with alternated loading[END_REF].

Bounded damage rate

For ductile failure, Suffis and Combescure developed, for the model proposed in [START_REF] Allix | Delay damage modelling for fracture prediction of laminated composites under dynamic loading[END_REF], an original and efficient way to derive an estimation of the characteristic length associated with the bounded plasticity rate model [START_REF] Suffis | Damage model with delay effect[END_REF].

A bounded rate version of the Lemaitre damage model is defined as:

Ḋ = 1 𝜏 𝑐 [1 -𝛷] , 𝛷 = 𝑒𝑥𝑝(-𝑎⟨𝐷 𝑠 -𝐷⟩), 𝐷 < 𝐷 𝑐 .
(5.1)

The standard static damage expression 𝐷 𝑠 serves as driving force to the actual damage 𝐷. The maximum rate is Ḋ𝑚𝑎𝑥 = 1 𝜏𝑐 . Where, 𝜏 𝑐 is the critical time necessary for the physical phenomena to take place: the nucleation of voids. The parameter 𝑎 defines the size of the damage zone. The form of the static damage law is independent of this formulation.

Bounded plasticity rate

In the context of finite strain and for the prediction of the deterioration of metallic parts subjected to bird impact, spurious mesh dependence was still observed using the above model (Court, 2006). Plasticity and damage induce two sources of instability, but only one is controlled by the use of a bounded rate damage model.

A detailed description of the physics of ductile failure in the case of aluminum can be found in [START_REF] Ghahremaninezhad | Ductile failure behavior of polycrystalline Al 6061-T6[END_REF]. The appearance of localized necking gives rise to the nucleation and coalescence of micro-voids and micro-cracks up to the formation of a macro-crack, cf. section 2.5. The onset of localized necking and damage can be predicted as in [START_REF] Chow | A Damage-coupled Criterion of Localized Necking Based on Acoustic Tensor[END_REF]. Nevertheless, if one wishes to predict the final state of the structure, a model is necessary, which allows mimicking the whole deterioration process.

In the framework of the bounded rate approach, a finite strain model was proposed, where the damage variable is a function of the equivalent plastic strain, whose rate is bounded and governs the damage evolution:

ṗ = 1 𝜏 𝑐 [1 -𝛷] , 𝛷 = 𝑒𝑥𝑝(- 𝑎 𝜎 𝑦 ⟨𝑓 ⟩), 𝐷 = 𝑝 𝑝 𝑐 𝑅 𝜈 . (5.2)
Where, 𝜏 𝑐 , the critical time, is associated to the movement of the dislocations. The driving force is the yield function 𝑓 , normalized with respect to the yield stress 𝜎 𝑦 . The parameter 𝑎 defines the more or less brittle character of the failure process. Only positive values of 𝑓 lead to a plastic strain.

Application of the bounded plasticity rate model in dynamics

The model has allowed preventing mesh dependence while restricting code developments to local modifications. Figure 5.1 shows the damage field of a plate under tensile loading. The plate is tapered so that the maximum load is situated on the left side. For the "classical" damage formulation the mesh refinement leads to a stronger localization of the ruptured zone (𝐷 = 1). With the bounded plasticity rate approach, the damaged zone remains constant. For dynamic tests ( ε = 25 𝑠 -1 ), which led to quite different failure scenarios than in this thesis, 2D simulations have been satisfactory compared to test results both on the prediction of the damage (equivalent crack path) and on the time to failure, cf. figure 5.2. For the classical formulation, the simulation depended on the mesh.

Implementation of the bounded plasticity rate models

The application of the bounded plasticity rate model to quasi-static phenomena poses a number of questions. First, the regularization using bounded rate model is achieved only if inertia forces are included in the analysis. Therefore, the simulations should be done in dynamics, at least for the failure phase where dynamic effects are locally important even for quasi-static loading. This is usually easily feasible as most commercial codes for static calculations provide an implicit dynamic solver. Moreover, in order to obtain sensible results, the characteristic time of the model should be low. This requires the use of very small time steps to properly follow the model. Finally, possible mesh distortions during the failure phase may lead to additional burdens.

Time step restriction

The limitation of the rate follows directly from the solution of the differential equation:

ṗ = 1 𝜏 𝑐 [1 -𝛷] , 𝛷 = 𝑒𝑥𝑝(- 𝑎 𝜎 𝑦 ⟨𝑓 ⟩).
(5.3)

In contrast to the previous application (Court, 2006), the model was implemented in a user-material, employing the theta-method. The solution of the constitutive model can be controlled as a mix of implicit and explicit solution. Details on the implementations are found in appendix J.

To effectively limit the plastic rate, the bounded plasticity rate model needs to be calculated at a time step inferior to the critical time 𝜏 𝑐 . This restriction can lead to long simulation times if global implicit calculations are used. The use of a global explicit simulation is a possible solution.

Coupling to damage and limitation of mesh distortion

As proposed in [START_REF] Allix | The bounded rate concept: A framework to deal with objective failure predictions in dynamic within a local constitutive model[END_REF], the damage variable directly depends on the equivalent plastic strain in the proposed formulation. In a first instance, the influence of triaxiality on the damage variable is neglected. The damage rate is then directly limited through the equivalent plastic strain rate. Due to the unclear and/or late onset of damage in the disk materials a threshold 𝑝 0 is included,

𝐷 = 𝑝 -𝑝 0 𝑝 𝑐 -𝑝 0 .
(5.4)

The common implementation of the Lemaitre damage model assumes that the effective stress σ is written together with the effective equivalent plastic strain p (Souza Neto, 2002). The yield function is defined as

𝑓 = 𝐽 2 (σ) -𝑅(p) = 0.
(5.5)

This means the macroscopic equivalent plastic strain 𝑝 becomes 𝑝 = p 1 -𝐷 .

(5.6)

If the critical damage is assumed close to one, the failure occurs at a multiple of the critical equivalent plastic strain 𝑝 𝑐 . This can lead to very large distortions of the mesh. This issue can be avoided by limiting the critical damage to lower values 𝐷 𝑐 < 0.3 (Souza Neto, 2002) or through remeshing.

In addition, we propose the combination of effective stress and real macroscopic equivalent plastic strain to reduce the macroscopic strain at failure,

𝑓 = 𝐽 2 (σ) -𝑅(𝑝) = 0
(5.7)

Figure 5.3 presents the evolution of total deformation and macroscopic stress for both formulations for linear hardening. The calculation was carried out by imposing a constant strain increment in one Gauss point. The material properties are Young's modulus 𝐸 = 100000 𝑀 𝑃 𝑎, plastic modulus 𝐶 = 200 𝑀 𝑃 𝑎, yield stress 𝜎 𝑦 = 800 𝑀 𝑃 𝑎 and the damage parameters 𝑝 0 = 0.08, 𝑝 𝑐 = 0.7, 𝐷 𝑐 = 0.99. The standard formulation leads to a total deformation greater than 150 %. The proposed formulation terminates close to 𝜀 𝑡 ≈ 𝑝 𝑐 , which limits the mesh distortion.

Application to axisymmetric cases with centrifugal loading

When studying the problem of localization a first common step is to examine the case of a bar under tensile loading. Nevertheless this case is not representative to the problem at hand and we seek for an equivalent case of a disk submitted to centrifugal loading, i.e. a disk of uniform thickness treated under axisymmetric hypothesis. It appears, as can be seen from the discussion in section 5.2.1, that this case is of limited interest. Therefore, a second test case has been studied in section 5.2.2, the one of a disk with a notch treated under axisymmetric hypothesis.

Simple disk with uniform axisymmetric cross-section Axisymmetric model of a circular disk with a uniform thickness

The disk has a radius of 100 mm and a thickness of 20 mm. Axisymmetry and mid-plane symmetry were used to simplify the model as much as possible, as is displayed with a coarse mesh of 𝑡 𝑒𝑙 = 1 𝑚𝑚 in figure 5.4. In contrast to a two dimensional bar, the volume of the elements is not identical in a uniform mesh, but varies along the radial axis. Consequently, if loaded under tension in radial direction, the element closest to the axis of rotation undergoes the highest stress. The model depicted in figure 5.4 is therefore similar to a bar of progressively growing cross-section. Another fundamental difference from the tensile bar is the influence of the centrifugal force. In the tensile bar, the internal force is constant for all cross-sections, whereas in the disk the axial force increases from the external rim to the axis of rotation. These geometric and loading conditions lead to a concentration of the stresses and consequently of the damage in the inner part of the disk. This example is therefore not very well suited to properly analyze localization properties of the models as will be obvious from the results.

The quasi-static character of the burst test was neglected and the loading was dynamic. The model was loaded with a constant rotational acceleration of θ = 10.5 × 10 6 𝑟𝑎𝑑/𝑠 2 . The results are therefore not representative of a burst test.

Material parameters

The density was 𝜌 = 4.5 𝑡𝑜𝑛 -9 /𝑚𝑚 3 and the elasto-plastic behavior was assumed non-linear using the extended Voce law: 𝐸 = 110000 𝑀 𝑃 𝑎, 𝜎 𝑦 = 800 𝑀 𝑝𝑎, 𝑅 ∞ = 1300 𝑀 𝑃 𝑎, 𝛾 = 3, 𝐶 = 200 𝑀 𝑃 𝑎. The damage parameters were taken as 𝑝 0 = 0.25, 𝑝 𝑐 = 0.5. Damage is linear and triaxiality independent for all models. The critical damage is set to 𝐷 𝑐 = 0.9, which defines a final plasticity value of 𝑝 𝑚𝑎𝑥 = 0.475.

For the bounded plasticity rate model, the critical time was chosen at approximately half of the maximum rate that occurred in non-limited calculations. The limit parameters are 𝑎 = 1 and 𝜏 𝑐 = 10 -4 𝑠.

Choice of the element size Three different meshes were used in the analysis with uniform element size of 𝑡 𝑒𝑙 = [0.25, 0.5, 1 𝑚𝑚], cf. figure 5.5(a-c). This lead to a time step in the order of 10 -8 𝑠 for the fine mesh. Standard integration with four Gauss points was used for the linear elements.

The convergence of the three meshes was verified for elasto-plasticity. Figure 5.5(d) shows the influence of the element size on internal energy, displacement of the external radius and equivalent plastic strain in the lower left corner. The values are normalized to the solution for the smallest element size. The internal energy is constant for all three meshes. The coarsest mesh leads to a variation of the equivalent plastic strain and displacement of less than 0.5 %. Computation including damage: results with a classical and bounded rate models Figure 5.6 shows the evolution of equivalent plastic strain for the two models for the three meshes in the element in the lower left corner, which is the most loaded. The curves were drawn up to the moment the critical damage in the element reaches 𝐷 𝑐 = 0.9, i.e. 𝑝 = 0.475.

As is clearly visible the strain approaches an instability. The use of the bounded plasticity rate model leads to a delay of this instability. In addition the slope of the evolution, the rate, is less strong at fracture of the element. This is a result of the limiting. Despite of these differences between the models, all three meshes lead to identical results for both models. No variation can be perceived at this global level. The evolution of the damage field for the classical model and the bounded plastic ones are shown in figure 5.7 and figure 5.8, starting from the moment the first element reaches the critical damage value. To highlight the precise spreading of damage the scale is fixed to 0.8 < 𝐷 < 0.9. The images are limited to the center of the disk, i.e. the left side of the model. The only clear distinction between the two models is slightly delayed occurrence of damage using the bounded rate model. In fact, the effect of the centrifugal forces leads to progressive damage gradient in the elements from the central axis to the external rim, whatever is the damage model.

Axisymmetric disk with a notch

In order to study a more meaningful case with respect to localization we consider now the case of a notched disk shown in figure 5.9 in the case of a coarse mesh (1 𝑚𝑚 mm element size). The main thickness of the disk is 30 𝑚𝑚 and near the radial position 𝑟 = 50 𝑚𝑚 its thickness is progressively reduced to 20 𝑚𝑚. A finer mesh with an element size of 0.5 𝑚𝑚 was used for the analysis of possible mesh dependence. All other parameters were unchanged compared to the case of the disk with a uniform thickness.

First analysis

This case led to a problem not encountered in the case of the disk of uniform thickness the fact that the fully damaged elements encounter large distortion which from a certain instant leads to non convergence issues. The damage fields are presented in figure 5.10 just before non-convergence for the fine meshes. For the classical damage model this corresponded to θ𝑟𝑒𝑓 = 9747 𝑟𝑎𝑑/𝑠 and for the limited plasticity model to θ𝑟𝑒𝑓 = 10046 𝑟𝑎𝑑/𝑠.

It can be seen that while for the limited plasticity rate damage model, the results are nearly independent of the mesh size this is not the case for the classical damage model. In this case, the coarsest mesh leads to a localization of damage, while the crack extends to half of the disk's thickness for the finest mesh. Such a large difference is a rather unexpected. It is the extent of the damage area perpendicular to the direction of the crack that usually poses problems with a classical model. This goes as far as a possible dependence on the mesh orientation. It is therefore suspected that the problem occurs due to the effect of local mesh distortion, which eventually causes the problem of convergence. Therefore the calculation was performed a second time including element deletion.

Simulations with element deletion

The previous simulations fail shortly after the rupture of the first elements, due to mesh distortion. Coarser meshes can accommodate more distortion and the calculation carries on longer. To avoid possible effect of mesh distortion the elements were deleted at the failure of the first Gauss point. The figure 5.11 shows the evolution of the damage field around the notch for the classical damage model at four different instants corresponding to the onset of the crack and different stages of its propagation up to the full fracturing of the disk. Let us note that here the color scale used for the figures covers the full ranges of damage 0 < 𝐷 < 0.9.

Final fracture occurs at θ = 9818 𝑟𝑎𝑑/𝑠 for the coarse mesh and at θ = 9734 𝑟𝑎𝑑/𝑠 for the fine mesh. This latter speed is only minimal above the speed for the initial fracture of the coarse mesh θ = 9723 𝑟𝑎𝑑/𝑠. In addition only a row of one or two elements in the fine mesh has undergone rupture, which is a result as expected. Such a mesh dependence is responsible to unrealistic dependence of the global results in other type of applications, as for example the applied forces at failure.

For the bounded rate model quite similar results are obtained for the different meshes at the same instants, cf. figure 5.12. Only the fracture of the first element varies considerably, θ = 10027 𝑟𝑎𝑑/𝑠 compared to θ = 9912 𝑟𝑎𝑑/𝑠. This seems to indicate that the coarse mesh is not fine enough. However, the results at the second instant are very similar. The fracture of the disk occurs for both meshes at θ = 10227 𝑟𝑎𝑑/𝑠. The different size of the crack in the images, can be explained with the very high level of damage in the intact elements. The actual size of the damaged zone with 𝐷 > 0.8 is quite equal in both meshes. Computation for even finer meshes should be performed to fully reach convergence.

These observations are confirmed by the fields of the equivalent plastic strain rate presented in figure 5.13 and figure 5.14. For the classical damage model the scale limit was fixed to the maximum ṗ = 1.7 × 10 5 /𝑠, while for the bounded rate model it cannot exceed ṗ = 1 × 10 4 /𝑠. The largest values for the classical model is reached just before fracture of the whole disk in the fine mesh at θ = 9713𝑟𝑎𝑑/𝑠. However, no other element attains equally high values. The limitation of the bounded rate model leads to a more homogenous distribution of the equivalent plastic strain rate.

Many of the elements in the crack path deform at the maximum rate which is required to obtain mesh independent results.

Use of explicit dynamics for quasi-static problems

The objective of this section is the study of the use of mass and time scaling in explicit simulations for the prediction of the burst speed. As already explained, explicit dynamics provides a gain in robustness. This concerns firstly elasto-plastic simulations including spin softening effects. Higher levels of equivalent plastic strain can be calculated. Another motivation for the application of dynamics is the possibility to switch from a quasi-static solution to a real dynamic by introducing failure to the calculation. Doing so, offers the possibility to benefit from the regularization effect of the bounded rate model. We first start to recall the main properties of the central difference scheme, cf. section 5.3.1. This is followed by a review of several methods proposed in the literature that enable the application of the central difference scheme to quasi static load case, cf. section 5.3.2. In order to examine the effects on the prediction of the critical burst speed and failure, three examples are considered. The first one concerns the simple case of a non-linear elasto-plastic mass spring system including the spin softening effect, cf. section 5.3.3. The second one concerns the application of these techniques to an axisymmetric disk model with tensile loading in order to check the potential of the approach on the structures of interest, cf. section 5.3.4. Finally, the disk with a notch is reexamined to discuss how quasi-static mesh independent rupture can be assured, cf. section 5.3.5.

Central Difference scheme

The central difference scheme is a particular version of the Newmark scheme with 𝛾 = 1/2 and 𝛽 = 0. This scheme is probably the most used algorithm for impact problems because the current nodal displacements 𝑛+1 U, cf. equation (5.8), depend only on nodal quantities computed in the previous time step at 𝑡 𝑛 . The current configuration does not change with the algebraic computations of the unknown accelerations 𝑛+1 A and the nonlinearities, as in contact problems, do not require iterative computations of the solution. System (5.8) reads ∀𝑛 ∈ {0, . . . ,𝑛 -1}:

M 𝑛+1 A = 𝑛+1 F ext -F int ( 𝑛+1 U) in 𝛺 ℎ | 𝑡 𝑛+1 , 𝑛+1 U = 𝑛 U + 𝑛 V d𝑡 + 𝑛 A d𝑡 2 2 in 𝛺 ℎ | 𝑡 𝑛+1 , 𝑛+1 V = 𝑛 V + 𝑛 A d𝑡 2 + 𝑛+1 A d𝑡 2 in 𝛺 ℎ | 𝑡 𝑛+1 , 𝑛+1 U = U over 𝜕𝛺 𝑢 | 𝑡 𝑛+1 , { 0 U, 0 V} = {U 0 ,V 0 } in 𝛺 ℎ | 𝑡 in .
(5.8)

In the central difference scheme the use of the lumped mass matrix M is preferable to the consistent one M, because it increases the stable time step and because it does not need any matrix inversion, i.e. LU factorization. This allows a very efficient resolution of the system (5.8), in which the single degrees of freedom are computed solving single independent equations. Such a system of independent equations can be solved with high performance architectures and shared-memory multi-processing loops. As described in [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF], the lumped mass matrix M can be obtained using:

• the nodal quadrature rules

• the row-sum technique

• the Hinton-Rock-Zienkiewicz (HRZ) technique: only the diagonal terms of the consistent mass matrix M are computed and then scaled to preserve the total element mass Each one of the mentioned techniques presents advantages and disadvantages, as summarized in table 5.1. Nodal quadrature rules and the row-sum technique applied to high-order elements determine negative or zero values at certain nodes. On the other hand, the HRZ technique avoids such values, but so far has not been supported by mathematical proofs. [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF]. Nevertheless, the application to simple elements, like linear triangles or bi-linear quadrilaterals, produce similar (if not identical) lumped mass matrices M, without negative or zero values. The Abaqus/Explicit element library, for instance, offers only linear shape functions with reduced integration.

Lumping technique

A major problem of explicit time integration schemes is their conditional stability. A stable time step 𝛥𝑡 𝑐𝑟𝑖𝑡 for an undamped problem is estimated as

𝛥𝑡 𝑐𝑟𝑖𝑡 ≤ 2 𝜔 max , (5.9)
where 𝜔 max is the maximum eigenfrequency of the spatial discretization 𝛺 ℎ . It can be demonstrated that (5.10) where 𝜔 el are the eigenfrequencies of the individual elements "el" of the mesh 𝛺 ℎ . For elasto-plastic and damage models, the effect of the constitutive relation leads to a reduction of those frequencies. This circumstance allows estimating an upper-bound for the time step considering elasticity only. Nevertheless, large displacements require its actualization.

𝜔 max ≤ max el {𝜔 el },
Only undamped systems are considered in this work. Even if common sense would suggest one to think differently, viscous damping reduces the stable time step, so that for elasticity in the geometrically linear regime and for linear damping .11) where 𝜔 max is again the maximum eigenfrequency of the undamped system and 𝜁 is the fraction of critical damping in the highest mode. Analytical expressions for bounding the maximum eigenfrequency 𝜔 max were given in [START_REF] Flanagan | Eigenvalues and Stable Time Steps for the Uniform Strain Hexahedron and Quadrilateral[END_REF] for the case of reduced-integration quadrilateral and hexahedral elements with undistorted or distorted shapes.

𝛥𝑡 𝑐𝑟𝑖𝑡 ≤ 2 𝜔 max ( √ 1 + 𝜁 2 -𝜁 ) , ( 5
Furthermore, a physical interpretation of the stability limit was introduced as Courant-Friedrichs-Lewy (CFL) condition [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF]. It prescribes that the stable time step 𝛥𝑡 𝑐𝑟𝑖𝑡 must be smaller than the time required by a dilatational stress wave (that is always bigger than an eventual shear wave) to traverse the smallest element of the given mesh 𝛺 ℎ :

𝛥𝑡 𝑐𝑟𝑖𝑡 ≤ 𝐿 el 𝑐 , (5.12)
where 𝐿 el is the characteristic length of the element and 𝑐 is the wave propagation speed. The characteristic length is estimated differently depending on the element technology. Some examples are summarized in table 5.2 from the LS-DYNA theoretical manual (LS-DYNA, 2014), where the estimations are defined considering the efficiency of their calculation.

Mass and time scaling techniques

The usual acceleration of the burst speed is about 10 𝑟𝑎𝑑/𝑠 2 . This results in total test times between 500 and 900 seconds, depending on material and geometry. For the meshes used in the dynamic calculations above the time step was of the order of 𝛥𝑡 < 1 × 10 -7 𝑠.

If the real test was to be simulated with this time step, this would lead up to 9 billion time steps.

Several techniques allow the reduction of the simulation time. First, it is possible to switch from the implicit solution to an explicit solution. The elastic phase and stable plastic phase can be covered in a relatively short time. However, today only few commercial codes allow this switching directly. An exception is the dynamics solver LS-DYNA, which was used for the axisymmetric calculations presented in this study. For other software packages (e.g. ABAQUS), the switching requires a considerable amount of data conversion between two codes.

Second, two widely used techniques are mass and time scaling. These techniques alter the model in such a way that the total simulation time decreases. However, the effect of these procedures remains to be fully studied, therefore recent attempts try to standardize the application procedure of explicit dynamics to quasi-static problems [START_REF] Gulavani | Explicit dynamic formulation to demonstrate compliance against quasi-static aircraft seat certification loads (CS25.561) -Part I: influence of time and mass scaling[END_REF][START_REF] Hughes | Explicit dynamic formulation to demonstrate compliance against quasi-static aircraft seat certification loads (CS25.561) -Part II: Influence of body blocks[END_REF]. This section proposes an overview of those techniques.

Mass scaling

From table 5.3 it is obvious that an increase in density leads to an increase in the critical time step, this approach being called mass scaling. The effect of mass scaling on the time is defined by a mass scaling factor 𝛼 𝑀 which modifies the density quadratically,

𝛥𝑡 𝑐𝑟𝑖𝑡 ≤ 𝑙 √ 𝜌𝛽 𝑀 𝑆 𝐸 , 𝛽 𝑀 𝑆 = 𝛼 2 𝑀 .
(5.13)

Several techniques exist for the application of the mass scaling. The simplest application is the uniform scaling of the mass. A more elaborate procedure is the scaling of every separate element to at least the target time step. The global mass scaling factor 𝛼 𝑀 is then not identical to the elemental factor 𝛼 𝑒𝑙 𝑀 . This procedure is known as selective mass scaling in the commercial code ABAQUS. To avoid confusion, we propose the term adaptive mass scaling. The term selective mass scaling is usually employed for non-uniform scaling of the two or three-dimensional elemental mass matrix. This procedure allows the decrease in the maximum eigenfrequency with a lower amount of artificial mass than in the constant scaling procedures [START_REF] Olovsson | Selective mass scaling for thin walled structures modeled with tri-linear solid elements[END_REF]. These techniques are also known as mass penalty methods [START_REF] Macek | A mass penalty technique to control the critical time increment in explicit dynamic finite element analyses[END_REF][START_REF] Plecháč | Implicit Mass-matrix Penalization of Hamiltonian Dynamics with Application to Exact Sampling of Stiff Systems[END_REF][START_REF] Askes | Increasing the critical time step: micro-inertia, inertia penalties and mass scaling[END_REF][START_REF] Lombardo | Lumped mass finite element implementation of continuum theories with micro-inertia[END_REF]. A major drawback of selective mass scaling is that the lumping of the mass matrix is no more possible. The solution of the system in equation (5.8) requires a special solver. For shell elements, due to their low aspect ratio, this method was further developed in [START_REF] Cocchetti | Selective mass scaling and critical time-step estimate for explicit dynamics analyses with solid-shell elements[END_REF][START_REF] Pagani | Computationally efficient explicit nonlinear analyses using reduced integration-based solid-shell finite elements[END_REF][START_REF] Tkachuk | Local and global strategies for optimal selective mass scaling[END_REF][START_REF] Tkachuk | Local and global strategies for optimal selective mass scaling[END_REF].

Time scaling

Another option is the increase in the loading rate. In literature this process is known as time scaling or load factoring [START_REF] Prior | Applications of implicit and explicit finite element techniques to metal forming[END_REF]. In the following, the time scaling factor 𝛼 𝑇 is defined as the increase in the static loading rate

𝑣 𝑙𝑜𝑎𝑑 = 𝑣 𝑙𝑜𝑎𝑑 𝑠𝑡𝑎𝑡𝑖𝑐 𝛼 𝑇 ⇒ θ = θ𝑠𝑡𝑎𝑡𝑖𝑐 𝛼 𝑇 (5.14)
For rate dependent materials the material coefficients should be scaled according to the time scaling factor [START_REF] Lorenz | Recent Advances and New Developments in hot Forming Simulation with LS-DYNA[END_REF]. Together with mass scaling the total scaling is then defined as 𝛼 = 𝛼 𝑀 𝛼 𝑇 .

(5.15)

Control of the validity of the solution through the energy ratio

The scaling techniques lead to an increase in the kinetic energy, through the augmented rate and mass. If the scaling is too high, the kinetic energy becomes very high in relation to the internal energy. In this case the obtained displacement differs too much from the result of a static computation.

In the case no quasi-static solution exists for the validation, several error criteria were proposed to control the solution. The first two criteria are the control of the kinetic forces compared with external forces, 𝜖 1 , and the control of the ratio of kinetic and internal energy, 𝜖 2 . A more complex proposition considers the kinetic energy increase from the initial yielding 𝑡 𝑦𝑖𝑒𝑙𝑑 and the plastic dissipation 𝐸 𝑝𝑙𝑎𝑠𝑡𝑖𝑐 𝑖𝑛𝑡 [START_REF] Chung | On the dynamic effects of explicit FEM in sheet metal forming analysis[END_REF][START_REF] Kim | A comparative study of implicit and explicit FEM for the wrinkling prediction in the hydroforming process[END_REF],

𝜖 1 = 𝐹 𝑘𝑖𝑛 𝐹 𝑒𝑥𝑡 100, 𝜖 2 = 𝐸 𝑘𝑖𝑛 𝐸 𝑒𝑥𝑡 100, 𝜖 3 = 𝐸 𝑘𝑖𝑛 -𝐸 𝑘𝑖𝑛𝑡 𝑦𝑖𝑒𝑙𝑑 𝐸 𝑝𝑙𝑎𝑠𝑡𝑖𝑐 𝑖𝑛𝑡 100.
(5.16)

The critical value of the energy criterion is often defined in the range of 5 -10 %, which seems arbitrary [START_REF] Wong | Incremental forming of solid cylindrical components using flow forming principles[END_REF][START_REF] Wang | Investigation of material deformation in multi-pass conventional metal spinning[END_REF][START_REF] Burkhart | Finite element modeling mesh quality, energy balance and validation methods: A review with recommendations associated with the modeling of bone tissue[END_REF][START_REF] Gulavani | Explicit dynamic formulation to demonstrate compliance against quasi-static aircraft seat certification loads (CS25.561) -Part I: influence of time and mass scaling[END_REF][START_REF] Huang | Study on plastic deformation behavior of hot splitting spinning of TA15 titanium alloy[END_REF]Kim et al., 2014b;[START_REF] Li | A numerical study of the effects of roller paths on dimensional precision in die-less spinning of sheet metal[END_REF][START_REF] Natário | Web crippling failure using quasi-static FE models[END_REF].

In fact this is one of the problems of the method: how to assess a priori the quality of the prediction associated with the use of scaling for a given quantity i.e. , in this study, the occurrence of failure or the determination of the critical rotational speed.

Further recommendations

In fact, many authors provide additional criteria and recommendations:

• The limitation of the total scaled mass 𝑚 𝑡,𝑠𝑐𝑎𝑙 to 𝑚 𝑡,𝑠𝑐𝑎𝑙 /𝑚 𝑡 = 500 (Kim et al., 2014b). For a uniform mesh, this means 𝛼 𝑀𝑚𝑎𝑥 = 22.3. However, for heterogeneous meshes, the scaling can be higher, as small elements contribute little to the total mass.

• No applications of time scaling to damage models have been reported. However, the deletion of elements from the simulation lead to a rise of the energy relation above 30 % [START_REF] Bekker | Introducing objective power law rate dependence into a visco-elastic material model of bovine cortical bone[END_REF]).

• The application of time scaling should be limited to ten times the maximal period 𝑇 𝑚𝑎𝑥 [START_REF] Prior | Applications of implicit and explicit finite element techniques to metal forming[END_REF], because the minimal eigenfrequency 𝜔 𝑚𝑖𝑛 dominates the mechanical response, 𝛼 < 10 𝑇 𝑚𝑎𝑥 = 10 𝜔 𝑚𝑖𝑛 .

(5.17)

• In case of non-linearities, the scaling should often be adapted. In [START_REF] Mattiasson | Solution of quasi-static, force-driven problems by means of a dynamic-explicit approach and an adaptive loading procedure[END_REF] an adaptive scheme that controls the loading rate depending on the velocity of a control node was proposed.

Other recommendations are only qualitative:

• Oscillations in the global solution curves should be avoided [START_REF] Natário | Web crippling failure using quasi-static FE models[END_REF].

• Loading shocks should be avoided through the use of adapted non-linear load functions (ABAQUS, 2011).

• The velocity should be minimal at the end of the computation [START_REF] Gulavani | Explicit dynamic formulation to demonstrate compliance against quasi-static aircraft seat certification loads (CS25.561) -Part I: influence of time and mass scaling[END_REF].

Modeling of the rotating disk problem as mass-spring system

This analysis was performed in order to acquire a first impression whether it is better to perform mass or time scaling and to determine the expectable errors considering the geometric non-linearity of the disk burst and plasticity.

Static solution

A mass 𝑚 rotates around a fixed point with the rate θ, cf. figure 5.15. The mass is a kind of effective mass chosen equal to half the mass of the equivalent disk, as would be obtained by using a single bar finite element. The position of the mass initially situated at the position 𝑟 0 is actualized during the computation to 𝑟 0 + 𝑢, u being the displacement.

The non linear measure of strain as a function of the displacement is given by 𝜀 = 𝑢 𝑟 0 +𝑢 . In order to provide a first analytical estimation of the limit speed, the exponential part of the extended Voce law is neglected. The model then reduces to linear hardening. As shown in figure 5.15(b) the stress can in this case be described as the sum of the yield stress 𝜎 𝑦 and the contribution of the product of the tangent modulus 𝐶 𝑡𝑎𝑛 and the delta between total deformation 𝜀 𝑡 and the deformation at the yield stress 𝜀 𝑦 . The tangent modulus is given by 𝐶 𝑡𝑎𝑛 = 𝐶.𝐸 𝐶+𝐸 . This model is of course debatable but it allows for a simple introduction of an instability associated to geometrical non-linearities combined with plasticity. One assumes that the instability occurs for a stress exceeding the yield stress 𝜎 𝑦 . The section being considered as constant, the force balance can then be simplified to

𝜌𝑟 0 2 𝑟 0 (𝑟 0 + 𝑢) θ2 = 𝜎 𝑦 + 𝐶 𝑡𝑎𝑛 ( 𝑢 𝑟 0 + 𝑢 - 𝜎 𝑦 𝐸 ) = 𝜎 𝑦 ( 1 - 𝐶 𝑡𝑎𝑛 𝐸 ) + 𝐶 𝑡𝑎𝑛 𝑢 𝑟 0 + 𝑢 , (5.18)
where 𝜌𝑟 0 2 represents the density per unit length. This quadratic equation in 𝑟 0 + 𝑢 has no solution if

[ 𝜎 𝑦 ( 1 - 𝐶 𝑡𝑎𝑛 𝐸 ) + 𝐶 𝑡𝑎𝑛 ] 2 -2𝜌𝐶 𝑡𝑎𝑛 (𝑟 0 θ) 2 ≤ 0. (5.19)
From this it follows that the limit value θ𝑙𝑖𝑚 of the rotational speed is

θ𝑙𝑖𝑚 = 𝜎 𝑦 ( 1 -𝐶𝑡𝑎𝑛 𝐸 ) + 𝐶 𝑡𝑎𝑛 ) √ 2𝜌𝐶 𝑡𝑎𝑛 𝑟 0 .
(5.20)

In order to predict a limit value in this simple case, both spin softening effect and non-linear strain displacement relation have to be taken into account. Moreover, this equation shows that as expected an increase in the density or in the radius leads to a decrease in θ𝑙𝑖𝑚 .

For the material parameters 𝜌 = 4.5 𝑡𝑜𝑛 -9 /𝑚𝑚 3 𝐸 = 100000 𝑀 𝑃 𝑎, 𝜎 𝑦 = 800, 𝐶 = 200 𝑀 𝑃 𝑎 and the radius 𝑟 0 = 100 𝑚𝑚, the limit is θ𝑙𝑖𝑚 = 7449 𝑟𝑎𝑑/𝑠. For these values the speed of the elastic limit is indeed obtained as

θ𝑒𝑙𝑎𝑠 = √ 2𝜎 𝑦 ( 1 - 𝜎𝑦 𝐸 ) 𝜌𝑟 2 0 < θ𝑙𝑖𝑚 . (5.21)
The value of the displacement at the instability is

𝑢 𝑙𝑖𝑚 = 𝐶𝑟 0 2(𝜎 𝑦 ( 1 -𝐶𝑡𝑎𝑛 𝐸 ) + 𝐶 𝑡𝑎𝑛 ) -𝑟 0 .
(5.22)

For the extended voce law, introducing an exponential term with 𝑅 ∞ = 1300 𝑀 𝑃 𝑎, 𝛾 = 3, the critical value of the rotational speed, computed numerically slightly changes to θ𝑙𝑖𝑚 = 7302 𝑟𝑎𝑑/𝑠, cf. appendix K.

Explicit solution using mass and time scaling

In this simple case, the time step that still ensures the stability of the central difference scheme is about 𝛥𝑡 = 1.6 × 10 -5 𝑠. To examine the effect of mass and time scaling, different combinations of the factors 𝛼 𝑀 and 𝛼 𝑇 were examined. The implicit simulation was carried out with an initial time step of 𝛥𝑡 = 0.01 𝑠, which was lowered automatically towards the instability down to the unscaled explicit time step. Figure 5.16(a) presents the solution for different amounts of mass scaling and no time scaling. The values are normalized to the static limit speed 𝜃 𝑟𝑒𝑓 and the associated equivalent plastic strain 𝑝 𝑟𝑒𝑓 . The application of explicit dynamics allows for an extension of the calculation to reach a higher plastic equivalent strain, than is possible in the implicit solution. The instability speed, θ𝑖𝑛𝑠𝑡𝑎𝑏𝑖𝑙𝑖𝑡𝑦 , is defined by the moment the calculation fails to calculate for higher speeds. In addition, we use the definition of the stable phase, where the increase in the equivalent plastic strain remains relatively low, i.e. 𝜃 𝑟𝑒𝑓 θ < 0.98, and a transition regime.

Until the static limit speed is reached the simulation results remain identical. The unscaled calculation 𝛼 𝑀 = 1 leads to almost exactly the same limit speed as the static calculation. The application of more mass scaling leads to a delay of the instability speed.

The limit speed error indicator is defined as

𝜖 𝜃 = ( θ𝑖𝑛𝑠𝑡𝑎𝑏𝑖𝑙𝑖𝑡𝑦 θ𝑟𝑒𝑓 - 1
) 100.

(5.23)

The application of mass and time scaling led to the same errors, regardless of the combination for the same total scaling 𝛼. Therefore figure 5.16(b) shows the limit speed error indicator in dependence of 𝛼. While the error remains initially low, it strongly increases for 𝛼 ≥ 100. As mass and time scaling have the same effect on the limit speed error it was expected that this should be the same for the kinetic energy criterion. The variation in the kinetic energy should follow, (5.24) where 𝑚 0 is the unscaled mass and 𝑣 0 the unscaled velocity. The same should also apply to the error criterion on the kinetic energy 𝜖 2 . Indeed, as visible in figure 5.17(a), the evolution of 𝜖 2 shows this behavior. The curve for 𝛼 𝑀 = 10 is equal to the curve of 𝛼 𝑇 = 10 and the one of 𝛼 𝑀 = 100 to the one obtained for 𝛼 𝑇 = 100. The curves are superposed. An augmentation of 𝛼 by the factor 10 leads to a rise of the average level of the error in the stable phase by the factor 100.

𝐸 𝑘𝑖𝑛𝛼 ≈ 1 2 𝑚 0 𝛽 𝑀 𝑆 𝑣 2 0 𝛼 2 𝑇 = 1 2 𝑚𝛼 2 𝑀 𝑣 2 𝛼 2 𝑇 ,
Considering the magnitude all calculations remain below an error of 5 % (dashed line) up to the static limit speed. One can note that the increase in the ratio starts for all scaling factors in the transition regime.

The evolution of the equivalent plastic strain rate was observed, cf. figure 5.17(b). For the unscaled calculation, the average strain rate level in the stable part of the calculation is low, for instance at θ/ θ = 0.95, ṗ ≈ 0.3 1/𝑠, compared with the one obtained at the moment the curve reaches the reference equivalent plastic strain 𝑝 = 𝑝 𝑟𝑒𝑓 , ṗ ≈ 10 1/𝑠. At the limit speed the rate is even higher.

In the stable phase, the time scaling leads to a proportional increase in the equivalent plastic strain rate and a change of the wavelength of oscillations. In addition, the scaling results in an increase in the strain rate in the transition phase, even though without a clear relation.

In contrast, for mass scaling, no increase in the strain rate level was observed in the stable regime. However, the rate which is eventually obtained at the instability speed, depends on the mass scaling. The rule applies that mass scaling leads to a decrease in the final rate. As for the time scaling, no unique relation was found.

As follows from the study of the mass spring system, only the total scaling seems to determinate the level of the errors. However, the time scaling already leads to a change of the rate in the stable phase. In the transition regime up to the instability speed mass scaling and time scaling have opposite effects on the equivalent plastic strain rate. Time scaling increases the rate and mass scaling decreases the rate. However, the relationship between the applied factors could not be established. We will discuss the relevance of this for the calculation of rupture later.

Application of scaling to an axisymmetric disk model with tensile loading

To create an equivalent load case, while using similar structures as they are used in the industrial application, a disk was loaded under tension. The displacement was applied on the external radius as depicted in figure 5.18. To confirm the reported limit of 5 -10 % for the criteria and to reveal possible decrease in simulation time, the model was calculated with implicit statics and explicit dynamics. Different mass and time scaling combinations were applied for the latter. The implicit step size was initially set to 𝛥𝑡 = 0.1 𝑠 and allowed to be automatically reduced to 𝛥𝑡 = 10 -6 𝑠 close to the unscaled explicit step, 𝛥𝑡 = 7.1 × 10 -7 𝑠.

Instead of the limit speed error, a strain based error can be defined. The elemental strain error between the equivalent plastic strain in the implicit simulation, 𝑝 𝑒𝑙 𝑖𝑚𝑝 , and the one in the explicit simulation, 𝑝 𝑒𝑥𝑝 , is calculated as

𝜖 𝑒𝑙 𝑝 = |𝑝 𝑒𝑙 𝑖𝑚𝑝 -𝑝 𝑒𝑙 𝑒𝑥𝑝 | 𝑝 𝑒𝑙 𝑖𝑚𝑝 100.
(5.25)

To provide a total measure of the quality of the scaling, which weights higher strains more than lower strains, a global error indicator is calculated as

ε𝑝 = 1 𝑛 ∑ |𝑝 𝑒𝑙 𝑖𝑚𝑝 -𝑝 𝑒𝑙 𝑒𝑥𝑝 | 𝑝 𝑚𝑎𝑥 𝑖𝑚𝑝 100,
(5.26)

where the maximum plastic equivalent strain that occurs in the implicit calculation 𝑝 𝑖𝑚𝑝𝑚𝑎𝑥 is used and 𝑛 are the number of elements. The amount of displacement applied was chosen in order for the maximum equivalent plastic strain in the implicit solution to reach 𝑝 𝑖𝑚𝑝𝑚𝑎𝑥 = 0.5. The plasticity was limited to the outer region of the disk, therefore figure 5.18 only shows a zoom of the external area of the disk.

Table 5.4 presents the used scaling factors and the obtained results. One combination for 𝛼 = 10 3 was tested and several for a total scaling of 𝛼 = 10 4 . The decrease in simulation time was calculated with respect to the simulation time of the implicit simulation 𝑡 𝑖𝑚𝑝 = 118 𝑠. The solution for 𝛼 = 10 3 takes twice the time of the implicit solution. For the faster scaling 𝛼 = 10 4 , the possible decreases of time are all above 300 % resulting in simulation times of less than a minute. The effect of other factors cannot be excluded at so small simulation times, hence it is difficult to compare the calculations against each other.

The scaling 𝛼 = 10 3 leads to a very low value of the error indicator ε𝑝 = 1.3 × 10 -2 %. For the three higher combinations with 𝛼 = 10 4 , the mean error remains below ε𝑝 = 0.8 %. However, there is a small advantage for the two scaling 𝛼 𝑀 = 10, 𝛼 𝑇 = 1000, 𝛼 𝑀 = 10, 𝛼 𝑇 = 1000 compared with 𝛼 𝑀 = 100, 𝛼 𝑇 = 100. The combination of equal factors leads to higher errors.

Regarding the kinetic error criterion 𝜖 2 , there is a clear advantage in the use of mass scaling 𝜖 2 = 2.3 × 10 -4 for 𝛼 𝑀 = 1000, 𝛼 𝑇 = 10 compared to 𝜖 2 = 2.3 × 10 -2 for 𝛼 𝑀 = 10, 𝛼 𝑇 = 1000. Nevertheless, all scaling combinations lead too a value of less than 5 % and should therefore be valid concerning the literature recommendation. However, as shown in the following there are very high local errors. The almost identical mean errors were confirmed by the analysis of the fields of the elemental error, which showed very little difference for all three scaling combinations for 𝛼 = 10 4 . Therefore, figure 5.19 only shows the variation of the elemental strain error for the scaling 𝛼 = 10 3 and one for 𝛼 = 10 4 .

The error remains below 1 % for the scaling of 𝛼 = 10 3 . Existence of very low plastic strains leads to a higher error, e.g. on the left side of the examined area. In contrast, for the higher scaling 𝛼 = 10 4 the local errors can be very high.

It needs to be remarked that the high errors do not always coincide with the location of the highest equivalent plastic strain values in the implicit solution, e.g. the lower edge of the neck of the disk, cf. figure 5.18. However, for the upper edge, error values close to 10 % exist in regions with very large equivalent plastic strains. The calculation of failure with a strain based failure criterion would then be significantly impacted. Little can be said of the relation between the change in scaling and the variation in the elemental strain error, as the spatial distribution of the errors changed and not only the magnitude. These results lead us to the conclusion that asymmetric scaling combinations perform marginally better and should be preferred for plastic simulations. As for the amount of scaling that can be used, one has to assume that 𝛼 = 10 4 is an upper limit, while scaling up to at least 𝛼 = 1000 can still produce acceptable results. For the lower scaling a decrease in simulation time was not obtained. However, as demonstrated later this is due to the large steps of the implicit calculation, which are not possible if regularization is applied.

Finally, the use of a limit of 5% for the simple energy criterion 𝜖 2 does not seem justified, when comparing to the large local strain errors. A next step would be the examination of the more elaborated error criterion 𝜖 3 .

Towards the quasi-static mesh-independent disk fracture calculation General considerations

The above considerations for scaling techniques should in general also apply to an axisymmetric disk model under centrifugal loading. Let us recall two important points. First, the use of unscaled explicit dynamics does not fundamentally increase the limit speed for the mass-spring model and the same is presumed for a larger model. Second, a strong variation of the limit speed was already observed for a scaling of 𝛼 = 100 in the case of the mass-spring model.

Three different cases for the burst calculation with local fracture must be distinguished. In the first case, the disk burst is solely determined by a criterion, which is applied in post-processing, i.e. disk burst equals failure of the first element. Then the constitutive relation remains elasto-plastic. Such criteria are for example the once defined in chapter 4. At large plastic strains the implicit computation can become problematic. However, if the criteria are not attained, a switch to explicit dynamics as described below can be applied.

The second and third case introduce damage into the constitutive model. The characteristic difference is the duration of the damage phase with respect to the failure strain. The second case is defined by an activation of damage at or very close to the final failure strain, damage is brittle. The third case defines a progressively growing damage, activated at for example half of the final failure strain. Failure of the disk occurs at the moment the disk breaks apart in several pieces, i.e. a continuous line of failed elements traverses the model. Due to the catastrophic propagation of any crack in the case of the high rotational loads, the second case only extends the first by the prediction of the crack pattern in the disk.

In both damage models, the regularization is necessary for the simulation of structures. If the bounded rate model is to be used, small time steps are required from the moment of initialization of rupture in the structure. With this small time steps, the time spent for the simulation depends mainly on the simulated time duration between the rupture initiation and the complete fracture of the disk. For the second failure case, this duration is most likely very short, due to the catastrophic propagation of a crack.

The constitutive models are in this cases given by the bounded rate model with damage as described above. For the proposed formulation of the damage law, the two cases are obtained by using different values for damage threshold more or less close to critical equivalent plastic strain.

The robust solution in all three cases can be assured with the same procedure. All models can be initialized with an implicit calculation for the stable phase and a switch to explicit dynamics just before rupture starts and/or non-convergence occurs. As explained above the use of scaling allows a decrease in the simulation time but introduces errors if applied in the transition regime. Therefore, the amount of scaling should be limited, if possible.

However, if simulation times are too long, this is not possible and the choice of the scaling parameter combination becomes important. Furthermore, an optimal procedure would decrease the scaling with respect to the simulation time in order to avoid the shift of the limit speed through the scaling. For the time scaling this is possible by a proper choice of the loading rate. However, none of the tested codes (ABAQUS/Explicit, LS-DYNA) has the option to reduce the mass scaling during the calculation. This means not only a switch from implicit to explicit is necessary, but several explicit simulations have to be executed consecutively. In consequence, a substantial amount of the calculation time would be used for the initialization and control of the explicit simulations.

Furthermore, the definition of a criterion that defines the moment of the switch from explicit to explicit simulation and how to realize such a scheme in an industrial software package requires further studies. At present time no technique exists to perform such an operation. Therefore, the following simulations were restricted to the switch from implicit to explicit.

Use of the bounded rate model

Besides the restriction of the time step, the use of the bounded rate model adds another requirement to the simulation. The mass scaling modifies the wave speed through the relation

𝑐 = √ 𝐸 𝜌 1 𝛼 𝑀 .
(5.27) Furthermore, the size of the localization zone depends proportional on the wave speed as was demonstrated in (Court, 2006),

𝑙 𝑙𝑜𝑐 = 2𝑐 𝜏 𝑐 𝑎 𝜎 𝑦 𝐸 + 𝐶 √ 𝐸(𝜀 + 𝑝 𝑐 -2𝑝) (𝐸 + 𝐶)(𝑝 𝑐 -𝑝 0 ) , (5.28)
where, 𝜀 is the total strain. For a given mesh size the localization zone should not change for two different scaling factors. As follows from the above equations, the critical time can be increased by 𝛼 𝑀 to guarantee this condition,

𝜏 𝑐 = 𝜏 𝑐 0 𝛼 𝑀 , (5.29)
where 𝜏 𝑐 0 is the unscaled critical time, as measured in dynamic experiments. Further complications arise from the fact that for a heterogeneous mesh the mass scaling factor 𝛼 𝑒𝑙 𝑀 is not constant. In addition, the application of mass scaling leads to decrease in the maximum rate. In contrast, the use of time scaling increases the maximum rate. Further studies need to address how a change in 𝜏 𝑐 influences these results.

At the time of writing of this thesis, there were still some issues in numerical integration 𝛼 𝑀 = 5, 𝛼 𝑇 = 500 𝜖 2 = 2 %. One has to note that for 𝛼 𝑀 = 25, 𝛼 𝑇 = 100 the kinetic error criterion reaches 𝜖 2 = 2.3 %.

Summary

This chapter explores the numerical possibilities to provide a mesh-independent fracture prediction of the burst speed. The application of a damage model, progressive or brittle, to finite elements leads to mesh dependent solutions. An overview on the techniques guaranteeing mesh-independent results is presented, namely, non-local and bounded rate models.

The bounded plasticity rate model limits the rate of the equivalent plastic strain. For dynamic tensile testing a comparison of mesh independent results with experiments exists in the literature. The influence of the centrifugal force was evaluated by studying an axisymmetric mesh of a simple uniform thickness in a highly dynamic simulation. However, this case always leads to the same kind of damage distribution regardless of the mesh size. Differences between a classical non-limited damage model and the bounded rate model are minimal concerning the damage distribution. The uniform thickness associated with the axisymmetric computation and the centrifugal loading force the onset of damage to the center of the disk and always lead to the same damage gradient regardless of the mesh size. In consequence, the case is not very interesting to study mesh dependence.

Therefore, a disk with a notch was studied. The strong localization made it necessary to use element deletion in order to calculate a mesh with finer resolution up to final fracture of the disk. The application of the bounded rate model leads to mesh-independent simulations, i.e. the size of deleted and strongly damaged elements is constant.

In order to apply the bounded rate model, time steps below the critical time of the model are necessary. This conflicts with the times of the burst test. If in addition explicit dynamics is used to avoid convergence problems, the time steps are even smaller and it becomes impossible to perform simulations in realistic times.

A review on mass and time scaling techniques is provided. These techniques allow to accelerate the computation in the case of explicit dynamics.

The effects of these acceleration techniques are first studied in detail for a mass-spring system and then for the simplified disk under tension. Three main points need to be remarked. First, the use of scaling in general leads to a shift of the instability speed towards higher speeds in the case of centrifugal loading. Second, time scaling leads to an increase in the rates in the stable part of the calculation, and in a more complex manner also to an increase in the achievable maximum rate. Mass scaling only leads to a decrease in the maximum rate. Third, it seems that high time scaling should be avoided together with high mass scaling.

With regard to the definition of a validation criterion, the tested criterion of the ratio of kinetic to internal energy predicts that all performed simulations are valid. Considering the existing errors in some of the simulations this does not seem correct.

Moreover, in the case of the bounded rate model as the critical time should be changed according to the mass scaling in the element. Finally, an examination of a classical damage model on the disk with a notch revealed that for quasi-static failure the expectable limit errors are small and the possible decrease in simulation time is large.

CHAPTER 6 Conclusion

The principal axis of this study was the identification of the material behavior of TA6V and U500 up to fracture under representative state of stresses with respect to the application in turbomachinery disks and the identification of the local condition at failure. For this, an integrated digital image correlation approach was employed. The other thematic axis of this thesis was the objective and robust prediction of rupture using the damage mechanics in the framework of explicit simulations for quasi-static loading and the study of the influence on the prediction of various scaling techniques.

Stress state sensitivity of the materials

The experimental campaign on the materials TA6V and U500 was motivated by earlier evidence of the stress state influence on plasticity and failure of a nickel alloy also used for disks. It was therefore decided to conduct, on dedicated test, a more detailed analysis regarding plasticity and failure of TA6V and U500. One of the main pillars of this approach was the application of the integrated digital image correlation technique. As this is a rather recent method, much of the work was dedicated to the methodological development.

The triaxiality state of the disk was analyzed. The main range of interest is between 𝜂 = 0.4 and 𝜂 = 0.75. Two different notched samples were designed, i.e. a thin sample with a range of triaxiality between 0.33 < 𝜂 < 0.5 and a thick sample with a range between 0.33 < 𝜂 < 0.8. The geometries were tested under quasi-static uniaxial tension.

The combination of large strains and notches led to frequent cracking of the paint used for the DIC speckle patterns. However, I-DIC, which can be seen as highly regularizing, grants a considerable robustness to the identification of the displacement fields and allows obtaining locally continuous strain fields even at the end of the test in presence of local cracks in the speckle pattern.

The use of I-DIC provides access to data that are not accessible for standard strain measurements. I-DIC allows the quality of a constitutive law and the model of the test to be checked against the recorded evolution of the speckle pattern and force data. The importance of using 3D models even for thin samples was noticed. Their use allows for a better identification as can be seen from the associated residuals. On the tested samples, the equivalent plastic strain locally reaches levels of 30 % for U500 and almost of 40 % for TA6V, which is about three times higher than the levels achievable in uniaxial tensile tests on smooth samples. Up to these levels of strain, the analysis of the test allowed to validate positive strain hardening model of plasticity. For the thick samples, the maximum equivalent plastic strain reaches roughly half of these values.

The difference of the maximum strain for the two sample types gave rise to a detailed analysis of the failure conditions on some selected samples. In order to gain insight into the possible cause of failure, fractographies were performed. For TA6V, the fracture surface indicates that void nucleation and growth are the primary reason of the found fracture. However, no indication of macroscopic damage was found before the occurrence of failure. For U500, fracture corresponds to the typical brittle failure phenomena. The fact that thick samples with high levels of triaxiality failed for levels of the equivalent plastic strain about half the ones observed for the thin samples, which have lower levels of triaxiality, led to the investigation into possible combined plasticity -stress criteria. The close integration of the experiments and simulations allowed the analysis of deformation and stresses of the elements near the failure surface at the instant just before failure. This allowed the testing of different criteria with a large amount of data. For this, a procedure was developed to analyze I-DIC mechanical fields in combination with the image of the fractured sample. The elements from an idealized fracture surface are extracted according to the crack on the recorded image, assuming the extrusion of the crack through the sample. The data of the extracted elements can then be analyzed to construct a fracture locus and to examine failure criteria. Several propositions for failure criteria from the literature were tested. It appears that the maximum principal stress criterion fits very well the extruded data. This should be checked on other test cases, on sample geometries involving other and various states of stresses.

A main improvement of the method would be the combination of all recorded data from all four cameras into one large I-DIC analysis. These cameras recorded the speckle pattern evolution on all four sides of the sample during the test. This multi-I-DIC should first lead to an improvement of the ROI boundary conditions. These boundary conditions could be interpolated from the side surface data or, even better, added to the identification in order to reconstruct more reliable boundary conditions. Such a procedure should improve the quality of the identification especially for thick samples. Moreover, it would be interesting to test the possibility to identify a unique constitutive model for all samples in a simultaneous identification. Such a method would require non-obvious developments but would allow defining in a clear way the best possible model from all the available information. In addition, it would allow validating or invalidating the current method, which defines the identified model as an average of the obtained hardening curves of several samples.

A step forward would be the analysis of the experimental fracture surface and not only an idealized one. This should be feasible, once the described improvements considering the boundary conditions in I-DIC are carried out. Another improvement of the method would be the introduction of the obtained fracture criteria in FE simulations. The erosion of elements from the simulation could lead to a virtual fracture surface. With additional procedures this virtual fracture surfaces could be compared to the recorded images of the experimental fracture surface. One could also check whether once the failure initiates, it propagates in a catastrophic manner. Further, with an adapted I-DIC procedure, the crack can be recreated in the computer-generated images. Such a procedure would enable the correlation of the fracture image with the I-DIC approach. The identification of plasticity and fracture could be performed simultaneously.

Together with the other improvements for I-DIC, the model identification for plasticity and fracture would then be one "big data" identification. Of course a practical aspect would be to find a procedure which allows avoiding the cracking of the paint used for the DIC speckle patterns. This would allow a better identification of the constitutive model close to failure.

Robust rupture simulations

The introduction of a local rate independent damage, if used to predict the failure of a structure, leads to spurious mesh dependence. In addition, the computation of failure using implicit static computations is very difficult and often leads to problems of convergence even when using path following techniques. A possible alternative is proposed in this thesis. It concerns the use of a bounded rate model combined with scaling techniques to simulate the failure of a quasi-static test in an objective manner. One of the points of interest is that such a model does not require complex developments, as necessary for non-local models, because it only implies local modifications of the constitutive law.

The bounded plasticity rate model limits the rate of the equivalent plastic strain. As originally only used for transient dynamics, the model was first tested for a dynamic loading of the disk. Two geometries were studied, namely, a disk of uniform cross-section and a disk with a notch. The disk of uniform cross-section induces the deterioration in the center from where it progressively spreads. This case appears not that relevant to study problems of localization. In contrast, for the disk with a notch, localization is effective and mesh dependence of the solution was encountered when using a classical damage model. The application of the bounded rate model leads to mesh-independent simulations; i.e. the size of deleted and strongly damaged elements is constant.

In order to apply the bounded rate model, the use of time steps below the stable time of the model is often necessary. For burst experiments this leads to simulations impossible to perform in realistic times. The idea was therefore to study, as often done for other applications, the effects of mass and time scaling techniques for the acceleration of explicit calculations.

Three main points need to be remarked if centrifugal loads are applied. First, the use of scaling in general leads to a shift of the predicted instability speed towards higher speeds. This effect should of course be controlled and limited. Second, time scaling leads to an increase in the rates and mass scaling leads to a decrease in the maximum rate. The relationship between these observations is complex. Third, it seems that high time scaling factors should be avoided in combination with high mass scaling factors.

If no implicit calculation exists a criterion is necessary to validate the procedure. The tested criterion of the ratio of kinetic to internal energy does not seem appropriate or at least not very sensitive to the occurring effects; the provided literature value of 5 % seems too high.

Considering these findings several conclusions can be drawn regarding the prediction of failure in the case of centrifugal loading. Due to problem of convergence in an implicit computation, the application of explicit simulations does allow the calculation of higher strains at the instability speed. Without scaling it does not change the predicted burst instability speed.

If brittle materials are considered and one is only interested in the burst speed, the application of a criterion (brittle damage model) is sufficient. Disk failure is almost immediate after the failure of the first element and up to this point static simulations can be used. This changes if the subsequent crack path is of interest for the design and a brittle damage model is applied. Similar to the case of a progressive damage model, regularization is necessary. Depending on the actual case of geometry and material the resulting simulation times can be improved with scaling techniques.

For a progressive damage model, regularization is obligatory and leads to considerably longer simulation times using the bounded rate approach together with implicit solution techniques and compared with the solution time of mesh-dependent static simulations. However, if explicit calculations are used, scaling techniques allow a decrease in the simulation time. The expected decrease in simulation times is large and the limit speed error remains less than 0.1 %.

In a future study the complete interaction of the bounded rate model and the scaling techniques should be studied. Once this is done the quasi-static burst prediction should be possible as well as the calculation of quasi-static rupture in tensile experiments.

No book can ever be finished. While working on it we learn just enough to find it immature the moment we turn away from it. (C.4)

The evaluation can be based on:

Hill criterion

The relationship between yield ratios and the coefficients of the

Titre:

  Vers la prévision de la rupture de disque de turbomachine : apport de la corrélation d'image numérique Mots clés : Turbomachines, DIC, Rupture, Triaxialité, Dépendance du maillage, Mass scaling

1. 1

 1 Schematic view of a turboshaft motor . . . . . . . . . . . . . . . . . . . . . 2.1 Material selection diagram for strength-temperature . . . . . . . . . . . . . 2.2 Microstructure of studied alloys . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Schematic stress-strain diagram of a brittle and a ductile material . . . . . 2.4 Disk fragments after burst . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Geometric non-linearity in disk simulations . . . . . . . . . . . . . . . . . . 2.6 Comparison of bust speed calculations and experiments for two different disks 2.7 Classification of fracture surfaces . . . . . . . . . . . . . . . . . . . . . . . . 2.8 The three stages of ductile rupture . . . . . . . . . . . . . . . . . . . . . . . 2.9 Formulations of the triaxiality function . . . . . . . . . . . . . . . . . . . . . 2.10 Selection diagram for samples with respect to triaxiality and Lode angle . . 2.11 Different types of fracture loci . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Explanation of correlation residuals . . . . . . . . . . . . . . . . . . . . . . . 2.13 Explanation of local and global digital image correlation . . . . . . . . . . . 3.1 Geometry, equivalent plastic strain field and triaxiality field for two disks . 3.2 Quantitative analysis of triaxiality distribution in disks . . . . . . . . . . . . 3.3 Sample geometries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Triaxiality analysis for TA6V . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Triaxiality analysis for U500 . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Uniaxial yield curves of TA6V and U500 . . . . . . . . . . . . . . . . . . . . 3.7 Experimental setup of cameras . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Examples of the Region of Interest (ROI) on the flat side of the samples . . 3.9 Image selection from all recorded images for I-DIC analysis . . . . . . . . . 3.10 Examples of meshes used for digital image correlation . . . . . . . . . . . . 3.11 Examples of the three different surface meshes for I-DIC . . . . . . . . . . . 3.12 Material orientation of samples . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Force-displacement curve as measured by the test machine . . . . . . . . . . 3.14 Definition of macro and meso measurement scales for DIC . . . . . . . . . . 3.15 Longitudinal displacement and strain fields for thin titanium alloy sample obtained from face camera images with DIC . . . . . . . . . . . . . . . . . . 3.16 Additional displacement fields for thin titanium alloy sample at final instant v 3.17 Normalized gray level residual images . . . . . . . . . . . . . . . . . . . . . 49 3.18 Dimensionless correlation residuals based on DIC . . . . . . . . . . . . . . . 49 3.19 Comparison between measured and identified I-DIC load levels . . . . . . . 50 3.20 Longitudinal displacement and strain fields on the surface for a thin titanium alloy sample obtained via I-DIC . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.21 Equivalent plastic strain fields for two different locations in plastic regime . 52 3.22 Von Mises' equivalent stress fields . . . . . . . . . . . . . . . . . . . . . . . . 52 3.23 Triaxiality fields for two different locations in elastic regime . . . . . . . . . 53 3.24 Triaxiality fields for two different locations in plastic regime . . . . . . . . . 53 3.25 Triaxiality history for four different points for TA6V . . . . . . . . . . . . . 54 3.26 Normalized gray level residual fields of TA6V thin sample for I-DIC . . . . 55 3.27 Variation of different residuals with the applied deformation . . . . . . . . . 56 3.28 Variation of the identification through mesh dimensionality and resolution . 57 3.29 Effect of sample alignment on identification . . . . . . . . . . . . . . . . . . 57 3.30 Comparison of stres/strain curves for different hardening models . . . . . . 58 3.31 Average hardening curves for TA6V . . . . . . . . . . . . . . . . . . . . . . 59 3.32 I-DIC identification of the Hosford and Hill criteria for TA6V . . . . . . . . 60 3.33 DIC fields at final test instant for a thick titanium alloy sample at final instant 61 3.34 Force -elongation curve for a thick sample for its identification and the verification for TA6V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.35 I-DIC results for thin samples made of U500 . . . . . . . . . . . . . . . . . . 62 3.36 I-DIC identification of the Hosford and Hill criteria for U500 . . . . . . . . 64 3.37 Mechanical fields for U500 in mid-section plane . . . . . . . . . . . . . . . . 64 3.38 Triaxiality history for four different points for U500 . . . . . . . . . . . . . . 65 3.39 Force -elongation curve for a thick sample for its identification and the verification for U500 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.40 Verification of boundary condition extrusion hypothesis . . . . . . . . . . . 67 4.1 Fractographies of notched samples made of TA6V . . . . . . . . . . . . . . . 70 4.2 Results from loading-unloading tests . . . . . . . . . . . . . . . . . . . . . . 70 4.3 Selection procedure for the fracture locus for a sample made of TA6V . . . 72 4.4 Superimposition of shifted fracture image and critical elements for the thick sample made of TA6V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.5 Proposition from literature for the fracture locus for elements in the critical section for TA6V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.6 Fracture locus for critical damage and triaxiality for elements in critical section for TA6V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.7 Fracture loci based on other mechanical quantities for elements in the critical section for TA6V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.8 Verification of the fracture for all elements in the last step for TA6V . . . . 76 4.9 Fractographies of notched samples made of U500 . . . . . . . . . . . . . . . 77 4.10 Selections of critical elements from different samples made of U500 . . . . . 77 4.11 Proposition from literature for fracture locus for elements in the critical section for U500 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.12 Fracture locus for critical damage and triaxiality for elements in critical section for U500 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Fracture loci based on other mechanical quantities for elements in the critical section for U500 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Application of the bounded plasticity rate model to a 2D tensile simulation 5.2 Comparison of dynamic experiment and mesh independent simulations . . . 5.3 Comparison of damage formulation with effective and with macroscopic equivalent plastic strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Model of a uniform disk with coarse mesh using axisymmetry and symmetry along the radial axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Mesh convergence analysis for an uniform axisymmetric disk model in plasticity 5.6 Evolution of the equivalent plastic strain in the most loaded element . . . . 5.7 Evolution of the damage field for the classical model . . . . . . . . . . . . . 5.8 Evolution of the damage field for the bounded rate model . . . . . . . . . . 5.9 Model of a notched disk with coarse mesh using axisymmetry and symmetry along the radial axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Variation of the equivalent plastic strain rate and damage fields . . . . . . . 5.11 Evolution of the damage field for the classical model . . . . . . . . . . . . . 5.12 Evolution of the damage field for the bounded rate model . . . . . . . . . . 5.13 Evolution of the equivalent plastic strain rate field for the classical model . 5.14 Evolution of the equivalent plastic strain rate field for the bounded rate model 5.15 Simplified model for the disk burst analysis . . . . . . . . . . . . . . . . . . 5.16 Effect of mass and time scaling on the limit speed of the mass spring system 5.17 Evolution of energy ratio and equivalent plastic strain rate for mass and time scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 Equivalent plastic strain field of disk under tension . . . . . . . . . . . . . . 5.19 Equivalent plastic strain error for the disk loaded under tension . . . . . . . A.1 Graphical representation of equivalent stress criterion . . . . . . . . . . . . D.1 Spatial analysis of triaxiality distribution in disks . . . . . . . . . . . . . . . D.2 Spatial analysis for most dominant triaxiality for samples made of TA6V . . D.3 Spatial analysis for most dominant triaxiality for samples made of U500 . . E.1 Variation of identification with regard to maximum strain in analysis . . . . F.1 I-DIC results for thick samples made of TA6V . . . . . . . . . . . . . . . . . F.2 I-DIC results for thick samples made of U500 . . . . . . . . . . . . . . . . . G.1 Comparison of identified yield curves for TA6V . . . . . . . . . . . . . . . . G.2 Comparison of identified yield curves for U500 . . . . . . . . . . . . . . . . I.1 Verification of all data points for fracture criteria . . . . . . . . . . . . . . . List of Tables 2.1 Alloy composition in weight percent . . . . . . . . . . . . . . . . . . . . . . 2.2 Triaxiality ranges for uniaxial sample geometries . . . . . . . . . . . . . . . 3.1 Selected models and identifiable variables . . . . . . . . . . . . . . . . . . . 3.2 Global identification residual for different mesh analyses and corresponding equivalent plastic strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Global identification residual and equivalent plastic strain for different hardening law analyses via 3D-I-DIC-f for TA6V . . . . . . . . . . . . . . . 3.4 Comparison between the von Mises, Hosford and Hill equivalent stress criterion for TA6V . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Hill criterion anisotropy coefficients for TA6V . . . . . . . . . . . . . . . . . 3.6 Global identification residual and equivalent plastic strain for different hardening law via 3D-I-DIC-f for U500 . . . . . . . . . . . . . . . . . . . . . 3.7 Comparison between the von Mises, Hosford and Hill equivalent stress criterion for U500 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Hill criterion anisotropy coefficients for U500 . . . . . . . . . . . . . . . . . 5.1 Mass lumping techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Characteristic length estimations for the main element technologies . . . . . 5.3 Wave speeds with the main model assumptions . . . . . . . . . . . . . . . . 5.4 Comparison of mass and time scaling results for the disk model with tensile loading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Comparison of explicit and implicit simulations for a damage model . . . . C.1 Average triaxialities in sample geometries calculated in different ways . . . . F.1 Characteristic values of thick samples made of TA6V . . . . . . . . . . . . . F.2 Characteristic values of thick samples made of U500 . . . . . . . . . . . . . H.1 Mean coefficients of extended Voce law . . . . . . . . . . . . . . . . . . . . . H.2 Coefficients for von Mises criterion and extended Voce hardening used for comparison to other criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . H.3 Coefficients for Hosford criterion and extended Voce hardening identification147 H.4 Coefficients for Hill criterion and extended Voce hardening identification . .

Figure 1 . 1 :

 11 Figure 1.1: Schematic design of a turboshaft motor, after[START_REF] Mazière | Eclatement des disques de Turbomachines -Burst of turboengine disks[END_REF]: the combustion gas flows out of the combustion chamber in the high pressure turbine, which is part of the first shaft and drives the compressor. The power turbine stage of the second shaft transforms the remaining stream energy into mechanical energy. This shaft also ensures the connection to the gearbox.

Figure 2 . 1 :

 21 Figure 2.1: Material selection diagram for strength-temperature: titanium and nickel alloys combine high strength and high temperature[START_REF] Lovatt | Material selection and processing[END_REF].

Figure 2 . 2 :

 22 Figure 2.2: Microstructure of the studied alloys.

Figure 2 . 3 :

 23 Figure 2.3: Schematic stress-strain diagram of a brittle and a ductile material: the ductile material features strain softening and softening of the elastic modulus.

Figure 2 . 5 :

 25 Figure 2.5: Geometric non-linearity in disk simulations: (a) axisymmetric model of disk; (b) evolution of equivalent plastic strain near the point 𝑁 and the global deformation of this point; (c, d) equivalent plastic strain fields before and after the marked instant.

( a )

 a Figure 2.6: Comparison of bust speed calculations and experiments for two different disks: (a) prediction with Hosford's criterion (𝑛 = 80) provides better results; (b) in the second case, the deformation at the experimental burst speed is lower; the von Mises criterion provides a better prediction.

( a )

 a Figure 2.7: Classification of fracture surfaces.

Figure 2 . 8 :

 28 Figure 2.8: The three stages of ductile rupture: A: undamaged material, B: micro voids form around inclusions, C: voids enlarge and coalescent, D: a micro crack is established, after[START_REF] Flatten | Lokale und nicht-lokale Modellierung und Simulation thermomechanischer Lokalisierung mit Schädigung für metallische Werkstoffe unter Hochgeschwindigkeitsbeanspruchungen[END_REF].

( a )Figure 2 . 9 :

 a29 Figure 2.9: Formulations of the triaxiality function 𝑅 𝜈 .

Figure 2 .

 2 Figure 2.10 presents an overview of the most used sample geometries to examine the effect of triaxiality and Lode angle.The analysis of the stress state relied for a long time on analytical solutions. For axisymmetric samples, formulas predict the strain fields in elasticity[START_REF] Neuber | Theory of Stress Concentration for Shear-Strained Prismatical Bodies With Arbitrary Nonlinear Stress-Strain Law[END_REF]. To correct the measured data, other propositions use stress concentration factors[START_REF] Pilkey | Peterson's stress concentration factors[END_REF] or geometrical considerations to define a global triaxiality value[START_REF] Bridgman | Studies in large plastic flow and fracture with special emphasis on the effects of hydrostatic pressure[END_REF]. Using FE calculations, a correction of these latter formulas was proposed for round notched bars and flat grooved plates[START_REF] Bao | On the cut-off value of negative triaxiality for fracture[END_REF][START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF]

Figure 2 . 10 :

 210 Figure 2.10: Selection diagram for samples in dependence of triaxiality and Lode angle (Bai et al., 2008).

Figure 2 .

 2 Figure 2.11 presents different possibilities to visualize the dependence of the fracture criterion. The results in figure 2.11(a, b) were obtained for the aluminum alloy Al-2024-T351.For increasing triaxiality, the equivalent plastic strain at fracture decreases. Earlier findings reported a cut off value at 𝜂 = -1 3[START_REF] Bao | On the cut-off value of negative triaxiality for fracture[END_REF], below which fracture never occurs. In contrast, more recent findings provide a fracture locus in the range -0.7 ≤ 𝜂 ≤ 0.7[START_REF] Papasidero | Ductile fracture of aluminum 2024-T351 under proportional and non-proportional multi-axial loading: Bao-Wierzbicki results revisited[END_REF]. In figure2.11(b) the Lode angle dependence was included. The locus for the fracture strain 𝜀 𝑓 is then three-dimensional.Tensional tests of TA6V at low temperatures led to brittle fracture[START_REF] Diiorio | Plastic deformation, damage and rupture of PM Ti-6Al-4V at 20K under monotonic loading[END_REF]. Figure2.11(c) presents the identified fracture locus. The valid domain is limited by the maximum principal stress and the maximum equivalent plastic strain. The interpretation of this locus is more complicated, namely, fracture at low strains depends only on the principal stress, while for higher strains it depends also on the maximum equivalent plastic strain.The central issue about these loci is how the data are obtained from the experiments. For the presented triaxiality -equivalent plastic strain locus, the FE averaging method from equation (2.34) was used. For the creation of the maximum principal stress -equivalent plastic strain locus, the elemental FE values at the centre of the sample were used.

( a )

 a Figure 2.11: Different types of fracture loci: (a, b) proposition for ductile fracture of Al-2024-T351; (c) proposition for brittle fracture of TA6V.

Figure 2 . 12 :

 212 Figure 2.12: Explanation of correlation residuals: (a) gray level image 𝑓 with the marked outline of speckles; (b) deformed image 𝑔, with 10 % stretch; (c) residual image 𝜙(𝑥) for zero displacement, high residuals indicate the locations of the largest discrepancy; residuals are growing from left to right, indicating the stretch; (d) 𝜙(𝑥) for the displacement field of the 10% stretch, hence 𝜙(𝑥) = 0.

  Figure 2.13: Explanation of local and global digital image correlation: (a,b) local DIC evaluates the gray level balance for every window independently;(c, d) for global DIC the entire mesh is considered[START_REF] Mathieu | Analyse de la tenue mécanique d'un liner en titane : apport des mesures de champs cinématiques[END_REF].

  33

  Figure 3.1: Geometry, equivalent plastic strain field and triaxiality field for two disks: (a) disk geometry 𝐴 (non-drilled); (b) disk geometry 𝐵 (drilled); (c, d) equivalent plastic strain fields at final simulation step for both disks were normalized to the maximum value; areas of high equivalent plastic strain exist at points (1-6); (e, f) triaxiality fields at final simulation step for both disks; areas of high triaxiality at points (2, 5, 7-9).

Figure 3 . 2 :

 32 Figure 3.2: Quantitative analysis of triaxiality distribution in disks: all elements are plotted in a triaxiality -equivalent plastic strain space to identify maxima.

Figure 3 . 3 :

 33 Figure 3.3: Sample geometries (dimensions in mm): both samples are displayed with examples of the speckle pattern; the clamp area is kept paint-free to avoid slip.

  Figure 3.5: Triaxiality analysis for U500: (a, b) triaxiality field from thick and thin sample; (c, d) distribution of elements in triaxiality equivalent plastic strain space.

Figure 3 . 6 :

 36 Figure 3.6: Uniaxial yield curves of TA6V and U500: both curves are normalized to their respective maximum equivalent plastic strain and ultimate strength.

  Figure 3.7: Experimental setup of cameras: (a) CAD scene; (b) real test arrangement; the sample (S) is observed on its flat sides by camera (1) and (4); camera (2) and (5) observe the notched sides; camera (3) uses a special device to observe the same side as camera (4).

( a )

 a Figure 3.8: Examples of the Region of Interest (ROI) on the flat side of the samples: the gray-marked area depicts the ROI.

Figure 3 . 9 :

 39 Figure 3.9: Image selection from all recorded images for I-DIC analysis.

  Figure 3.10: Examples of meshes used for digital image correlation: (a) thin sample with 376 DOF and 314 elements; (b) thick sample with 216 DOF and 174 elements.

  3.2.4 Parameters for I-DIC Mesh resolution for I-DICThe mesh refinement determines the finite element discretization error and the correlation quality. To analyze the sensitivity of these effects, mesh size as well as dimensionality (2D vs. 3D simulations) were analyzed. This study addresses both the numerical convergence and the modification of the residuals. For higher resolutions, the correlation can be better. The mesh sizes are coarse, intermediate and fine (I-DIC-c, I-DIC-i and I-DIC-f), cf. figure3.11. The boundary nodes remain the same for all three meshes. The nodes of the coarse mesh are identical to the DIC mesh.

  Figure 3.11: Examples of the three different surface meshes for I-DIC: (a) 60-pixels coarse mesh with 157 elements; (b) 25-pixels intermediate mesh with 549 elements; (c) 10-pixels fine mesh with 2660 elements.

Figure 3 . 12 :

 312 Figure 3.12: Material orientation of samples.

Figure 3 .

 3 Figure 3.13: Force-displacement curve as measured by the testing machine for selected images.

Figure 3 .

 3 Figure 3.14(a) shows the strain field at the final instant measured with DIC. To exclude any displacement from the testing machine, one defines a virtual extensometer strain, 𝜀 = ū𝑦 /𝑙 0 , (3.3)

  figure 3.14(b).

Figure 3

 3 Figure 3.14: (a) Longitudinal dic strain field and definition of the measurement scales A: macroscopic and B: mesoscopic; (b) mean engineering axial stress versus mean engineering axial strain for two gauge lengths 𝑙 0 .

( a )

 a Figure 3.15: Longitudinal displacement and strain fields for thin titanium alloy sample obtained from face camera images with DIC.

  Figure 3.16: Additional displacement fields for thin titanium alloy sample at image 47, 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.072

( a )Figure 3 . 17 :

 a317 Figure 3.17: Normalized gray level residual images 𝜙 for DIC analysis.

Figure 3 . 18 :

 318 Figure 3.18: Dimensionless correlation residuals 𝜙 𝑐 (𝑡) based on DIC; instants 𝜀 𝑚𝑎𝑐𝑟𝑜 = {0.034, 0.061, 0.072} at image 28, 37 and 47 are marked with dash-dotted lines.

Figure 3 . 19 :

 319 Figure 3.19: Global I-DIC results, characteristic instants are marked with vertical lines: (a) comparison between measured and identified I-DIC load levels; (b) identified hardening curve.

( a )

 a Figure 3.20: Longitudinal displacement and strain fields on the surface for a thin titanium alloy sample obtained via I-DIC.

Figure 3 .

 3 Figure 3.21 displays the equivalent plastic strain field 𝑝(𝑥,𝑡) when 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.072.

Figure 3 . 21 :

 321 Figure 3.21: Equivalent plastic strain fields for two different locations at 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.072 (i.e. plastic regime).

Figure 3 . 22 :

 322 Figure 3.22: Von Mises' equivalent stress fields: (a) in the elastic regime; (b) for fully developed elasto-plasticity.

  Figure 3.23: Triaxiality fields for two different locations at 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.003 (i.e. elastic regime)

Figure 3 . 24 :

 324 Figure 3.24: Triaxiality fields for two different locations when 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.072 (i.e. plastic regime).

Figure 3 . 25 :

 325 Figure 3.25: Triaxiality history for four different points located either at the notch tip or in the center of the considered surfaces (i.e. external plane or midsection plane).

  Figure 3.26: (a-c) Normalized gray level residual images for the I-DIC-3D-f analysis; (d-f) difference of normalized gray level residuals between DIC and I-DIC-3D-f.

( a )

 a Figure 3.27: Variation of different residuals with the applied strain: (a) mean dimensionless correlation residuals 𝜙 𝑐 (𝑡) for entire analyzed sequence; (b) total identification 𝜒 𝑡𝑜𝑡 , correlation 𝜒 𝑐 and equilibrium residuals 𝜒 𝐹 for I-DIC-3D-f.

Figure 3 . 28 :

 328 Figure 3.28: Variation of the identification through mesh dimensionality and resolution.

Figure 3 . 29 :

 329 Figure 3.29: Effect of sampe alignment on identification results: (a) the two face cameras and the sample provide the same identified hardening curve; (b) sample were less well aligned, which resulted in different hardening laws.

Figure 3 . 30 :

 330 Figure 3.30: Comparison of stress/strain curves for different hardening models for TA6V.

( a )

 a Figure 3.31: (a) Averaging of hardening curves from I-DIC analysis for TA6V (L); (b) average hardening curves for both directions.

  Figure 3.32: I-DIC identification of the Hosford and Hill criteria: (a) Hardening curves for von Mises, Hosford and Hill criterion; (b) variation of normalized total residuals for exponent n; the minimum is at 𝑛 = 64.

Figure 3

 3 Figure 3.33: DIC fields for a thick titanium alloy sample at final instant: (a) longitudinal displacement on face fits well to recorded image; (b, c) longitudinal and transverse displacement fields of side camera show necking.

Figure 3 .

 3 Figure 3.34: Force -elongation curve for a thick sample for its identification and the verification with the average model of the thin samples for TA6V.

Figure 3 .

 3 Figure 3.35: I-DIC results for thin samples made of U500: (a) examples of force -elongation curves for two thick samples in both directions; (b) average hardening curves for both directions.

Figure 3 . 37 :

 337 Figure 3.36: I-DIC identification of the Hosford and Hill criteria: (a) Hardening curves for von Mises, Hosford and Hill criterion; (b) variation of normalized total residuals for exponent n; the minimum is at 𝑛 = 74.

Figure 3 . 38 :

 338 Figure 3.38: Triaxiality history for four different points located either at the notch tip or in the center of the considered surfaces (i.e. external plane or midsection plane) for U500.

Figure 3 .

 3 Figure 3.39: Force -elongation curve for a thick sample for its identification and the verification with the average model of the thin samples for U500.

( a )

 a Figure 3.40: Verification of boundary condition extrusion hypothesis: (a) the sample is modeled completely and as a model of the ROI; (b) force -elongation curves from the two models.

  69

Figure 4 . 1 :

 41 Figure 4.1: Fractographies of notched samples made of TA6V: (a) entire fracture surface; (b) a zoom reveals larger voids (1) and dimples (2).

  figure 4.2(a).

Figure 4 . 2 :

 42 Figure 4.2: Results from loading-unloading tests: (a) force -elongation curve and (b) I-DIC maximum equivalent plastic strain.

Figure

  Figure 4.2(b)shows the maximum equivalent plastic strain that was reached at every macroscopic elongation in I-DIC. The instants of the last two unloadings are marked in both diagrams with 𝜀 𝑚𝑎𝑐𝑟𝑜 𝐼 = 0.032 and 𝜀 𝑚𝑎𝑐𝑟𝑜 𝐼𝐼 = 0.053. At these instants, the equivalent plastic strain is 𝑝 𝐼 = 0.163 and 𝑝 𝐼𝐼 = 0.303. The absence of significant deviation of I-DIC from the experiment means that the global effect of a possible damage is not noticeable at 𝑝 𝐼𝐼 . Unfortunately, due to the strong heterogeneity of the test, this does not imply that no

( a )

 a Shifted fracture image (b) Superposition with last I-DIC solution mesh (c) Selected elements in I-DIC mesh (d) Superposition of extracted elements

Figure 4 . 3 :

 43 Figure 4.3: Selection procedure for the fracture locus for a thin sample made of TA6V: (a) the broken sample; (b) superposition of broken sample and I-DIC result before fracture; (c) selection of elements in the critical zone; (d) superposition of the equivalent plastic strain in the selected zone with the fracture image.

Figure 4 . 4 :

 44 Figure 4.4: Superposition of shifted fracture image and critical elements for the thick sample made of TA6V.

( a )Figure 4 . 5 :

 a45 Figure 4.5: Proposition from literature for the fracture locus for elements in the critical section for TA6V.

Figure 4 . 6 :

 46 Figure 4.6: Fracture locus for critical damage and triaxiality for elements in critical section for TA6V.

( a )Figure 4 . 7 :

 a47 Figure 4.7: Fracture loci based on other mechanical quantities for elements in critical section for TA6V.

Figure 4 . 8 :

 48 Figure 4.8: Verification of the fracture for all elements in the last step for TA6V.

Figure 4 . 9 :

 49 Figure 4.9: Fractographies of notched samples made of U500: (a) entire intergranular fracture with complex surface structure (1); (b) zoom with distortion bands (2) and microvoids on the grain boundary (3).

Figure 4 . 10 :

 410 Figure 4.10: Selections of critical elements from different samples made of U500.

( a )Figure 4 . 11 :

 a411 Figure 4.11: Proposition from literature for fracture locus for elements in the critical section for U500.

Figure 4 . 12 :

 412 Figure 4.12: Fracture locus for critical damage and triaxiality for elements in critical section for U500.

  (a) Von Mises stress -triaxiality for elements in critical section (b) Maximum principal stress -von Mises stress for elements in critical section (c) Maximum principal stress -triaxiality for elements in critical section (d) Verification of the fracture locus maximum principal stress -equivalent plastic strain for all elements

Figure 4 . 13 :

 413 Figure 4.13: Fracture loci based on other mechanical quantities for elements in the critical section for U500.

Figure 5 . 1 :

 51 Figure 5.1: Application of the bounded plasticity rate model to a 2D tensile simulation (Court, 2006).

Figure 5 . 2 :

 52 Figure 5.2: Mesh independent simulations for dynamic testing: (a) comparison of simulation results to experiment; (b) objectivity of the results with respect to the size of the mesh (Allix, 2012).

Figure 5 . 3 :

 53 Figure 5.3: Comparison of damage formulation with effective and with macroscopic equivalent plastic strain: the usage of the effective equivalent plastic strain can lead to very large total deformations.

Figure 5 . 4 :

 54 Figure 5.4: Model of a uniform disk with coarse mesh using axisymmetry and mid-plane symmetry.

  Figure 5.5: Mesh convergence analysis for an uniform axisymmetric disk model in plasticity: (a-c) the mesh is refined homogeneously; (d) the values of internal energies, displacement of the external radius and equivalent plastic strain in the lower left corner are compared to the solution of the finest mesh.

Figure 5 . 6 :

 56 Figure 5.6: Evolution of the equivalent plastic strain in the most loaded element.

Figure 5 . 7 :

 57 Figure 5.7: Evolution of the damage 𝐷 field for the classical model for (a-d) 𝑡 𝑒𝑙 = 1𝑚𝑚, (e-h) 𝑡 𝑒𝑙 = 0.5𝑚𝑚 and (i-l) 𝑡 𝑒𝑙 = 0.25𝑚𝑚 (only the center of the disk is shown).

Figure 5 . 8 :

 58 Figure 5.8: Evolution of the damage 𝐷 field for the bounded rate model for (a-d) 𝑡 𝑒𝑙 = 1𝑚𝑚, (e-h) 𝑡 𝑒𝑙 = 0.5𝑚𝑚 and (i-l) 𝑡 𝑒𝑙 = 0.25𝑚𝑚 (only the center of the disk is shown).

Figure 5 . 9 :

 59 Figure 5.9: Model of a notched disk with coarse mesh using axisymmetry and mid-plane symmetry.

( a )

 a Figure 5.10: Variation of the equivalent plastic strain rate ṗ and damage 𝐷 fields for two models at the instant the fine mesh computation fails.

  Figure 5.11: Evolution of the damage 𝐷 field for the classical model for (a-d) 𝑡 𝑒𝑙 = 1𝑚𝑚 and (e-h) 𝑡 𝑒𝑙 = 0.5𝑚𝑚.

  Figure 5.12: Evolution of the damage 𝐷 field for the bounded rate model for (a-d) 𝑡 𝑒𝑙 = 1𝑚𝑚 and (e-h) 𝑡 𝑒𝑙 = 0.5𝑚𝑚.

  Figure 5.14: Evolution of the equivalent plastic strain rate ṗ field for the bounded rate model for (a-d) 𝑡 𝑒𝑙 = 1𝑚𝑚 and (e-h) 𝑡 𝑒𝑙 = 0.5𝑚𝑚.

  Figure 5.15: Simplified model for the disk burst analysis.

( a )

 a Figure 5.16: Effect of mass and time scaling on the limit speed of the mass spring system: (a) the limit curves for 𝛼 𝑇 𝑆=1 ; (b) the increase in the instability error for the total scaling 𝛼.

( a )

 a Figure 5.17: (a) Evolution of the error criterion in kinetic energy: the curves with the same values for 𝛼 𝑀 and 𝛼 𝑇 are superposed; (b) evolution of the equivalent plastic strain rate varies with time scaling.

Figure 5 . 18 :

 518 Figure 5.18: Equivalent plastic strain field of disk under tension: the displacement is applied to the external radius nodes and leads only to plastic strain in the external region.

( a )Figure 5 . 19 :

 a519 Figure 5.19: Equivalent plastic strain error 𝜖 𝑒𝑙 𝑝 in % for the disk loaded under tension for two different total scaling.

  equivalent plastic strain is included 𝑝 𝑚𝑎𝑥 .

  Figure D.1: Spatial analysis of triaxiality reveals most dominant triaxiality value in term of a intensity function 𝛤 ; higher values are more important.

Figure D. 3 :

 3 Figure D.3: Spatial analysis for most dominant triaxiality for U500. For U500, significant amounts of very high strains exist for 𝜂 = 0.5, cf. figure D.3(a). For 𝑝 ′ = 0.4, the maximum triaxiality is 𝜂 = 0.7. For the thin samples, most part of the volume is in a state around 𝜂 = 0.4, cf. figure 3.5(b).

F

  Figure F.2: I-DIC results for thick samples made of U500: (a) examples of force -elongation curves for two thick samples in both directions; (b) average hardening curves for both directions.

Table 2 . 2 :

 22 Triaxiality ranges for uniaxial sample geometries.

	Type	𝜂 𝑎𝑣	Source
	Round, smooth	0.4	(Wierzbicki et al., 2005)
	Round, notch	0.49-0.95 (Wierzbicki et al., 2005; Bao et al., 2005)
	Flat, smooth	0.36-0.38

Table 3 . 1 :

 31 Selected models and identifiable variables.

	Model	Variables	Unknowns Equation
		Plasticity analysis
	Ludwik	𝐸, 𝜎 𝑦 , 𝐾, 𝑚	4	eq. (2.16)
	Voce	𝐸, 𝜎 𝑦 , 𝑅 ∞ , 𝛾	4	eq. (2.17)
	ext. Voce	𝐸, 𝜎 𝑦 , 𝑅 ∞ , 𝛾, 𝐶	5	eq. (2.18)
		Stress state analysis
	ext. Voce + D.Prager 𝐸, 𝜎 𝑦 , 𝑅 ∞ , 𝛾, 𝐶, 𝛽	6	eq. (2.18 & 2.20)
	ext. Voce + Hill	𝐸, 𝜎 𝑦 , 𝑅 ∞ , 𝛾, 𝐶, 𝑅 𝑖𝑗	10	eq. (2.18 & 2.21)
	ext. Voce + Hosford	𝐸, 𝜎 𝑦 , 𝑅 ∞ , 𝛾, 𝐶	5	eq. (2.18 & 2.19)

Table 3 . 2 :

 32 Global identification residual 𝜒 tot for different mesh analyses and corresponding equivalent plastic strain 𝑝 𝑚𝑎𝑥 .

Table 3 . 3 :

 33 Global identification residual 𝜒 tot and equivalent plastic strain for different hardening law analyses via 3D-I-DIC-f.

	Constitutive law 𝜒 tot 𝑝 max
	Ludwik	9.26 0.37
	Voce	9.23 0.39
	ext. Voce	9.13 0.37
	Elastic	339	-

Table 3 . 4 :

 34 Comparison between the von Mises, Hosford and Hill equivalent stress criterion with extended Voce law via 3D-I-DIC-i.

		𝜒 tot Relative gain 𝑝 𝑚𝑎𝑥 𝑅(𝑝 = 0.1) MPa
	von Mises 9.13	0	0.36	1147
	Hosford	8.62	0.056	0.35	1236
	Hill	9.06	0.008	0.29	1346

Table 3 .5: Hill

 3

criterion anisotropy coefficients.

Table 3 . 6 :

 36 Global identification residual 𝜒 tot and equivalent plastic strain for different hardening law via 3D-I-DIC-f.

Constitutive law 𝜒 tot 𝑝 max Ludwik 9.56 0.32 ext. Voce 9.44 0.31

Table 3 . 7 :

 37 Comparison between the von Mises, Hosford and Hill equivalent stress criteria with extended Voce law via 3D-I-DIC-i 𝜒 tot Relative gain 𝑝 𝑚𝑎𝑥 𝑅(𝑝 = 0.1) MPa

	von Mises 9.5	0	0.28	1239
	Hosford	9.21	0.031	0.29	1315
	Hill	8.77	0.076	0.34	1052

Table 3 . 8 :

 38 Hill criterion anisotropy coefficients.

Table 5 . 1 :

 51 Mass lumping techniques specified in

Table 5 . 2 :

 52 Characteristic length estimations for the main element technologies.For elastic materials with constant Young's modulus 𝐸, Poisson's coefficient 𝜈 and density 𝜌, the dilatational wave propagation speeds are summarized in table 5.3 in relation to the main model assumptions.

	Element technology	Characteristic length 𝐿 el
	2-node truss and beam	length
	3-node triangle and 4-node tetrahedron	minimum altitude
	4-node quadrilateral and 8-node hexaedron	𝑉 el /𝐴 el,max	a

Table 5 .3: Wave

 5 speeds with the main model assumptions.

	Model assumption	Dilatational wave speed 𝑐
	Truss and beam		√ 𝐸/𝜌
	Plane stress and shell	√	𝐸/(𝜌(1 -𝜈) 2 )
	Solid and plane strain	√ (𝐸(1 -𝜈))/(𝜌(1 + 𝜈)(1 -2𝜈))

Table 5 . 4 :

 54 Comparison of mass and time scaling results for the disk model with tensile loading.

	𝛼 𝑇 10 10 100 1000	𝛼 𝑀 100 10 3 𝛼 1000 10 4 100 10 4 10 10 4	Decrease (%) -50 440 360 380	ε𝑝 1.3 × 10 -2 4.1 × 10 -8 𝜖 2 (%) 0.7 2.3 × 10 -4 0.8 2.3 × 10 -2 0.7 2.3

Table 5 . 5 :

 55 Comparison of explicit and implicit simulations for a damage model

	𝛼 𝑇	𝛼 𝑀	𝛼	𝛥𝑡 (s)	θ𝑚𝑎𝑥 (rad/s) 𝜖 2 (%) 𝑡 𝑟𝑒𝑎𝑙 (h)
	25	100	2500	9 × 10 -06	8801	2.30	0h06
	500	5	2500 4.5 × 10 -07	8801	2	0h04
	5	500	2500 4.5 × 10 -05	8801	1.6	0h07
	25	200	5000 1.8 × 10 -05	8807	3.3	0h04
	25	IMPLICIT	25	10 -5	8778	0.03	9h40
	25	IMPLICIT	25	Auto	8790	0.001	0h01

  • An averaging in space, volume and equivalent plastic strain weighted for triaxiality, 𝜂 2 , and volume weigthed for equivalent plastic strain 𝑝 𝑎𝑣𝜂 𝑎𝑣 4 = 𝑝 𝑎𝑣 𝑖 -𝑝 𝑎𝑣 𝑖-1 𝑝 𝑎𝑣𝑛𝑠 (𝜂 𝑎𝑣 2 ,𝑖-1 + 𝜂 𝑎𝑣 2 ,𝑖 ) 2 . (C.5)• First an averaging in time for triaxiality and then weighted in space with final volume and final equivalent plastic strain weighted𝜂 𝑒,𝑎𝑣 5 (𝑖 𝑒 ) = 𝑝 𝑒 𝑖 -𝑝 𝑒 𝑖-1 𝑝 𝑒𝑚𝑎𝑥 (𝜂𝑒,𝑖-1 + 𝜂 𝑒,𝑖 ) 2 𝜂 𝑎𝑣 5 = ∑ 𝑛𝑒 𝑘 𝜂 𝑒,𝑎𝑣 5 ,𝑘 𝑉 𝑒,𝑓,𝑘 𝑝 𝑒,𝑚𝑎𝑥,𝑘 ∑ 𝑛𝑒 𝑘 𝑉 𝑓,𝑘 𝑝 𝑒,𝑚𝑎𝑥,𝑘

	𝑛𝑠
	∑
	𝑖 𝑝𝑙.

𝑛𝑠 ∑ 𝑖 𝑝𝑙. (C.6) Results for the tensile samples in ch. 3.1 are

Table C . 1 :

 C1 Average triaxialities in sample geometries calculated in different ways.Type𝜂 𝑎𝑣 1 𝜂 𝑎𝑣 2 𝜂 𝑎𝑣 3 𝜂 𝑎𝑣 4 𝜂 𝑎𝑣 5

	APPENDIX D
	Detailed volume based analysis of triaxiality
	7 mm TA6V	0.38 0.47 0.46 0.44 0.45
	0.7 mm TA6V 0.37 0.4 0.42 0.42 0.41
	7 mm U500	0.4 0.47 0.49 0.46 0.47
	0.7 mm U500	0.38 0.41 0.44 0.41 0.41

  Hill criterion are given by (𝑝) is the yield stress, defined by the reference yield curve. 𝜏 0 = 𝜎 0 / √ 3 and 𝜎 𝑦 𝑖𝑗 , 𝜏 𝑦 𝑖𝑗 define the yield stresses of the stress component.

	𝐹 = 𝐺 = 𝐻 = 𝐿 =	𝜎 2 0 2 𝜎 2 0 2 𝜎 2 0 2 3𝜏 2 [ 1 (𝜎 𝑦 22 ) 2 + [ 1 (𝜎 𝑦 33 ) 2 + [ 1 (𝜎 𝑦 11 ) 2 + 0 2 (𝜏 𝑦 23 ) 2 , 𝑀 = (𝜎 𝑦 1 33 ) 2 -1 (𝜎 𝑦 11 ) 2 -1 (𝜎 𝑦 22 ) 2 -3𝜏 2 0 2 (𝜏 𝑦 13 ) 2 , 𝑁 = ] 1 (𝜎 𝑦 , 11 ) 2 ] 1 (𝜎 𝑦 , 22 ) 2 (𝜎 𝑦 33 ) 2 , 1 ]	0 2 (𝜏 𝑦 12 ) 2 , 3𝜏 2	(H.1)
	where 𝜎 0

Table H . 4 :

 H4 Coefficients for Hill criterion and extended Voce hardening identification Sample 𝜎 𝑦 (MPa) 𝑅 ∞ (MPa) Calculate stress 𝜎 𝑛+1 from 𝜎 𝑛 ,𝛥𝜀 𝑛+1 with the constitutive model 5. Calculate internal force and stiffness 𝐹 𝑖𝑛𝑡,𝑛+1 = 𝐴𝜎 𝑛+1 . If 𝑡 𝑚𝑎𝑥 is not reached, go to 1

	6. Find new displacement
	𝑢 𝑖+1 𝑖𝑛𝑐,𝑛+1 = -𝐾 𝑖+1 𝑛+1	-1 𝑅 𝑖+1 𝑛+1 ,	𝑢 𝑖+1 𝑛+1 = 𝑢 𝑖 𝑛+1 + 𝑢 𝑖+1 𝑖𝑛𝑐,𝑛+1 .	(K.5)
	7. If residual is too large, go to 2
	8. If 𝑡 𝑚𝑎𝑥 is not reached, go to 1
	Explicit algorithm
	The principal steps of the explicit calculation are
	1. Increase time and loading
	𝑡 𝑛+1 = 𝑡 𝑛 + 𝛥𝑡,				θ𝑛+1 = θ.𝑡 𝑛+1 .	(K.6)
	2. Calculate new displacement
	𝑢 𝑛+1 = 𝑢 𝑛 + 𝛥𝑡𝑣 𝑛 +	1 2	𝛥𝑡 2 𝑎 𝑛	𝛾	𝐶 (MPa)	(K.7)
	TA6V (L)		1163			163	23	350
	U500 (L) 3. Increase external force and strain 641	177	44	2355
	Sample Ta6V (L) 𝐹 𝑒𝑥𝑡,𝑛+1 = 𝑚(𝑟 0 + 𝑢 𝑛+1 ) θ2 𝑅 11 1.00 𝑛+1 , 𝑅 22 0.87	𝑅 33 0.91 𝛥𝜀 𝑖 𝑛+1 =	𝑅 12 0.82 𝑢 𝑖 𝑛+1 -𝑢 𝑛 𝑅 13 𝑅 23 0.50 0.49 𝑟 0 + 𝑢 𝑛	(K.8)
	U500 (L) 4. (K.9) 1.00 1.13 1.02 1.25 0.13 0.96
	6. Calculate new acceleration
	𝐹 𝑘𝑖𝑛,𝑛+1 = 𝐹 𝑒𝑥𝑡,𝑛+1 -𝐹 𝑖𝑛𝑡,𝑛+1 ,	𝑎 𝑛+1 = 𝑚 -1 𝐹 𝑘𝑖𝑛,𝑛+1	(K.10)
	7. Calculate new velocity
	𝑣 𝑛+1 = 𝑣 𝑛 +	1 2	𝛥𝑡(𝑎 𝑛 + 𝑎 𝑛+1 )	(K.11)
	8.

State of the art

Identification of plasticity at high strain levels and for heterogeneous samples using I-DIC

Perspectives for the simulation of material failure

𝑉 el is the volume of the element "el", 𝐴 el,max = max𝑖{𝐴 𝑖,facet } is the biggest facet area of the element "el".

of the bounded rate model in connection with low loading rates and large time steps. Therefore we limit the following to a discussion of the decrease in simulation time, which is to be expected in the case of the use of a progressive damage model. To accomplish this the classical damage model, without the regularization, was used.

Expected decrease in simulation time

The model of the disk with a notch was simulated under quasi-static loading with the previously mentioned model parameters. To take advantage of the switch from implicit to explicit all simulations were initialized with a dynamic implicit calculation. It was verified that the used time steps lead to a converged solution. At the moment of the switch at 8470 𝑟𝑎𝑑/𝑠, the maximum equivalent plastic strain was 𝑝 = 0.16. This was done to ensure that in the case any temporary errors occur at the instant of the switch, these take place before the onset of damage at 𝑝 0 = 0.25 at 8680 𝑟𝑎𝑑/𝑠. The coarse mesh presented above was used and the element deletion was not necessary, as element distortion only occurred when the damage band had traversed the model.

Table 5.5 provides an overview on the simulation parameters and the corresponding results. Two implicit simulations were carried out. As it is possible to apply the time scaling to the implicit calculations and to not penalize the implicit solution time this was done.

For the current case the critical time is 𝜏 𝑐 = 10 -4 𝑠. Therefore a reference implicit dynamic simulation was performed. In order to be compatible with the critical time used in regularized computations, the time step used was decreased to 𝛥𝑡 = 10 -5 𝑠 at the moment of the switch. The implicit calculation respecting this time step takes about 10 hours of computation. In addition, the results of an implicit simulation, which does not take into account the time step restriction and therefore does not allow mesh-independent simulations with the bounded rate model, are given for comparison. For the explicit simulations, four combinations of mass and time scaling were analyzed.

The values of the limit speed obtained for two different total scaling 𝛼 𝑇 = [2500, 5000] values differ about 6 𝑟𝑎𝑑/𝑠. This corresponds to the delay effect observed for the massspring system. The limit speed for the implicit reference calculation is only θ = 8790 𝑟𝑎𝑑/𝑠. At this point, the maximum damage value is 𝐷 = 0.45. For the explicit calculation with 𝛼 = 2500, this value is reached at θ = 8799 𝑟𝑎𝑑/𝑠, very close to the limit speed of θ = 8801 𝑟𝑎𝑑/𝑠. One can therefore assume that θ = 8790 𝑟𝑎𝑑/𝑠 is the very close to the real limit speed of an unscaled solution. Then the limit speed error is only 0.1% for the scaling of 𝛼 = 2500. This is remarkable considering the decrease in simulation time; all explicit simulations remain below 10 minutes time.

In addition the final value of the error in kinetic energy 𝜖 2 is provided. For the implicit reference simulation 𝜖 2 = 0.03, while for the scaled simulations values are in the range between 1.6-3.3%. This is a difference by two orders of magnitude. This factor corresponds to the order of magnitude in the difference of the total scaling factor between implicit and explicit calculations i.e. 𝛼 = 25 and 𝛼 𝑇 = [2500, 5000]. The relation in equation (5.24) is confirmed by this result.

Examining the errors for the same total scaling, one finds that it is better to use mass scaling than time scaling as for 𝛼 𝑀 = 500, 𝛼 𝑇 = 5 𝜖 2 = 1.6 % compared to

The angle 𝜑 defines the size of the yield surface in relation to the hydrostatic pressure. For a dependence on the triaxiality the states circle becomes smaller, if the principal stresses are larger. The Lode angle 𝜃 describes a variation of the criterion in the deviatoric plane. From the definition of θ it can be recognized that θ = -1 corresponds to an axial symmetry loading in compression, θ = 0 to a plane strain state or a generalized shear one and θ = 1 represents an axial symmetry loading in tension [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF].

APPENDIX B

Relationship between maximum principal stress and triaxiality

Triaxiality is defined as

The maximum principal stress can then be written as

APPENDIX C

Calculation of average triaxiality values

The calculation of the average quantities from FE data is not a well-defined procedure.

First, the ROI needs to be defined. In the following this always is the entire model. Second, the averaging of the values depends on the way the values are calculated from the FE data. The following definitions are used: element volume 𝑉 𝑒 , element triaxiality 𝜂 𝑒 , element equivalent plastic strain 𝑝 𝑒 and the step number 𝑖 𝑠 . Where 𝑛 𝑠 = 𝑖 𝑝𝑙. ..𝑖 𝑒𝑛𝑑 only contains the steps where equivalent plastic strain occurs.

• Volume weighted

𝜂 𝑎𝑣 1,𝑖 .

(C.1)

• Volume and equivalent plastic strain weighted in two steps

𝜂 𝑎𝑣 2,𝑖 .

(C.2)

• Volume and equivalent plastic strain weighted in one step

∑

APPENDIX E

Variation of the identification with regard to maximum deformation

The following study concerns the test of a thin longitudinal sample, whose results are presented in section 3.3.2. It is studied how the identification behaves if less than these 47 images are used for the identification, i.e. if the test data is not analyzed up to fracture. Image sequence of 17, 27, 37 of a total of 47 images at fracture were analyzed. Each series is identified with the maximum macroscopic strain: 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.01, 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.03, 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.06 and for the total series 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.07. The optimal material model for each series was obtained, cf. figure E.1(a). The identified model was then applied to the entire test data. For the first series, the material behavior is mainly elastic. The predict forces are only correct for elasticity. The identification 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.03, predicts a correct yield stress but wrong hardening. The forces are too low. The hardening curve for 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.06 is marginally above the curve for 𝜀 𝑚𝑎𝑐𝑟𝑜 = 0.07, but both curves led to almost identical forces.

APPENDIX H

Identified parameters for plasticity

Von Mises criterion Verification of all data points from all steps (red) and last step (blue) for two of the criteria

APPENDIX J

Implementation of the bounded rate model

The bounded rate model was implemented using the theta scheme. The constitutive equation is

The solution for a time step 𝑛 + 1 is

The solution to this non-linear equation is obtained with the newton method. For 𝜃 = 0 the scheme becomes explicit. For 𝜃 = 0.5 it is implicit/explicit. For a linear hardening law the bounded rate model reduces to

3)

The analytical solution to this equation is possible. Supposing 𝑝(𝑡 = 𝑡 𝑛 ) = 𝑝 𝑛 and 𝑝(𝑡 = 0) = 0

APPENDIX K

Calculation schemes for a mass spring system

To model an elasto-plastic behavior the stiffness of the spring is calculated by assuming a beam element of cross section 𝐴 0 and length 𝑟 0 . The mass of the beam 𝑚 = 𝐴 0 .𝑟 0 𝜌 is equally distributed to each end. The non-linearity requires the iteration 𝑖 of the solution for the instant 𝑢 𝑛+1 . The effects of the large displacements intervene in the calculation of the strain increment 𝛥𝜀, the calculation of the stiffness 𝐾, and the calculation of the external force. The effect on the external force is associated with spin softening. The calculation of the stress 𝜎 and the tangent modulus follows from the constitutive model.

Implicit algorithm

The main steps of the implicit calculation are