2019troy0010 Année

Docteur

Kai Zhang

J.-M Bergheau

Professeur Des

Universites Président

P.-O Bouchard

Professeur -Hdr Rapporteur

PROFESSOR DOKTOR Mme S Reese Universitat

Ingegnieur -Hdr Rapporteure

M C Labergère

Universites Examinateur

M A Roos Docteur

Ingenieur -Hdr Examinateur

M H Badreddine

Maitre De Conferences

M K Saanouni

Modeling of Initial and Induced Anisotropies in Hot Sheet Metal Forming in Presence of Ductile Damage Champ disciplinaire : Sciences pour l'Ingénieur

Keywords:

Third invariant of the second-rank tensor det( )

. The total energy is assumed to be the sum of the elastic (reversible) energy We, the energy dissipated in the kinematic hardening Wα and the energy stored in the isotropic hardening Wr. This total energy is assumed to be the same for both configurations. It enables the systematic definition of pairs of effective state variables that can be indifferently used both in stress and in strain spaces with the Legendre-Fenchel transformation. It not only achieves a strong coupling between damage and all other phenomena, but also can keep the interaction between the various phenomena.

For the isothermal problems the observable variables are the total strain and stress tensors   ,   , and the internal variables represent different dissipative phenomena as the inelastic flow ( , ) e   , the isotropic hardening   , r R ,the kinematic hardening   , X  and the isotropic ductile damage   , d Y . This suggests the existence of a fictitious undamaged configuration (where d = Y = 0), in which the state space is defined by the couples of effective variables   ( , ) e   ,   ( , ) X  ,  ( , ) r R 

. The assumption of total energy equivalence postulates that the two configurations have the same total energy so that the effective and real variables are related by the following equations:

( , ) 1/ 2( : ) 1/ 2( : ) ( ) ( ) , ( )

IV

Abstract

This work is dedicated to the development of an advanced behavior constitutive model coupling various physical phenomena, for the precise prediction of the behavior and ductile damage of thin sheets during their forming at high temperature. Interest is mainly focused on hexagonal-close-packed (HCP) metals such as magnesium and titanium alloys. The proposed constitutive equations take into account: (i) initial and induced anisotropies, (ii) tension-compression asymmetry, (3i) isotropic and kinematic mixed nonlinear hardening, (4i) yield surfaces distortion induced by the loading path, (5i) loading rate dependence within a unified viscoplastic formulation, and (6i) strong coupling between the various phenomena including isotropic ductile damage. In particular, the coupling with ductile damage takes into account the microcracks closure with a novel approach based on a dual dependence on the stress triaxiality and the Lode angle. The formulation of the model is carried out within the framework of the thermodynamics of irreversible processes with state variables and considering a unified theory of non-associative viscoplasticity in finite transformations. Associated numerical aspects are developed as part of the finite element (FE) method and implemented in the ABAQUS/Explicit FE calculation code via the VUMAT user routines. A parametric study is systematically performed to show the predictive capabilities of the proposed model. The procedure for identifying model parameters is then discussed. This procedure is based on minimizing the difference between the experimental measurements and the numerical simulation results by considering the responses of the materials for different tests conducted until the fracture occurrence. Finally, FEM numerical simulations of various high temperature thin sheet forming processes are performed and their comparisons to the experimental results are presented and analyzed.

V

Résumé

Ce travail est dédié au développement d'un modèle de comportement avancé couplant divers phénomènes physiques pour la prédiction fine du comportement et de l'endommagement en simulation des procédés de mise en forme de tôles minces à hautes températures. L'intérêt est principalement porté aux métaux à structure hexagonale compacte (HC) comme les alliages de Magnésium et de Titane. Les équations constitutives proposées prennent en compte : (i) les anisotropies initiales et induites, (ii) l'asymétrie de traction-compression, (3i) un écrouissage mixte non linéaire isotrope et cinématique, (4i) la distorsion de la surface de charge induite par le trajet de chargement, (5i) la dépendance à la vitesse de chargement dans le cadre d'une formulation viscoplastique unifiée, et (6i) le couplage fort entre les différents phénomènes incluant l'endommagement ductile isotrope. Particulièrement, le couplage avec l'endommagement ductile tient compte de la fermeture des microfissures avec une nouvelle approche originale basée sur une double dépendance à la triaxialité des contraintes et à l'angle de Lode. La formulation du modèle est réalisée dans le cadre de la thermodynamique des processus irréversibles considérant une théorie unifiée de viscoplasticité non-associée en transformations finies. Les aspects numériques associés, sont développés dans le cadre de la méthode des éléments finis (EF) et implémentés dans le code de calcul par EF ABAQUS/Explicit via des routines utilisateurs VUMAT. Une étude paramétrique est systématiquement réalisée pour montrer les capacités prédictives de la modélisation proposée. La procédure d'identification des paramètres du modèle est ensuite discutée. Cette procédure est basée sur la minimisation de la différence entre les mesures expérimentales et les résultats des simulations numériques en considérant les réponses des matériaux dans différents essais menés jusqu'à rupture. Finalement, des simulations numériques par EF de quelques procédés de mise en forme de tôles minces à hautes températures sont réalisées et leurs comparaisons aux résultats expérimentaux sont présentées et analysées.

VI

Notations

First-rank tensor or vector

Introduction

The increasing demand on lightweight engineering materials in automotive and aerospace industries has increased dramatically due to the urgent need for improving fuel efficiency and reducing CO2 emissions.

With a low density and a high strength to weight ratio, high strength (HS) metallic materials (e.g., Magnesium (Mg) alloys, Titanium alloys…) become ideal candidates. However, their formability is highly decreased when performed at room temperature. The simplest way to increase the formability properties of these HS metallic materials is to perform their forming processes under adequately high temperature usually called hot forming. On the other hand, it is well known that this kind of materials exhibit complex microstructure leading to highly complex mechanical behavior at room temperature. Note that, performing hot forming processes is not a new task by itself since it is widely used for bulk metal forming (forging, stamping, extrusion…). However, for sheet metal forming, mainly performed at room temperature, it is still a novel task. This requires the use of advanced and highly predictive constitutive models together with associated numerical methods.

For the case of Mg alloys, the industrial applications are limited by different kinds of pronounced initial and induced anisotropies [START_REF] Agnew | Plastic anisotropy and the role of non-basal slip in magnesium alloy AZ31B[END_REF][START_REF] Khan | Mechanical response and texture evolution of AZ31 alloy at large strains for different strain rates and temperatures[END_REF][START_REF] Shi | On the macroscopic description of yield surface evolution by means of distortional hardening models: application to magnesium[END_REF], the strong Strength Differential (SD) effect (tensioncompression asymmetry) [START_REF] Kelley | Plane-strain compression of magnesium and magnesium alloy crystals[END_REF][START_REF] Khan | Mechanical response and texture evolution of AZ31 alloy at large strains for different strain rates and temperatures[END_REF][START_REF] Steglich | Modelling direction-dependent hardening in magnesium sheet forming simulations[END_REF], the hardening asymmetry [Khan11, [START_REF] Kabirian | Plastic deformation behavior of a thermo-mechanically processed AZ31 magnesium alloy under a wide range of temperature and strain rate[END_REF][START_REF] Lee | Distortional hardening concept for modeling anisotropic/asymmetric plastic behavior of AZ31B magnesium alloy sheets[END_REF] and the poor formability at room temperature [Chen03, [START_REF] Kim | A temperature-dependent elasto-plastic constitutive model for magnesium alloy AZ31 sheets[END_REF][START_REF] Rodriguez | Effect of strain rate and temperature on fracture of magnesium alloy AZ31B[END_REF]. These properties are related to their Hexagonal-Close-Packed (HCP) crystalline microstructure. Consequently, the modeling of these complex material behaviors is still today a highly challenging task.

Efforts have been made to describe the material behavior of these microstructurally complex materials at various scales. Crystal plasticity theory [Jain07; Argon08; Proust09, Yoshida16] provides a naturally convenient framework to deduce the macroscopic (i.e. valid for an aggregate of single crystals) material behavior from the individual behavior of each crystal. This aims to relate the single crystal constitutive behavior with that of the overall aggregate with full field approaches or mean field approaches [START_REF] Lebensohn | Self-consistent modelling of the mechanical behavior of viscoplastic polycrystals incorporating intragranular field fluctuations[END_REF][START_REF] Perdahcıoğlu | Constitutive modeling of two phase materials using the mean field method for homogenization[END_REF] as the self-consistent approach [Lebensohn93; [START_REF] Saanouni | Micromechanical modeling of low cycle fatigue under complex loadings-PartI. Theoretical formulation[END_REF][START_REF] Boudifa | A micromechanical model for inelastic ductile damage prediction in polycrystalline metals for metal forming[END_REF] which are widely used to model the interaction of a grain with the surroundings. The initial anisotropy and tension-compression asymmetry of magnesium alloys were investigated by many researchers [Agnew01;

Agnew05; Lou07; Guo15; Zhou16] at the microscale. It is noted that the microscopic models are useful for capturing the macroscopic mechanical behavior of these highly heterogeneous metals. However, to fulfill the aim of modeling large-scale forming processes with acceptable computational efficiency, a macroscopic phenomenological approach is more suitable. On the other hand, the simple associative plasticity theory with a single yield function is not efficient to describe the various strong initial and induced anisotropies of the advanced metallic materials [Stoughton02]. Accordingly, the use of non-associative flow rule to describe the anisotropic plastic flow of metals is more convenient. In addition, the accurate prediction of the inelastic flow localization giving rise to ductile damage initiation and macroscopic cracks propagation under various complex loading paths require advanced models accounting for the strong coupling between the behavior and the damage occurrence [START_REF] Cao | A comparative study of three ductile damage approaches for fracture prediction in cold forming processes[END_REF]. The framework of continuum damage mechanics (CDM) has been developed and applied to describe the ductile damage with including complex material behaviors in the last decades [Lemaitre85, Saanouni03, Bonora05, Besson09, [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF][START_REF] Murakami | Continuum damage mechanics: A continuum mechanics approach to the analysis of damage and fracture[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. Due to its convenience in coupling damage with different material phenomena, the CDM approach is adopted in this work in order to describe the damage occurrence in metallic materials formed by large anisotropic inelastic strains.

With this motivation, the aim of this thesis is dedicated to develop a thermodynamically consistent phenomenological model to well describe the complex behaviors which can be applied to hot sheet metal forming processes with higher accuracy. The modeling of this problem involves many strong couplings with large inelastic strains, kinematic and isotropic hardenings, thermal effects, and the ductile damage. Based on the theoretical and numerical framework developed by our research team (LASMIS), an anisotropic thermo-elasto-viscoplastic model with non-linear mixed hardening (isotropic, kinematic and distortional hardening) strongly coupled with isotropic ductile damage is proposed in this thesis. Microcracks closure effect and Lode angle dependence are carefully taken into account. Asymmetry in yield stresses is captured by introducing the J3 invariant into Hill yield function, and the description of hardening asymmetry is realized by Lode angle dependence function of hardening parameters.

The formulations of the model are performed in the framework of thermodynamics of irreversible processes using generalized non-associative theory under finite transformations. The associated numerical aspects are developed in the framework of finite element method and implemented in ABAQUS/Explicit FE code via the users' developed subroutines VUMAT.

The calibration of the anisotropic parameters is conducted by minimizing the objective function between the numerical results and experimental observations [START_REF] Yoon | Asymmetric yield function based on the stress invariants for pressure sensitive metals[END_REF].To identify the behavior and damage parameters, an inverse methodology is applied which is based on minimum error value between simulation and experimental responses. This inverse procedure was realized throughout Python script which combines the ABAQUS/Explicit FE software with the MATLAB-based minimization code [START_REF] Yue | Identification of fully coupled anisotropic plasticity and damage constitutive equations using a hybrid experimental-numerical methodology with various triaxialities[END_REF][START_REF] Souto | Material parameter identification within an integrated methodology considering anisotropy, hardening and rupture[END_REF].

An exhaustive parametric study is conducted to check the ability of the proposed model in capturing the complex phenomena characterizing the behavior and ductile damage of the materials. Finally, applications are investigated to some sheet metal forming processes under high temperature with comparison to some available experimental results.

The thesis is organized as following:

Chapter I starts by a review of sheet metal forming processes, then presents the complex phenomena exhibited by the innovative lightweight materials in both a typical hot and cold sheet metal forming processes.

Chapter II gives the fundamentals of finite strain viscoplasticity theory, as the kinematics of elastoviscoplastic decomposition, strain and stress measures as well as the main conservation laws. A review of the existing ductile damage models is presented as well as the description on physical aspects of ductile damage. A literature survey about the development of time-independent plasticity and viscoplasticity theories to capture complex material behaviors mentioned in Chapter 1 is given based on some published works. The proposed anisotropic thermo-elasto-visco-plastic constitutive equations fully coupled with isotropic ductile damage with embracing many phenomena observed in metal viscoplasticity (initial and induced anisotropies, SD effect, hardening asymmetry, etc.) are given in Section 2.6. A new formulation of the microcracks closure damage effect is developed and discussed. Finally, a short discussion of the contact and friction in metal forming are given.

Chapter III is dedicated to the numerical aspects for solving the equilibrium problems with fully coupled thermo-mechanical equations under large inelastic strains with damage effect. First, the strong and weak forms of the initial and boundary value problem (IBVP) are defined. Then the time and space discretization of the IBVP leading to obtain of highly nonlinear and strongly coupled algebraic system is given. The dynamic explicit resolution scheme is used to solve the IBVP while an iterative local integration scheme is used to perform the local integration of fully coupled ordinary differential equations (ODEs) associated with thermo-elasto-viscoplastic constitutive equations with damage.

In Chapter IV, a systematic parametric study is conducted with RVE. The effects of anisotropy, tensioncompression asymmetry, temperature, strain rate, hardening asymmetry parameters are carefully investigated. The combined effect of distortion of the yield surface and tension-compression asymmetry is studied. The triaxiality and Lode angle effect on the ductile fracture locus are carefully discussed.

In chapter V, the parameters determinations and applications of fully coupled CDM model are given. Based on the available experimental results, the material parameters are calibrated. Finally, the applications of the proposed damage model to three point bending test (TPB), circular cup deep drawing test (CCD) and crossshaped cup deep drawing test (CSD) are performed and their results are discussed.

Finally, the main conclusions and some perspectives of the present work are presented.

Introduction

Sheet metal forming is widely used for automotive production among many other industrial sectors. In these processes of sheet metal forming, the final shape of a part is made from a flat metal sheet. First, sheet plate is cut into pieces by different operations like shearing, slitting, cutting, sawing or produce from coil; then the desired shape can be achieved by punching, stamping, embossing, bending, stretch forming, deep drawing and a variety of other processes. During these forming processes, the changes in the shape of the sheet are mainly obtained by large plastic deformation. In this chapter, the sheet metal forming processes are presented, a general introduction to sheet materials as well as the major physical phenomena involved during their forming by large inelastic strains including material behaviors and different failure mode are given.

General presentation of sheet metal forming

Metal forming processes are suitable for quantity production of components with high production rates and a minimum waste of material. Among the metal forming processes, sheet metal forming processes are commonly used in many industries, such as automotive industry (car doors, hoods, fenders, etc.), household products (cooker hoods, refrigerators, sinks, etc.), food industry (cookware, canned goods, etc.), aerospace, ships, etc. In recent years, the demand of lightweight engineering materials in automotive and aerospace industries has increased significantly due to the urgent need for reducing energy consumption as well as CO2 emissions. Accordingly, nowadays lightweight construction design concept has become a popular term in new industrial development. However, the formability of certain lightweight metals (i.e. Mg alloys) at room temperature are not satisfied for large-scale industrial application. To avoid this limitation, the forming processes are performed under high temperature. Since the ductility of the material is expected to increase with the temperature increase, the formability of these metals is quitely improved at elevated temperature.

Meanwhile, forming forces and applied pressure are reduced as a result of the material ''softening'' induced by the temperature increase. These processes are commonly called warm or hot forming. From literature, there are two kinds of definition for material forming at elevated temperature: (1) when the heating temperature is above the temperature of recrystallization, the forming process is hot forming, warm forming means the heating temperature is between the room and recrystallization temperatures.

(2) The forming process is classified to cold forming which has no pre-heating, otherwise, it is called hot forming if there is heating process. In this thesis, we will focus on the second way to define the hot forming. Typical hot Chapter 1: Major physical phenomena exhibited in hot sheet metal forming 6 forming process is shown in Figure 1.1, the metal sheet is processed in heating, transporting to forming process, followed by quenching and ageing, finally, the desired form is obtained.

Materials and major physical phenomena 1.3.1 Target materials of hot sheet metal forming

The innovative lightweight construction materials, such as magnesium alloys, titanium alloys, aluminum alloys with excellent material properties are very promising materials in recent years, since these materials are widely demanded in many industrial fields, as shown in Figure 1.3. Aluminum alloys have been already

widely used in automotive manufacturing (frames, car body panels and wheels, etc), ship building, as well as aircraft and other aerospace structures. Titanium alloys are often used in aerospace applications such as for turbine blades engines, airframe applications and space applications, due to their superior strength to weight ratio. They are also used in offshore petroleum industry, such as drilling risers, because of their high flexibility (low ratio of modulus over strength), excellent corrosion and fatigue resistance.

However, their application in automotive industry is impeded by their high price. Magnesium alloys, with a high ratio of strength versus density, are excellent candidate for use in automotive applications aiming to enable lightweight design and improve fuel economy, which also makes them a good alternative to traditional steel and aluminum alloys. For instance, magnesium alloys can save 20%-25% weight over aluminum alloys in the car of full-size "Savana" and "Express" vans of General motors [START_REF] Elektron | College Street, P.O. Box[END_REF][START_REF] Jia | Experimental study on the mechanical properties of AZ31B-H24 magnesium alloy sheets under various loading conditions[END_REF].

In addition to these materials, the high-strength and ultra-high strength steel materials, with a high specific strength in comparison with conventional steel materials, also demonstrate their capability in reducing the weight of car component. However, these innovative lightweight construction materials and the highstrength/ultra-high strength steel materials have limitations in terms of forming capability and some of them are classified as hard-to-deform materials at room temperature. As a result, they become the target materials of hot sheet metal forming.

Material behaviors

Temperature and strain rate dependent material behavior

A fatal drawback of the innovative lightweight materials is that they often exhibit very poor formability compared with traditional steels. For example, the poor formability of Magnesium alloys at room temperature is due to the limited number of active slip systems in their hexagonal close-packed (HCP)

crystal structure [START_REF] Agnew | Plastic anisotropy and the role of non-basal slip in magnesium alloy AZ31B[END_REF]Wu15,[START_REF] Trang | Designing a magnesium alloy with high strength and high formability[END_REF]. It is important to recall that the plastic strain of HCP materials is divided into slip and twinning mode. The activation of these modes is highly dependent on both critical stress and loading directions. Figure 1.4 shows the loading directions to activate preferentially one of these Even though twinning and non-basal slip are activated to some extent at room temperature, the improvement with respect to formability is very small. In order to form these kind of materials with complex component geometries, additional glide planes should be activated. The use of hot forming can be applied to overcome this shortcoming of Mg alloys, since pyramidal slip system can be thermally activated at elevated temperature, which improve their ductility and formability obviously [START_REF] Mekonen | Experimental and numerical investigation of Mg alloy sheet formability[END_REF]. Figure 1.5a shows the temperature dependency of inelastic flow of Mg AZ31B, at a given strain rate, significant thermal softening caused by the temperature increase can be observed, unlike the room temperature behavior, ductility and formability are greatly improved with elevated temperatures [START_REF] Khan | Mechanical response and texture evolution of AZ31 alloy at large strains for different strain rates and temperatures[END_REF][START_REF] Tari | Elevated temperature constitutive behavior and simulation of warm forming of AZ31B[END_REF][START_REF] Rodriguez | Effect of strain rate and temperature on fracture of magnesium alloy AZ31B[END_REF]. The forming force and the spring back effect are decreased with the temperature increases. The strain rate has very important influence on the flow curves, as illustrated in Fig. 1.5b, the flow stress increased obviously with the strain rate increase. With the temperature increase, the strain rate sensitivity becomes more evident, as shown in Figure 1.5c, about 85% decrease in flow stress and 100% increase in elongation to failure at temperature of 300°C over the same range of strain rate can be found. For aluminum and titanium alloys, the same temperature dependent material behavior can also be observed [START_REF] Abedrabbo | Forming of aluminum alloys at elevated temperatures-Part 1: Material characterization[END_REF][START_REF] Tabei | Constitutive modeling of Ti-6Al-4V at a wide range of temperatures and strain rates[END_REF]. strain rate at temperature of 100°C; (c) Different strain rate at temperature of 300°C [START_REF] Rodriguez | Effect of strain rate and temperature on fracture of magnesium alloy AZ31B[END_REF].

Chapter 1: Major physical phenomena exhibited in hot sheet metal forming 9

Initial plastic anisotropy and tension-compression asymmetry

Since the final sheet metals are obtained by successive rolling operations, a clear initial anisotropy is conferred to these sheets [START_REF] Banabic | Anisotropy of sheet metal[END_REF]. This initial anisotropy has a strong effect on the formability of sheet metals [START_REF] Vladimirov | On the influence of kinematic hardening on plastic anisotropy in the context of finite strain plasticity[END_REF]. As shown in Figure 1.6, the initial yield stress and R-value of titanium alloys Ti-6Al-4V in tension varying significantly with the material orientations [START_REF] Gilles | Experimental characterization and elasto-plastic modeling of the quasistatic mechanical response of TA-6 V at room temperature[END_REF]. Another pronounced anisotropic behavior has been observed in some Mg alloys (e.g. AZ31), as shown in Figure 1.7, where the yield stress in tension is much higher than that in compression. This behavior is often called tension-compression asymmetry or Strength Differential (SD) effect [Kelly68, [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF][START_REF] Cazacu | Orthotropic yield criterion for hexagonal closed packed metals[END_REF]. This special feature of magnesium alloys results from the activation of deformation twinning under compression stress state. For titanium Ti-6Al-4V, the compressive yield stress is higher than the tensile yield stress [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], but the asymmetric behavior is not pronounced. The tension-compression asymmetry in yield stress is also observed in other materials (i.e. DP980 steel), as reported in [START_REF] Holmen | Strength differential effect in age hardened aluminum alloys[END_REF][START_REF] Maeda | Measurement of the strength differential effect of DP980 steel sheet and experimental validation using pure bending test[END_REF].

Induced anisotropy

In practice of sheet metal forming, the sheet metals are formed with large plastic strains under complex loading paths, the initial anisotropies evolve due to the textural changes inside the material. This affects the hardening evolution and consequently causes the yield surface distortion. As shown in Figure 1.8, the evolution of the yield surface of aluminum AL1100 is highly dependent on the loading path [START_REF] Khan | Evolution of subsequent yield surfaces and elastic constants with finite plastic deformation. Part II: A very high work hardening aluminum alloy (annealed 1100 Al)[END_REF]. The initial yield surface evolves to subsequent yield surfaces through translation, distortion, extension or shrinkage, represented by a sharp nose in the loading direction and flat in the opposite direction, depending on the applied loading path. Other cross-effects were correlated to the different hardening properties. The induced anisotropic behavior is also reported by [START_REF] Shi | Loading path dependent distortional hardening of Mg alloys: Experimental investigation and constitutive modeling[END_REF] for magnesium alloy AZ31 and by [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] for titanium alloy Ti-6Al-4V.

Figure 1.8: Distortion of the yield surface depending on loading path [START_REF] Khan | Evolution of subsequent yield surfaces and elastic constants with finite plastic deformation. Part II: A very high work hardening aluminum alloy (annealed 1100 Al)[END_REF].

Hardening asymmetry

In addition to the anisotropy and asymmetry of yielding, Mg alloys also exhibit unusual hardening evolution [START_REF] Agnew | Plastic anisotropy and the role of non-basal slip in magnesium alloy AZ31B[END_REF][START_REF] Guo | On twinning and anisotropy in rolled Mg alloy AZ31 under uniaxial compression[END_REF] compared to other materials, the stress-strain flow evolution has a sigmoidal shape [START_REF] Khan | Mechanical response and texture evolution of AZ31 alloy at large strains for different strain rates and temperatures[END_REF][START_REF] Kabirian | Plastic deformation behavior of a thermo-mechanically processed AZ31 magnesium alloy under a wide range of temperature and strain rate[END_REF], which is referred as hardening asymmetry [START_REF] Lee | Distortional hardening concept for modeling anisotropic/asymmetric plastic behavior of AZ31B magnesium alloy sheets[END_REF].

Figure 1.9: Typical stress-strain evolution curve of Mg alloy AZ31B under uniaxial tension and compression at room temperature [Zhang 19].

The evolution of compressive stress-strain response of Mg alloy AZ31, as shown in Figure 1.9, follows a sigmoidal evolution with a lower flow stress plateau and upper flow stress plateau with an inflection point.

At lower strain stage, the plastic flow is governed by the twinning mechanism, where a low stress plateau was observed, with the plastic strain increase, the slip-dominated deformation becomes active, with a high hardening rate, and finally the upper flow stress plateau appears.

Damage and fracture 1.4.1 Ductile damage

It is vital to understand the mechanisms that cause the inelastic strain localization which eventually leads to ductile damage occurrence in a sheet metal component. Damage is regarded as micro-defects density, which accumulates to a certain level before the initiation of macroscopic cracks. From the physical point of view, ductile damage is described as nucleation, growth and coalescence of microcavities or microcracks induced by large inelastic strains in the vicinity of inclusions in material. In metal forming by large irreversible inelastic strains, the damage mechanism is often ductile fracture occurring in areas where large inelastic strains strongly localize, rather than brittle fracture occurring without appreciable plastic strain. In sheet metal applications three failure phenomena are usually observed [START_REF] Björklund | Ductile Failure in High Strength Steel Sheets[END_REF]: (i) Ductile tensile fracture,

Creep damage

In addition to ductile damage, there are other damage types like creep damage under high temperature, fatigue damage under cyclic loading, etc. Normally, for metallic materials, the yield stress decreases at elevated temperature, meanwhile, the material behavior is more sensitive to strain rate. Creep damage takes place at high temperature, from mid to elevated temperatures that are above about one third of the absolute melting temperature for metals. The creep damage mechanisms are thermally activated phenomena and the microdefects initiate and evolve along the grain boundaries often, starting from initial defects as triple-joint. be taken into account for life assessment during design procedure of the components which are used at high temperature and subjected to static loading paths (i.e. loading paths with relatively low velocity). However, in this work dealing with sheet metal forming where the loading velocity is high enough, the creep damage is neglected and only ductile damage is considered, since the targeted sheet metal forming applications for the lightweight engineering materials are performed at elevated temperature subject to rapidly evolving loading paths which cannot result in creep damage.

Example of defects of target material in metal forming

For the sake of simplicity, the forming tools are assumed perfectly made without any defect and only the defects of the deforming sheet or blank in the forming process are considered. For example, in the deep drawing process, one of the most widely used forming process, three currently observed defects that occur in deep drawing operations for sheet metals are:

Chapter 1: Major physical phenomena exhibited in hot sheet metal forming 13

(1) Earing

The light-weight material with pronounced initial anisotropy due to its crystallographic texture induced by rolling process exhibits earing. The planar anisotropy in the sheet metal drawing leads to the formation of irregularities in the upper edge of a deep drawn cup, which is called earing, as shown in Figure 1.12d.

However, if the plastic flow in the material is perfectly isotropic, earing will not appear.

(

Conclusions

In this chapter, we have presented the major physical phenomena exhibited by the thin sheet materials in both typical hot and cold sheet metal forming processes. In order to predict the damage occurrence of these kind of metallic materials, advanced mechanical models will be developed accounting for the full coupling between all these phenomena including the ductile damage. It should bear in mind that ductile damage is the natural consequence of large hardened inelastic strains localization in narrow (shear) bands. Inside these highly localized zones, inelastic flow, hardening, damage and thermal exchanges are highly active.

Consequently strong interactions between all these phenomena cannot be ignored. It should be noted that Chapter 1: Major physical phenomena exhibited in hot sheet metal forming 14 the reliability of the simulation of the forming processes is highly dependent on the constitutive equations prediction capabilities. The use of fully coupled thermo-mechanical behavior with damage in numerical simulations can be:

(1) In order to avoid the damage initiation during some bulk and sheet metal forming processes and ensure formed parts without defects.

(2) Or contrarily, in order to accelerate the damage initiation and growth in sheet metal cutting, blanking or metal machining by chip formation among others.

In this work, we focus mainly on the prediction of ductile damage which may occur during metal forming in order to avoid its occurrence. Even without final fracture, the ductile damage can locally take place in formed parts leading to strong changes in the material properties. The formulation of an advanced thermoelasto-viscoplastic model fully coupled with ductile damage and accounting for the main aspects characterizing HCP materials will be developed in more detail in next chapter.

Introduction

This chapter is dedicated to the theoretical aspects related to the formulation of the fully coupled constitutive equations in the framework of thermodynamics of irreversible processes with pairs of state variables. In order to present the developed model, some fundamentals of the finite viscoplasticity theory and some basic elements of mechanics and thermodynamics of the materially simple continua are given. These details are taken from literature mainly in [Callen60, Lemaitre85, Lemaitre92, Khan95, Saanouni12], among others.

Next, the framework of Continuum Damage Mechanics (CDM) is presented and used to develop the full coupling with the isotropic ductile damage [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF][START_REF] Lemaitre | A Course on Damage Mechanics[END_REF][START_REF] Murakami | Continuum damage mechanics: A continuum mechanics approach to the analysis of damage and fracture[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. Some developments on anisotropic elasto-visco-plasticity theory are reviewed. Finally, an anisotropic thermoelaso-visco-plastic constitutive equations fully coupled with isotropic ductile damage for metallic materials under large strains are formulated.

Fundamentals of finite viscoplasticity

Homogeneous transformation and deformation gradient

In line with the standard notations in continuum mechanics, a deformable solid body is assumed to occupy at time 0 t  the initial undeformed configuration 0 C , while t C is the current configuration (deformed) at time 0 t  , as shown in Figure 2.1 [START_REF] Sidoroff | The geometrical concept of intermediate configuration and elastic-plastic finite strain[END_REF][START_REF] Lemaitre | A Course on Damage Mechanics[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. Every material point 0 P in the reference configuration 0 C can be identified by its position vector

1 1 2 2 3 3 X X e X e X e        , it can be transformed into t Pof coordinates ( , ) x X t   in the current configuration t C by deformation mapping function ( , ) X t    according to : ( , ) ( , ) x X t X u X t           (2-1)
In which ( , )

u X t  
is the displacement vector of the material point 0 P at time t. The gradient of the homogeneous transformation between 0 C and t C , is the second-rank operator F (or ij F ) which is also called deformation gradient. In the classical local (or Cauchy) continuum, the transformation gradient F allows the complete description of the changes in shape, size and orientation of the continuum.

( ) ( , ) / ( ( , )) / 1 ( , ) / 1 ( ) F Grad x X t X X u X t X u X t X Grad u                           (2-2)
It bears to emphasize that the value of det( ) J F  must verify 0 J    to define correctly the motion of the local continuum.

Polar decomposition of deformation gradient

According to the polar decomposition theorem, the deformation gradient F can be decomposed into a pure rotation (orthogonal tensor) and a pure stretch (symmetric tensor) in the following form:

F R U V R     (2-3)
where R is the rigid body orthogonal rotation tensor, U is the right (Lagrangian) symmetric stretch tensor of second-rank defined with respect to 0 C , V is the left (Eulerian) stretch tensor (second-rank) defined with respect to t C (as shown in Figure 2.2).

Strain and strain rate tensors

In order to define and quantify the change of distance between two material points between original and deformed configuration, we define the right Cauchy-Green stretch tensor C defined in 0 C and left Cauchy-

Green stretch tensor B defined in t C by: 2 2 . . . . . . . . T T T T C F F U R RU U B F F V R R V V            (2-4)
A simple way to obtain the material deformation in the homogeneous transformation is to calculate the difference between the scalar products of the elementary vectors .

.

dx dx dX dX       
, as given below:

  1 . . . 1. . 1 . dx dx dX dX dX C dX dx B dx                       (2-5)
Accordingly, the Green-Lagrangian strain tensor E and the Euler-Almansi strain tensor A are defined respectively as:

2 1 1 ( ) ( ) 2 2 E U I C I     (2-6) 2 1 1 ( ) ( ) 2 2 A I V I B     (2-7)
Numerous strain measures can be obtained by using the Lagrangian and Eulerian stretch tensors. For instance, the overall Lagrangian definitions can be rationalized in the forming form (with m is an integer).

1 0 ln( ) 0 m U I if m m U i f m               (2-8)
The time derivatives of the transformation gradient from Eq. (2-1) gives:

1 dx FdX FF dx Ldx             (2-9)
where the spatial velocity gradient Lwhich can also be decomposed into a symmetric tensor D (total strain rate tensor) and a skew symmetric tensor W (material spin or rotation rate tensor), defined as:

1 1T T L FF RUU R RR D W           ( 2 -1 0 ) 1 1 . . [ . ]. S S S T D L F F R UU R            (2-11) 1 1 . . [ . ]. . A A A T T W L F F R UU R R R             
(2-12)

Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 20

Stress tensors

The Cauchy stress is the most widely used measure of stress to describe the surface forces. As illustrated in

Figure 2.3:

Representation of all quantities used to define stress measures.

The second-rank and symmetric Cauchy stress is defined using the measure of the elementary internal force in the material point t P according to:

t t dF dS     (2-13)
In large strains, numerous definitions for the stress tensors can be found in the literature [Khan95, Mase99, Saanouni12, Voyiadjis15], among others. Unfortunately, most of them do not have physical sense, these various stress tensors can be easily expressed by each other and the transformation gradient. The relationship between the different stresses tensors are summarized in Table 2.1, readers can refer to [START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF] for more details.

The formulation of mechanical models consists of the definition of some relations (i.e. constitutive equation) between stress and strain tensors. The use of Eulerian tensors has the advantage to get simple relations but poses the problem of the objectivity when anisotropic media are considered. However, the Lagrangian tensors avoid the objectivity problem because they are defined with respect to the unchanged initial configuration but they induce complex relationships without physical meaning. The formulation in a rotating frame which consists of the use of stress and strain quantities having Eulerian eigenvalues and Lagrangian orientation and leads to avoid the problems of the pure Eulerian and Lagrangian formulations.

 Kirchhoff  Piola-Lagrange  Piola-Kirchhoff S Cauchy  / J  . / T F J  ( . . ) / T F S F J Kirchhoff  J . T F  . . T F S F Piola-Lagrange  1 .( ) T J F   1 .( ) T F   . F S Piola-Kirchhoff S 1 1 . .( ) T JF F    1 1 . .( ) T F F    1 . F  

Kinematics of elasto-viscoplastic decomposition

Based on the concept of intermediate configuration and in order to study the elastoplastic media numerous works have focus on the decomposition of the total gradient F into elastic e F and plastic p F parts, i.e.

. .

e v p e v p F F F QV F Q F      (2-14)
This kinematics is illustrated in Figure. By comparison with time independent elastoplasticity this last configuration shall vary in time in the context of viscoplasticity. This time dependency of this configuration can be related physically, to some recovery phenomena [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF]. However, these phenomena of recovery need long time to take place and induce significant effects [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF]. For metal forming processes the loading velocity is high enough so that the elastic unloading can be supposed time independent as for classical time independent elastoplasticity.

Accordingly, the kinematics defined for time independent elastoplasticity, can be extended to the elasto viscoplastic framework and the objectivity can be ensured by the use of the rotating frame formulation defined by the rotation tensor Q representing the orientation of the isoclinic triad with respect to the current configuration. The main idea is that all the tensors defined in the current configuration are transformed (mapped) by the rotation tensor Q leading to the use of the Eulerian eigenvalues and Lagrangian orientations of the isoclinic configuration [START_REF] Sidoroff | The geometrical concept of intermediate configuration and elastic-plastic finite strain[END_REF][START_REF] Dogui | Plasticité anisotrope en grandes déformations[END_REF][START_REF] Sidoroff | Some issues about anisotropic elastic-plastic models at finite strain[END_REF][START_REF] Badreddine | On non-associative anisotropic finite plasticity fully coupled with isotropic ductile damage for metal forming[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. So any tensor with upper bar ( )  is mapped to the isoclinic configuration by Q according to:

   

(second rank tensor)

: :

(fourth rank tensor)

T T T T Q T Q T Q Q T Q Q              (2-15)
In the above equations, e V is the right elastic stretch tensor.

D L W D D                 ( 2 -1 6 ) Q T v p W Q Q W W      (2-17) where 1 [ . ] vp vp vp S D F F   
is the viscoplastic strain rate. eJ   is the rotated Jaumann rate of the small elastic strain tensor. vp W is the viscoplastic spin rate which will rule the rotating frame evolution according to Eq.(2-17). Different possible choices for this tensor have been discussed in [START_REF] Badreddine | On non-associative anisotropic finite plasticity fully coupled with isotropic ductile damage for metal forming[END_REF], in the present work we assume a kinematical definition for this tensor giving ( ):

vp vp vp W K V D  where ( ) vp K V
is a fourth-rank tensor which is function of the viscoplastic stretch tensor vp V , if the corotational (or Jaumann) rotating frame is chosen, this leads to a general equation of type: ( ) 0 vp K V  (refer to [START_REF] Badreddine | On non-associative anisotropic finite plasticity fully coupled with isotropic ductile damage for metal forming[END_REF] for the different forms of ( ) vp K V depending in the rotating frame).

In the following parts, all the developed constitutive equations are made with respect to the locally rotated configurations vp C or C (Figure 2.4) and for the sake of conciseness the upper bar is not used.

Fundamental laws of thermodynamics

Any formulation of a behavior model shall be done within the framework of thermodynamic of irreversible processes with appropriate number of state variables pairs. So it is necessary to recall the basic concepts of thermodynamics of continuum mechanics since they are very important in governing representative state variables of material continuum in material modeling. In this section, the main conservation laws of the physics of material continuum are presented.

Mass conservation law

The first conservation law postulates that the mass of the domain is kept unchanged for every transformation, this law can be easily expressed in the following differential form:

0 x div u        (2-18)
where

x div u      
denotes the spatial divergence (with respect to Euler coordinates) of the velocity field.

Momentum conservation law or the principle of virtual power

The momentum balance is obtained when the sum of the virtual power of internal force int ( 2 -1 9 )

The local form momentum balance can be expressed by the following partial differential equation with associated Neumann boundary condition:

( ) 0 div f u in F n on                     (2-20)
where  is the Cauchy stress tensor, the surface forces vector and the body forces vector are F  and f  respectively, u   is the acceleration vector, and n  is the outward vector normal to the boundary surface  of the solid.

The first law of thermodynamics or energy conservation

The first law of thermodynamics, also called energy conservation law, states that the internal energy rate per unit volume of the isolated system with volume  , must be equal to the sum of the external stress power and the heat flux received by  . The local form of the first principle of thermodynamics can be mathematically expressed as below:

: () 0 D e div q          (2-21)
where the product of total strain rate tensor D and Cauchy stress tensor stands for the stress power per unit volume in the deformed configuration, e is the specific internal energy (per mass unit ),  is the internal or body heat source, and q  is the heat flux vector received by the system throughout its boundaries.

The second law of thermodynamics

The second law of thermodynamics states that the total entropy of an isolated system can never decrease over time. Therefore, the rate of the entropy production is always greater than or equal to the amount of heat received divided by the absolute temperature. Its local form expressed in the current configuration:

( ) 0 q s div T T        ( 2 -2 2 )
where s is the specific (per unit of mass) entropy and T is the absolute temperature.

The Clausius-Duhem inequality

The Clausius-Duhem inequality is the combination of the first and second laws of the thermodynamics stated above. After eliminating the internal body heat quantity, the new inequality called Clausius-Duhem inequality is expressed by: :

( ) ( ) 0 q D s T g r a dT T           ( 2 -2 3 )
where the specific free energy or Helmholtz free energy per unit mass  is defined by:

e Ts    ( 2 -2 4 )
This inequality can be used as a measure for the thermodynamic admissibility, which plays an important role in the formulation of the constitutive equations of continua in the framework of the thermodynamics of irreversible processes.

Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 25

Ductile damage

Physical aspects of ductile damage

Damage in metallic materials is associated with a physical process of creation and evolution of micro-voids and microcracks following the well-known three stages of microdefects nucleation, growth and coalescence.

Ductility is regarded as an intrinsic ability of materials to undergo a certain amount of plastic strains before the final fracture occurs. The fracture of ductile material occurs in areas where large inelastic strains strongly localize prior to the initiation of macroscopic cracks. The damage occurring under large plastic or viscoplastic strains is called ductile damage (different from brittle damage, creep damage or fatigue damage) which is frequently observed in metal forming failure. As illustrated in Figure 2.5, the representation of ductile damage evolution in different length-scale is established at the characteristic time of the deformation process. At the smallest length scale, damage is related to the same processes of inelastic strains. At macroscale, damage is represented as the progressive degradation of a material, which exhibits a decrease in material properties as stiffness and strength.

Classification of ductile damage models

In order to describe the ductile damage, many damage models can be found in literature. They can be classified into three types: (1) uncoupled failure criteria; (2) weakly coupled models; (3) fully coupled models often formulated in the framework of continuum damage mechanics (CDM).

(1) The uncoupled approaches are based on specific failure criteria written in terms of stress or strain invariants giving the final fracture when their critical values are reached [START_REF] Freudenthal | The Inelastic Behavior of Engineering Materials and Structures[END_REF][START_REF] Johnson | Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures[END_REF][START_REF] Bao | On fracture locus in the equivalent strain and stress triaxiality space[END_REF][START_REF] Ebnoether | Predicting ductile fracture of low carbon steel sheets: Stress-based versus mixed stress/strain-based Mohr-Coulomb model[END_REF]. The advantage of these approaches are their simplicity by using the so-called failure indicators (criteria) which are easy to implement in FE software. Due to their uncoupled nature, the presence of micro-defects has no effect on mechanical fields (i.e. no interactions between the microdefects and the material behavior or no coupling). So that these models do not consider stress and stiffness softening within the material caused by damage occurrence during deformation.

This leads to an easy to use calibration of the material parameters. However, their major weakness relies in application to complex loading paths and large plastic strains. Recently, the research team (2) In the second approach, the influence of ductile damage in the yield condition is taken into account only on the yield function, by a porosity fraction and failure is predicted to occur when the porosity reaches a given critical value. The elastic stiffness is then not affected by damage (no decrease of the stiffness due to the damage occurrence). The first damage model within this approach was proposed by Gurson [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]. Later, Tvergaard and Needleman [START_REF] Needleman | An analysis of ductile rupture in notched bars[END_REF][START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] extended the Gurson model to include the void coalescence mechanism, this is called GTN model).

The starting point of Gurson's theory is the microscopic idealization of porous metals as aggregates containing voids of simple geometric shapes embedded in a metallic matrix whose behavior is governed by a rigid plastic von Mises yield function. The damage variable is considered as the local voids volume fraction. Extensions of the GTN model based on micromechanical studies (e.g.

[Acharya00, Tvergaard04, Bonfoh04]) have been made. However, there are still some limitations of this kind of models, where a large number of material constants needs to be identified and their physical meaning clarified, making their calibration procedure very difficult. Readers can refer to the work by [START_REF] Cao | Modeling ductile damage for complex loading paths[END_REF] for more details.

(3)

. The CDM approach with local formulation is adopted in this work in order to describe the isotropic ductile damage occurrence in metallic materials formed by large anisotropic visco-plastic strain at elevated temperature.

The concept of effective state variables in CDM

Ductile damage occurrence results from a progressive deterioration process of each material point including nucleation, growth and coalescence of microvoids and microcracks. The macroscopic behavior of the material is highly affected by these micro-defects at a given point of the material. The typical force-displacement curve with ductile damage effect is given in Figure 2.7. Clearly, from this figure, the ductile damage has deep effect on both elasticity and plasticity behaviors as the ductile damage grows. In the CDM framework, the damage interacts strongly with all the thermo-mechanical fields according to the appropriate coupling theories [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF][START_REF] Lemaitre | Engineering Damage Mechanics[END_REF]Besson10 and Saanouni12].

According to the effective stress concept of Kachanov [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF], two configurations can be defined: (i) the real damaged configuration including the local micro-discontinuities due to the micro-cracks and microvoids, (ii) the continuum undamaged fictitious configuration, as shown in Figure 2.8. In the damaged configuration the RVE with the cross-section area A is subjected to the uniaxial tensile force: Since in the damaged configuration the cross-section area is highly discontinuous and the distribution of the micro-defects is unknown, it is not easy to define the state variables at any damaged area. To make this easier, a fictitious undamaged configuration (Figure 2.8b) fully continuous and free from any defect in which the state variables can be easily defined.

In this fictitious undamaged configuration, effective state variables can be defined as function of the damage variable based on appropriate equivalence assumptions. In the literature, various equivalence principles are defined, such as: strain equivalence, stress equivalence and energy equivalence [Besson10, Murakami 12].

W d W d d                        ( 2 -2 6 ) ( , ) 1/ 2( : ) 1/ 2( : ) ( ) ( ) , ( ) k k X W d X X W d X d (2-27) ( , ) 1/ 2( ) 1/ 2( ) ( ) ( ) , ( ) r r i i R W r d r R r R W r r d r R d               
(2-28)

There are many choices for the damage effect functions  e (d),  k (d)

Time-independent plasticity and viscoplasticity theories without damage effects

In order to model the complex material behaviors introduced in Chapter 1, many works have been made to describe these complex behaviors at various scales. Crystal plasticity models [Jain07, [START_REF] Argon | Strengthening mechanisms in crystal plasticity[END_REF]Proust09] are provided by relating the single crystal constitutive behavior with that of the overall aggregate basing on full field approaches or mean field approaches [START_REF] Lebensohn | Self-consistent modelling of the mechanical behavior of viscoplastic polycrystals incorporating intragranular field fluctuations[END_REF][START_REF] Besson | Continuum models of ductile fracture: a review[END_REF][START_REF] Perdahcıoğlu | Constitutive modeling of two phase materials using the mean field method for homogenization[END_REF] like the selfconsistent approach [START_REF] Lebensohn | A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys[END_REF][START_REF] Saanouni | Micromechanical modeling of low cycle fatigue under complex loadings-PartI. Theoretical formulation[END_REF][START_REF] Boudifa | A micromechanical model for inelastic ductile damage prediction in polycrystalline metals for metal forming[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF] which are widely used to model the scale transition from grains scale to the macroscale (the RVE scale). For example, the initial anisotropy and tension-compression asymmetry of magnesium alloys were investigated by many works [Agnew 01, Agnew05, Lou07, Guo15 and Zhou16] at the microscale. It is noted that the microscopic models are useful for capturing the mechanical behavior of these highly heterogeneous metals. However, to fulfill our aim of modeling macro-scale forming processes with acceptable computational efficiency, a macroscopic phenomenological approach is more suitable. The following sections are focused on the modeling of the material behaviors at macro-scale based on phenomenological approach.

Time-independent plasticity

Yield condition, hardening model and flow rule are the basic plasticity concepts in phenomenological modeling. In order to define the limits which determines the transition from reversible and irreversible strains, the yield criterion is needed. The flow rule is a link of the stress and plastic strain components. The hardening models describe the evolution of the yield function during the inelastic strains. In time independent plasticity, all the above aspects are not sensitive to the rate effect i.e. to the velocity of the applied loading path.

Yield function

The yield functions are usally formulated in the stress space as a convex scalar valued function of the stress invariants. Generally, the stress invariants-based yield functions are defined using the stress invariant 1 I , 2 J and 3 J expressed as follows:

1 1 1 2 2 3 3 ( ) (2-31)

I tr         (2-29) 2 2 2 2 1 1 2 2 2 2 3 3 3 3 1 1 1 2 2 3 1 3 1 : 2 J S S
where S is the deviatoric stress tensor, defined as

1 ( ) I 3 S t r    
, with I is the second-rank identity tensor.

The classical well-known von Mises yield criterion assumes that the plasticity of the metals occurs when the 2 J reaches a yield value shown as follows:

2 3 3 : 0 2 y y f J S S        ( 2 -3 2 )
Hill's criterion [START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF] represents an extension of von Mises criterion by using the generalization of invariants to anisotropy.

Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 31 : : 0

y f S H S     (2-33)
where H is the forth-rank tensor which defined by six anisotropic parameters (F, G, H, L, M and N) to define various orthotropies. This criterion widely used for his simplicity to describe the initial anisotropic flow, it cannot describe accurately the plastic flow under complex (non-monotonic and non-proportional) loading paths, including some SD (Strength Differential) effects.

By including the effect of the third stress invariant 3 J which enables proper relative weighting of all shear stresses, Drucker [START_REF] Drucker | Relation of experiments to mathematical theories of plasticity[END_REF] proposed a yield function capable to describe the experimental data points located between von Mises and Tresca yield surfaces:

3 2 6 2 3 0 y f J cJ      ( 2 -3 4 )
To predict the SD effects in the anisotropic models, Blarlat and his co-authors [START_REF] Cazacu | Generalization of Drucker's yield criterion to orthotropy[END_REF][START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF][START_REF] Cazacu | Orthotropic yield criterion for hexagonal closed packed metals[END_REF] have developed two kinds of yield functions: (1) by introducing the third stress invariant based on the Drucker's criterion [START_REF] Drucker | Relation of experiments to mathematical theories of plasticity[END_REF], named CB04 given by Eq.(2-35); (2) by introducing a new material parameter to control the asymmetry in tension and compression and with extension to anisotropy using a linear stress translation based on Balart's non-quadratic yield criterion, named CPB06 in Eq.(2-36).

3 2 6 2 3 ( ) ( ) 0 y f J c J        (2-35)       1 1 2 2 3 3 0 a a a y f k k k                ( 2 -3 6 )
where 1  , 2  and 3  are the principal values of

: C S  
, with C the operator of plastic anisotropy. The stress invariants in Eq.(2-34) are replaced by 2 1 : 2 J     and 3 det( ) J    . These two yield functions have been extended by others to describe the pronounced in-plane anisotropy [START_REF] Plunkett | Orthotropic yield criteria for description of the anisotropy in tension and compression of sheet metals[END_REF].

Yoon [START_REF] Yoon | Asymmetric yield function based on the stress invariants for pressure sensitive metals[END_REF] proposed an anisotropic and asymmetric yield function based on CB04 for pressure sensitive metals given by:

3/2 1/3 1 2 3 ( ) 0 y f bI J cJ         (2-37)
These models have been successfully implemented into finite element codes and validated through numerical simulations of metal forming processes for different materials [START_REF] Li | Anisotropic and asymmetrical yielding and its distorted evolution: modeling and applications[END_REF][START_REF] Barros | Study on the effect of tension-compression asymmetry on the cylindrical cup forming of an AA2090-T3 alloy[END_REF]. Other methods to account for initial plastic flow anisotropy based on non-quadratic criteria, consist of using linear transformation stress tensors in isotropic yield functions to take into account the orthotropic behavior, as proposed for aluminum sheet metals [Barlat89, [START_REF] Barlat | A six-component yield function for anisotropic materials[END_REF][START_REF] Karafillis | A general anisotropic yield criterion using bounds and a transformation weighting tensor[END_REF][START_REF] Barlat | Plane stress yield function for aluminum alloy sheets-part 1: theory[END_REF][START_REF] Bron | A yield function for anisotropic materials application to aluminum alloys[END_REF][START_REF] Aretz | Unconditionally convex yield functions for sheet metal forming based on linear stress deviator transformation[END_REF][START_REF] Grilo | On the modelling of complex kinematic hardening and nonquadratic anisotropic yield criteria at finite strains: application to sheet metal forming[END_REF].

Various orthotropic non-quadratic yield functions can also be found among others [START_REF] Yoshida | A model of large-strain cyclic plasticity describing the Bauschinger effect and workhardening stagnation[END_REF][START_REF] Banabic | An improved analytical description of orthotropy in metallic sheets[END_REF] Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 32 Banabic10].

Flow rule

It is also important to describe correctly the evolution of all the internal variables as well as the evolution of yield surfaces. Indeed, the appropriate flow rule and hardening models play an important role in constitutive modeling.

To define the evolutions of dissipative phenomena with several strongly nonlinear and fully coupled dissipative phenomena (plasticity, viscoplasticity, hardening, thermal exchange, damage, etc.), the concept of flow potential which has the same properties of convexity in the stress space with yield function is defined.

The flow potential F can be taken as the same with the yield function f ( F f  ), this case is called the associative plasticity theory leading to Associative Flow Rules (AFR). However, the general case of the Non-Associative Flow Rules (non-AFR) is defined by considering

F f  with / / f F        and even with / / f F        .
It is important to recall that the amount of plastic flow increment is governed, for a given stress state and history, by the yield function f , while the direction of flow is collinear to the plastic potential normal. In crystalline plasticity, the resolved shear stress on several intersecting slip systems determines the yield function, the activation of any slip system determines the flow potential. The choice of these two flow rules depends on whether the stress states have the same effect on the plastic flow. In practice, the use of the non-associative theory is more convenient to a better prediction of experimental results of a wide class of metallic materials with multiple dissipative phenomena. In this work, the non-associative plasticity is used.

Hardening models

The proportional expansion of the yield surface is described by the isotropic hardening, which is represented by a pair of scalar state variables in continuum mechanics framework. However, the sliding of the yield surface center is classically described by the kinematic hardening, which is represented by a couple of symmetric and deviatoric second-rank tensor state variables called back stress and its associated strain variable, as shown in Figure 2.9. The kinematic hardening plays an important role in describing the Bauschinger effect in reversal loading. Several linear [START_REF] Prager | A new methods of analyzing stresses and strains in work hardening plastic solids[END_REF][START_REF] Ziegler | A modification of Prager's hardening rule[END_REF] and non-linear [Armstrong66, Chaboche86, Yoshida02, Dettmer04] kinematic hardening model can be found in the literature. However, the texture evolution in polycrystals due to rotation of the atomic lattice in single grains during large inelastic strains can lead to complex macroscopic mechanical behaviors, and the complex interplay at the microscale of a polycrystal leads to an evolving macroscopic anisotropy of the yield surface. The anisotropic behavior evolves with the material texture change during large plastic strain under complex loading paths. This induced anisotropy cannot be characterized only by classical isotropic and kinematic hardening without changing the shape of the yield surface. Meanwhile, unusual hardening evolution at compression loading of Mg alloys also need to be described by specific constitutive model. The macroscopic modeling of these aspects are discussed in the following: [START_REF] Teodosiu | Microstructure in the continuum modelling of plastic anisotropy[END_REF] can capture the cross hardening effect, but it cannot describe the evolving distortion of the yield function. An alternative method, motivated from physical point of view, has been used based on an evolving plastic fourth-rank tensor to describe the distortion of the yield surface [START_REF] Voyiadjis | Anisotropic distortional yield model[END_REF][START_REF] Haddadi | Towards an accurate description of the anisotropic behaviour of sheet metals under large plastic deformations: modelling, numerical analysis and identification[END_REF]Feigenbaum07,[START_REF] Pietryga | A finite deformation model for evolving flow anisotropy with distortional hardening including experimental validation[END_REF]. Feigenbaum and Dafalias [Feigenbaum07] proposed a thermodynamicallyconsistent framework with strong coupling between kinematic hardening and distortional hardening, but this model could only strictly be complying with the second law of thermodynamics in the case of small strains. However, the subsequent works have been made to extend the model to large strains framework [START_REF] Feigenbaum | Directional distortional hardening at large plastic deformations[END_REF]. François model [START_REF] François | A plasticity model with yield surface distortion for non-proportional loading[END_REF] is the first work that includes distortional hardening using a distorted deviatoric stress tensor instead of the classical deviatoric stress tensor in the yield criterion. Within other authors as [START_REF] Shutov | A viscoplasticity model with an enhanced control of the yield surface distortion[END_REF], the distortion was described by introducing a so-called distortional back stress in the equivalent stress. Barlat and coworkers [Barlat11,[START_REF] Barlat | Extension of homogeneous anisotropic hardening model to cross-loading with latent effects[END_REF][START_REF] Barlat | Enhancements of homogenous anisotropic hardening model and application to mild and dualphase steels[END_REF][START_REF] Jeong | A comparative study between microand macro-mechanical constitutive models developed for complex loading scenarios[END_REF] proposed a Homogeneous Anisotropic Hardening (HAH) model to describe the plastic flow during strain path changes.

It can describe the yield surface distortion using an internal structure tensor in the yield stress. Based on these distortional hardening models, extensions by coupling with other physical phenomena have been made. [START_REF] Shi | On the macroscopic description of yield surface evolution by means of distortional hardening models: application to magnesium[END_REF] analyzed three distortional hardening models with application to Magnesium alloys. Bartels and Mosler [Bartels17] extended the distortional hardening model modified by [START_REF] Shi | On the thermodynamically consistent modeling of distortional hardening: A novel generalized framework[END_REF] to the thermomechanical coupled setting. Badreddine et al. [Badreddine17] extended the François [START_REF] François | A plasticity model with yield surface distortion for non-proportional loading[END_REF] model to the non-associative plasticity including the full coupling with isotropic ductile damage. Rokhgireh [START_REF] Rokhgireh | Application of a new distortional yield surface model in cyclic uniaxial and multiaxial loading[END_REF] proposed a new distortional yield surface model which can describe the egg-shaped distorted yield surface, and the model shows better prediction of uniaxial and multiaxial ratchetting for steels. Qin et al. [Qin18] combined distortional hardening with isotropic and kinematic hardening in their model for complex loading paths, in which a scalar variable was introduced to expand the yield surface only in the orthogonal direction associated to the dislocation microstructure orientation, in such a manner that the yield surface is distorted.

Shi and Mosler

The fundamentals of some distortional hardening models are briefly discussed in the following:

(1)François' model [START_REF] François | A plasticity model with yield surface distortion for non-proportional loading[END_REF] Francois' model [START_REF] François | A plasticity model with yield surface distortion for non-proportional loading[END_REF] The yield function has the following form:

( , , ) 0

d y f X R S X R        (2-38)
where d S is defined by: : 2 ( )

o o d sat y S S S S X X R     (2-39)
with sat X is a scalar parameter denoting the saturation value of the kinematic hardening, the orthogonal part

of d S to the kinematic hardening X defines 2 : o S X S S X X   .
(2) Teodosiu's model [START_REF] Teodosiu | Microstructure in the continuum modelling of plastic anisotropy[END_REF] The Teodosiu's model [START_REF] Teodosiu | Microstructure in the continuum modelling of plastic anisotropy[END_REF][START_REF] Haddadi | Towards an accurate description of the anisotropic behaviour of sheet metals under large plastic deformations: modelling, numerical analysis and identification[END_REF] is based on a modified Hill-type yield function according to:

( ): :( ) ( ) 0 y f S X H S X R f          (2-40)
where ( ) f  represents the contribution of the organized dislocation structures on the isotropic hardening.

The evolutions of isotropic and kinematic hardenings are given as follows:

( ) R s a t R C R R      and ( ) X s a t X C X n X      (2-41)
where R C and x C are isotropic and kinematic hardening modules, n is the current direction of the plastic strain rate tensor, sat R and sat X are saturation values of R and X depending on the fourth-rank tensor . There are two parts of dislocation structures contributing to the strength, one is associated with current slip systems given as : :

D n n   
, and the other is related to the latent slip systems denoted by

L D n n     
 , their evolution equations are postulated as: (3)Feigenbaum and Dafalias' model [Feigenbaum07] The distinguished character of Feigenbaum and Dafalias' model [Feigenbaum07] is that a fourth-rank tensor is used to describe the evolution of distortional hardening and derived all evolution equations on the basis of the fulfillment of the dissipation inequality.

( ) nl L S L L s a t L C          and   1 2 ( ) D S D s a t D D C g g              (2-
    2 0 : ( ( : ) ) : 0 d r H f S X H n X A S X R         (2-43)
The fourth-rank tensor d H defines the shape of the yield surface which evolves with the strain through the fourth-rank anisotropic tensor A, the evolution function of this tensor is given in Eq.(2-44). The term : r n X is responsible for the directionality of the distortion while 0 H is a constant fourth-rank tensor representing the initial plastic anisotropy.

  2 1 2 3 : 2 r r r A A S X n X n n A A               (2-44)
where 1 A and 2 A are material parameters.

(4) HAH model [START_REF] Barlat | An alternative to kinematic hardening in classical plasticity[END_REF] The homogeneous yield function-based anisotropic hardening (HAH) model [START_REF] Barlat | An alternative to kinematic hardening in classical plasticity[END_REF] describes a partial distortion of the yield surface under plastic loading which can describe the Bauschinger effect without the concept of kinematic hardening. This yield function f is given by:

  1 1/ 1 2 ˆ( , ) [
] ( ) : : : : ( ) 0 q q q s s s s p q h q p q q q p eq y eq y e q

f s f h s h s f h s h s

                                   
(2-45)

The yield function combines the stable component q  and the fluctuating component h  , ˆs h is the deviatoric microstructure tensor. The material parameters 1 f and 2 f can be represented by two new parameters 1 g and 2 g [START_REF] Barlat | An alternative to kinematic hardening in classical plasticity[END_REF], according to:

1 1 1 1 1 q q f g         and 1 2 2 1 1 q q f g        
(2-46) f s s

with
      .
-

Tension-compression asymmetry in hardening

In addition to the tension-compression asymmetry in yielding, Mg alloys also exhibit unusual hardening evolution due to the twinning effect [START_REF] Guo | On twinning and anisotropy in rolled Mg alloy AZ31 under uniaxial compression[END_REF], the flow stress-strain evolution has a sigmoidal shape [START_REF] Khan | Mechanical response and texture evolution of AZ31 alloy at large strains for different strain rates and temperatures[END_REF][START_REF] Kabirian | Plastic deformation behavior of a thermo-mechanically processed AZ31 magnesium alloy under a wide range of temperature and strain rate[END_REF], which have referred as hardening asymmetry [START_REF] Lee | Distortional hardening concept for modeling anisotropic/asymmetric plastic behavior of AZ31B magnesium alloy sheets[END_REF]. Lee et al. [START_REF] Lee | Constitutive modeling for anisotropic/asymmetric hardening behavior of magnesium alloy sheets[END_REF] extended the twosurface hardening model for Mg alloys to include the unusual hardening asymmetry, the gap function is used to control the gap between the two surfaces in such a manner that the hardening differential effect can be captured. (1)Two-surface model proposed by Lee et al. [START_REF] Lee | Constitutive modeling for anisotropic/asymmetric hardening behavior of magnesium alloy sheets[END_REF] There are two yield surfaces in this model: a loading surface inside and a bounding surface outside, as shown in Figure 2.10. The current stress on the loading surface is at point a and the corresponding stress on the bounding surface is at point A. The correspondence between these two points is defined by the common yield surface normal directions. The hardening rate is determined by the gap between the current and corresponding stresses. The two surfaces can be contacted at points with the same normal direction. The yield functions for the loading and bounding surfaces are written as:

    1 2 0 0 iso iso f loading surface f A bounding surface                     ( 2 -4 7 )
where  is the Cauchy stress andis the back stress of the loading surface,  and  are the stress and back stress of the bounding surface. iso  and iso  represent the size of the loading surface and bounding surface respectively. The corresponding stress  on the bounding surface shares the same normal direction with the current stress at the loading surface, since the shapes of these two surfaces are the same. (1 exp( ( ))) 1 exp( (( ) / c ))

S y S T y T U y U J R f J R J R                   ( 2 -4 8 )
(1 exp( ( ))) 1 exp( (( ) / c ))

(1 exp( ))

P T T T T T T P T T T P U D U U U U P U U U P S S S S a R Q b a R Q b R Q b                                      ( 2 -4 9 )
The subsequent work in the same lab have been done to take into account the temperature effect and distortional hardening [START_REF] Kim | A temperature-dependent elasto-plastic constitutive model for magnesium alloy AZ31 sheets[END_REF][START_REF] Lee | Cyclic behavior of AZ31B Mg: experiments and non-isothermal forming simulations[END_REF][START_REF] Lee | Distortional hardening concept for modeling anisotropic/asymmetric plastic behavior of AZ31B magnesium alloy sheets[END_REF] with some changes in the equivalent stress and slight modification of hardening corresponding to different deformation modes. The idea remains the same, which is to apply three different yield criteria at three different deformation modes. To achieve this, a deformation criterion determining deformation modes under the plane stress condition should always be defined previously [START_REF] Kim | A temperature-dependent elasto-plastic constitutive model for magnesium alloy AZ31 sheets[END_REF][START_REF] Lee | Cyclic behavior of AZ31B Mg: experiments and non-isothermal forming simulations[END_REF]. The slip-dominant mode which is equivalent to uniaxial tension, the twinningdominant mode which is equivalent to uniaxial compression, and the detwinnning-dominant mode which is equivalent to tension after compression are defined as follows:

                       (2-50)
where 1   and 2   are the major and minor in-plane principal strain increments respectively, p T  is the accumulated plastic strain at twinning-dominant mode. This three-surface model was modified by [START_REF] Muhammad | An elastoplastic constitutive model for evolving asymmetric/anisotropic hardening behavior of AZ31B and ZEK100 magnesium alloy sheets considering monotonic and reverse loading paths[END_REF] to include the anisotropic yield function of CPB06 for three different deformation modes.

Each of the three deformation modes has been assigned to both equivalent stress and isotropic hardening.

No location translation of the yield surfaces was considered (i.e. no kinematic hardening). Another disadvantage is that a large number of material parameters need to be identified, since the anisotropic parameters are evolving with the accumulated plastic strain.

(3) Steglich's model [START_REF] Steglich | Mechanism-based modelling of plastic deformation in magnesium alloys[END_REF] Steglich et al. [START_REF] Steglich | Mechanism-based modelling of plastic deformation in magnesium alloys[END_REF] proposed a phenomenological plasticity model in which the slip and twinning are treated separately. Two yield surfaces are applied in this model, an anisotropic yield criterion g f proposed by Barlat together with the isotropic and kinematic hardening law is applied to describe the dislocation induced deformation, the CPB06 criterion t f was chosen to describe the deformation by twinning. , here the operator of plastic anisotropy C  and C are different. The interaction between the glide and twinning mechanism is realized through the isotropic hardening evolution shown as follows:

        1 2 3 3 1 1 2 1 1 2 2 3 3 1 ( , ) 0 2 ( ) 0 g g g g t t t a a a a g g g t a a a t t t f f k k k                                             (2-51)
    1 1 2 2 1 1 2 2 ( , ) 1 exp( ) 1 exp( ) ( ) exp( ) 1 1 exp( ) g g t g g t g g g g g g t t t t t t t t t t t R H Q b Q b R H Q b Q b                                         ( 2 -5 2 )
This model accounts for the transition behavior from twining-dominated to crystallographic slip-dominated deformation. It is different from the two-surface model proposed by Lee et al. [START_REF] Lee | Constitutive modeling for anisotropic/asymmetric hardening behavior of magnesium alloy sheets[END_REF] and the three-surface model proposed by Nguyen et al. [START_REF] Nguyen | A practical constitutive model for AZ31B Mg alloy sheets with unusual stress-strain response[END_REF], it provides phenomenological descriptions which take into account the information on the micro-mechanisms. Since the "detwinning" mechanism is not included, the model is unable to represent reverse loading paths.

Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 40

Thermo-viscoplasticity (time-dependent plasticity)

Viscoplasticity potential

With the thermodynamic approach a viscoplastic potential vp  should be included. By using the effective state variables, various forms of viscoplastic potential with strain hardening and damage can be introduced [START_REF] Besson | Mécanique non linéaire des matériaux[END_REF]. Two viscoplastic potentials are often used, the first is the power function of the yield criterion f which is called Norton-Hoff potential as give below: ,,,,, 1

    1 ,
v m v vp v v f R X T K R X T m K             (2-53)
The second viscoplastic potential is a hyperbolic function of the yield criterion f with the following form:

    1 2 2
, , , , , , cosh

v v vp v f R X T R X T K K K           (2-54) where v K , v m , 1 v K and 2 v K
are the material parameters. The evolution equations can be obtained by deriving the viscoplastic potential vp  with respect to the state variables.

Thermo-mechanical coupling

Within the hot sheet metal forming processes, the temperature rise is mainly caused by the heating operation.

The complexity of the mechanical behavior including anisotropy and hardening of the material is also increased by the thermomechanical coupling. It bears emphasis that the anisotropic material parameters have a very high non-linear change with the corresponding temperatures, for instance, high order polynomial fitting functions was adopted as a compromise method to predict the nonlinear variation of material parameters with temperature in the work of [Abedrabbo06] which investigates the anisotropic yield behavior of AA3003-H111 alloy at elevated temperatures. Khan and his coauthors [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] established a strain rate and temperature dependent anisotropic yield criterion based on Khan-Huang-Liang (KHL) phenomenological model.

Note that the dissipation phenomena caused by inelastic strains leads to temperature increase inside the material. The experimental investigations in [START_REF] Bednarek | The heating rate impact on parameters characteristic of steel behaviour under fire conditions[END_REF] proved that not only the temperature itself but also the heating rate has a significant effect on carrying capacity of the structure. Many thermomechanical coupled elasto-viscoplasticity theories can be found in the literature [Simo92, Chaboche97, Voyiadjis04, Chaboche08, Saanouni12]. Lestriez [START_REF] Lestriez | Numerical prediction of ductile damage in metal forming processes including thermal effects[END_REF] focused on the 'strong' coupling between the thermal and the elasto-viscoplastic behavior. A detailed description of viscoplasticity was presented in [START_REF] Chaboche | A review of some plasticity and viscoplasticity constitutive theories[END_REF],

in which the necessity of temperature rate terms in the context of hardening rules was discussed. Ganczarski and Skrzypek [Ganczarski09] considered the temperature dependence of all material functions that characterize viscoplasticity and damage with the additional temperature rate terms in all evolution equations of thermodynamic conjugate forces. A thermodynamic framework for constitutive modeling of ratedependent materials was derived in [Al-Rub12]. In the work of [START_REF] Egner | On the full coupling between thermo-plasticity and thermo-damage in thermodynamic modeling of dissipative materials[END_REF], a thermo-mechanical coupled formulation with non-associative plasticity and non-associative damage is developed, in which not only accounting for temperature-induced softening but also damage-induced softening. Soyarslan and Bargmann

[Soyarslan16] formulated a consistent-thermodynamic framework for finite multiplicative thermo-plasticity coupled to damage along the same lines with Simo and Miehe [START_REF] Simo | Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation[END_REF]. Based on all these works, the proposed model is presented in the following section.

Formulation of the proposed fully coupled model

State variables and effective state variables

In this work, the model is formulated in the framework of the thermodynamics of irreversible processes with state variables [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF][START_REF] Skrzypek | Modeling of material damage and failure of structures: theory and applications[END_REF]. The observable variables and internal variables are introduced, as shown in Table 2.2. Two pairs of observable variables:

 ( , B ) for total strain tensor and stress tensor or ( , B ) if total incompressibility is accepted.

 ( , e T s ) for absolute temperature and specific entropy.

Five pairs of internal variables: (

 ( ,
1 e e d      , 1 d      ) (2-55) ( 1 d      , 1 X X d    ) (2-56) ( 1 r d r     , 1 R R d     ) (2-57)
The parameter  represents a coupling parameter governing different damage effects on the isotropic hardening [START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. In order to take into account the microcracks closure effect which consists of rigidity recovery as well as a lower damage rate under compression(negative stress) than that in tension(positive stress), Lemaitre and Ladeveze [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF] proposed a spectral decomposition of the elastic damage force to a negative (or compressive) part and positive (or tensile) part. This straightforward spectral decomposition is applied to both strain and stress tensors to get positive and negative parts [START_REF] Lemaitre | Engineering Damage Mechanics[END_REF][START_REF] Lemaitre | Mécanique des matériaux solides[END_REF][START_REF] Bouchard | An enhanced Lemaitre model formulation for materials processing damage computation[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. Any symmetric second-rank tensor T can be additively decomposed in unique way into positive and negative parts T T

Table 2.2 Internal variables and their associated variables

T    
, the trace of tensor T can also be decomposed into positive and negative parts : 1

: 1 : 1 T T T    
. The application to stress and strain tensors can be done as follows:

1 1 d h d               ( 2 -5 8 ) 1 1 e e e e d h d               (2-59)
Accordingly, the effect of the ductile damage on the elastic behavior described by the effective variables can be decomposed into deviatoric and hydrostatic parts as following:

                                                                          (2-60)
where eD  is the deviatoric part of the small elastic strain tensor e  , eH  is the hydrostatic part. S is the deviatoric part of Cauchy stress and H  is the hydrostatic part. The microcracks closure parameter h has a range from 0 to 1, when 0 h  the microcracks close as soon as the stress tensor is negative, while 1 h  , the microcracks closure effect is neglected (see discussion in [START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]). The disadvantage of this approach is that it includes some discontinuity of the potential functions due to the spectral decomposition and the effect of the microcracks closure is only applied to the elastic part of damage energy release rate. To avoid the problem of discontinuity in all the previous works [Saanouni12, Badreddine10,

this decomposition is only used to define the evolution of damage and the rigidity recovery is ignored, so that in the yield surface and plastic potential the spectral decomposition is ignored to conserve their continuities.

In this work, in order to avoid the problem induced by the spectral decomposition, and to include the effect of microcracks closure effect on all the considered phenomena, a new approach is proposed. The basic idea consists in considering a continuous function of parameter h depending on triaxiality in the following form:

1 1 ( ) tanh( ) 2 2 c c h h h h        (2-61) where 2 ( ) ( ) tr J    
is the triaxiality, c h and h  are two material parameters which control the evolution of h , the objective of this function is to ensure h =1 when  is greater than a certain positive value, meanwhile when  is smaller than a certain negative value, h should be equal to c h (smaller than 1). As illustrated in The definition of the effective state variables is rewritten in the following form:

1 ( ) e eH eH h d        and 1 ( ) e H H h d        ( 2 -6 2 ) 1 ( ) eD e h d       and 1 ( ) S S h d     (2-63) 1 ( ) h d       and 1 ( ) X X h d     (2-64) 1 ( ) r r h d r      and 1 ( ) r R R h d      (2-65)
The parameters e  and r  are coupling parameters governing different damage effects on the hydrostatic stress and isotropic hardening respectively.

State potential and state relations

The specific Helmholtz free energy     , , , , , , ,

e e r T d r T          
, defined in the fictive undamaged configuration as a convex function of strain-like state variables in the effective strain space and concave of the temperature, is taken as a state potential [START_REF] Germain | Mécanique, vols. I and II, Ellipses[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. It can be decomposed into a thermoelastic part the  and a thermo-inelastic part thin  considering the assumption that the inelastic strain and hardening do not affect the elastic properties of the material.

                    (2-66) 2 0 0 0 2 0 0 0 1 ( )( : ) ( )( : ) 3 ( ) ( )( :1) ( ) 2 2 1 (1 ( ) ) ( )( : ) (1 ( ) ) ( )( : ) 2 3 1 ( ) ( ) ( )( :1 ) ( ) 2 1 1 ( ) : ( ) 3

T T T T h d T h d T C h d T T T T T T

C T Q T                                                           2 2 1 1 (1 ( ) ) ( ) : (1 ( ) ) ( ) 3 2 r r h dCT h d QTr                        (2-67)
In these equations,

( ) e T

 and

( ) e T



are the Lame's constants (

  ( ) ( ) / (1 )(1 2 ) e T E T        and   ( ) ( ) / 2(1 ) e T E T     ), ( ) e T 
is the compressibility modules defined by the following equation ratio  are assumed to be constant values, assuming their insensitivity to the temperature. T and 0 T are the absolute temperature and its reference value. The state potential can be written using the effective damaged material properties as following:

( ) (3 ( ) 2 ( )) / 3 ( ) / (3(1 2 ))
2 0 0 0 2 1 ( )( : ) ( )( : ) 2 3 ( ) ( )( :1) ( ) 2 1 1 ( ) : ( ) 3 2 eH eH eD eD the e e eH e thin T T C T T T T T T C T Q T r                                        (2-68)
where the effective damaged material properties are expressed as: The combination between the first and second principle of thermodynamics supplies the Clausius-Duhem inequality which expresses a statement concerning the irreversibility of natural processes. In the local form, the Clausius-Duhem inequality can be written as follows:

  ( ) ( ) 1 ( ) e e T T h d       ,   ( ) ( ) 1 ( ) e e e T T h d        , 1 ( ) e h d            ( ) ( ) 1 ( ) C T C T h d     and   ( ) ( ) 1 ( ) r Q T Q T h d      .
: ( ) ( ) 0 e q D sT g r a dT T           (2-69)
Taking the time derivative of Helmholtz free energy defined by Eq. (2-65) gives

: :

e e r d T r d T                             
( 2 -7 0 ) Substituting Eq. (2-70) into Eq.(2-69), we obtain the following local form of Clausius-Duhem inequality:

( ): : ( ) ( : ) ( ) 0 e v p e e q D s T r d grad T T r d T                                    (2-71)
This inequality is very important in the formulation of constitutive equations. It gives the expression of natural irreversibility processes which involved the energy dissipation.

By cancelling some terms in this inequality, a classical thermodynamic state law for the overall stress-like variables can be easily obtained:

  0 (1 ( ) ) ( ) 3 1 ( ) ( ) 1 2(1 ( ) ) ( ) H e e
eH e e e

S eD e h d T h d T T T h d T

                            ( 2 -7 2 ) 0 0 3 1 ( ) ( :1) (T T ) e eH e e C s h d T T                 ( 2 -7 3 )
Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage

47 2 (1 ( ) ) ( ) 3 X h dC T          ( 2 -7 4 ) (1 ( ) ) ( ) r R h d QTr r         ( 2 -7 5 ) e r Y Y Y Y d          (2-76)   1 1 0 1 ( ) ( )( : ) ( ) ( ) : 2 3 ( ) ( ) (T T ) ( :1) 2 1 ( )
h d T h d                          (2-76a) 1 ( ) ( ) : 3 Y h C T      (2-76b) 1 2 1 ( ) ( ) 2 r r r Y h d QT r     
(2-76c)

Dissipation analysis

It is assumed that the mechanical and the thermal dissipation are separately positive, so the dissipation analysis can be split into two parts: mechanical (or intrinsic) part M  and the thermal part T  .

0 M T       (2-77) ( ) 0 T q grad T T        (2-78) : : vp M D X Rr Yd          
( 2 -7 9 )

2.6.3.1Thermal dissipation analysis

The heat flux vector q  can be obtained from Fourier potential using the classical linear heat theory [START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF] in which k is the heat conduction coefficient.

( ) q k grad T     ( 2 -8 0 )
The generalized heat equation can be obtained by using this equation in conjunction with the first law of thermodynamics [START_REF] Lestriez | Modélisation numérique du couplage thermo-mécaniqueendommagement en transformations finies[END_REF][START_REF] Saanouni | Computational damage mechanics: Application to metal forming simulation[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. It will be used for solving the weak variational function associated with the thermal problem.

                                           
(2-81)

Mechanical dissipation

In order to define the flux variables ( , , , )

vp D r d    
using the non-associative plasticity theory, the intrinsic dissipation should satisfy the residual inequality defined in Eq.(2-77). For this end a yield criterion ,,,, )

( , , , , ) f X R d T  and dissipation potential ( ,
F X R Y d T 
are defined as follows:

  assumed to be function of the normalized Lode angle L  and equivalent plastic strain p according to:

1/3 3 3 ( ) 0 1 ( ) 1 ( ) r c c H J y Z w Z R f T h d h d            (2-82)                  1/3 3 3 ( ) 1 2 0 ( ) 3 ( : ) 4 1 ( ) 1 ( ) 1 ( ) , ( ) , 2 1 ( )
( ) 1 1 ( ) r r p p H J s T L T L Z w Z R aX X F C T h d h d h d S T Y Y T bR S T Q T h d s T h d                         (2-83) where 0 ( ), ( ), ( ), ( ), ( ) r S
1 1 2 2 3 2 1 ( , ) ( ) [1 ( )][ ( ) ( ( ) ( )) ( ( ( )))] 2 L L s a p a T Tanh a T a T a T Tanh p p T           (2-84) 1 1 2 2 3 2 1 ( , ) ( ) [1 ( )][ ( ) ( ( ) ( )) ( ( ( )))] 2 L L s b p b T Tanh b T b T b T Tanh p p T           (2-85) where 1 2 3 1 2 3 1 2
( ), ( ), ( ), ( ), ( ), ( ), , , ( )

s a T a T a T b T b T b T p T  
are the material parameters, p is the equivalent plastic strain, the normalized Lode angle L  ( 1

1 L     ) is given by: 3 3/ 2 2 3 3 6 2 1 1 arccos( ) 2 L L J J        
(2-86)

The unusual stress-strain evolution appears only in compression, while in tension the behavior is normal.

The idea is that the hardening parameters a and b shall keep constant in tensile loading while their values evolve under compressive loading. The variation tendency of the two parameters in compression takes very

T T T Q S Q S Q Q S Q QQ S S        (2-87)
where Q S is the rotated deviator of the stress tensor. A detailed parametric study of these material parameters is performed in Chapter 4.

  ( , ) ( ) ( ) ( ) Tanh( ) s L s h t e n s h L S T S T S T S T      
(2-88)

The parameter ( )

i i i d i d H Z S X H T S X i c p     ( 2 -8 9 )       3 det det ( ) : ( ) ; ( , ) i i i i d J Z Z H T S X i c p       ( 2 -9 0 )
where H and ' H are two anisotropic fourth-rank symmetric operators characterized by six material constants respectively, all these parameters are considered to be varying with temperature. In the nonassociative theory, each of these two operators is independently used in the criterion (Eq. 2-82) and in the flow potential (Eq. 2-83). Note that by taking 1 w the equivalent stress of [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF] is obtained, by taking 0 w  the orthotropic Hill's quadratic equivalent stress is recovered, meanwhile, the equivalent stress can be restricted to classical von Mises one when 0 w  and

3 3 1 2 2 3 1 1) ( dev H I I     . 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 2 i i i i i i i i i i i i i i i i G H H G H H F F G F F H H N M L                               and ' ' ' ' 0 0 0 ' ' ' ' 0 0 0 ' ' ' ' 0 0 0 0 0 0 2 ' 0 0 0 0 0 0 2 ' 0 0 0 0 0 0 2 ' i i i i i i i i i i i i i i G H H G H H F F G F F H H N M L                                (2-91)
Following the work of François [START_REF] François | A plasticity model with yield surface distortion for non-proportional loading[END_REF], a simple way to describe the distortion hardening within the classical thermodynamics framework is using the kinematic hardening variable to control the distorted yield surface. In this work, the François model is modified the by using three adequate material parameters to control the yield surface distortion [START_REF] Yue | Ductile damage prediction in sheet metal forming processes[END_REF]Badreddine17]. The classical deviatoric stress tensor S is replaced by a 'distortion stress' d S as shown hereafter:

    0 0 1 0 2 : 2(1 ( ) ) ( ) / 1 ( ) ( ) : 2(1 ( ) ) ( ) / 1 ( ) ( ) r r c d c l y l y S S S S X h d X T R h d T X X S h d X T R h d T                  (2-92)   0 0 1 : 2(1 ( ) ) ( ) / 1 ( ) ( ) r p d p l y S S S S X h d X T R h d T         
(2-93) 0 : :

x x S X S S S and S X X X    
(2-94)

The distortion hardening parameters 1 1 2 ( ), ( ), ( )

c p c l l l
X T X T X T are dependent on temperature. 1 ( ) c l X T and 1 ( )

p l X T
control the distortion ratio of the yield surface and plastic potential, 2 ( ) c l X T affect the size of the yield surface in the direction orthogonal to loading directions. It bears to emphasize that the distortional hardening parameters should be greater than the saturation value of the kinematic hardening (i.e. ,

1 s a t ( ) X C/ a c p l X T   )
to keep the convexity of the yield surface [Yue14, Badreddine17, Zhang18].

In the context of viscoplasticity, the following fluxes variables which define the evolution of the dissipative phenomena can be postulated as follows [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF][START_REF]Mechanics of anisotropic materials[END_REF]:

Eulerian viscoplastic strain rate tensor:

vp p f f F D n          (2-95)
Kinematic hardening strain rate tensor:

( ( , ) ) x f f F n a p X              (2-96) Isotropic hardening strain rate:   ( , ) r f f F r n b p r R            (2-97)
Isotropic ductile damage rate:

( ) 0 ( ) ( ) ( , ) (1 ( ) ) s T f f T Y Y T F d Y S T h d                       (2-98) with   0 1 1 : 1 ( ) (1 ( ) ) ( ) / 1 ( ) ( ) r p D A p l y X S n n I h d h d X T R h d T                                
(2-99)

Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 52

    0 0 1 0 1 ( : ) 2(1 ( ) ) ( ) / 1 ( ) ( ) 1 ( : )( : ) 1 ( ) (1 ( ) ) ( ) / 1 ( ) ( ) ( : ) r r A A p l y x A p l y S S n n h d X T R h d T n S X X n h d S h d X T R h d T X X                                    (2-100)   3 2 0 0 2 1 ( : )( : ) 1 1 1 ( ) 2(1 ( ) ) ( ) / 1 ( ) ( ) r r A r p l y S S n X n h d h d X T R h d T                       (2-101)   3 2/3 3 3 ( ): : ( ) 3 T H J A H J w Z H T Z Z Z H T n Z w Z           (2-102)
where the amount of viscous strain rate f   can be defined with the following Norton-Hoff form

[Lemaitre09]:   ( ) v n T f v f K T    (2-103)
The positive scalar f   can be assimilated to a viscoplastic 'multiplier' similarly to the time independent plasticity. But this viscoplastic 'multiplier' is not defined by a consistency condition like in time independent plasticity because the yield function given by Eq.(2-82) can be greater than zero for the context of viscoplasticity.

From Eq.(2-103) we can deduce the following viscoplastic "yield function" given by:

  1 ( ) ( , , , , , ) ( , , , , ) ( ) ( , , , , ) 0
v vp v n T v f v f X R d T f X R d T K T f X R d T             (2-104)
where the viscous stress

  1 ( ) ( ) v v n T v f K T    
is defined by a scalar quantity where ( )

v K T and ( ) v n T are
material parameters characterizing the viscous effect. Eq.(2-104) is the yield function in viscoplasticity and will be solved numerically at each integration point of each finite element in next chapter.

Thermodynamical admissibility

The thermodynamical admissibility consists in proving the positivity of the Clausius-Duhem inequality, the viscoplastic intrinsic dissipation is expressed as:

 : : m d M v p D X Rr Yd               
(2-105)

Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 53 With the evolution equations in Eqs. (2-95) to (2-98), the inelastic dissipation can be rewritten as below:

            0 : :( ) ( ) 0 (1 ( ) ) s T p x r f T Y Y T Y n X n a T R n b T r S T h d                       (2-106) Since 0 f   
, the equation becomes:

            0 : :( ) ( ) 0 (1 ( ) ) s T p x r T Y Y T Y n X n a T R n b T r S T h d             (2-107) Since O S S 
especially when the loading path is linear at least by part, we assume that 1 ( )

A x p n n n h d     and 1 1 ( ) r r n h d    
. The equation above can be rewritten as:

            0 ( ): : 0 1 ( ) (1 ( ) ) 1 ( ) r s T A T Y Y T X n R Y a T X b T rR S T h d h d h d                  (2-108) Let Z X   
, the first term in Eq. (2-108) can be written as:

  1/3 3 3 3 2/3 3 3 : ( ): :( ( ): ) ( ): 3 1 ( ) 1 ( ) 1 ( ) T A H J H J H J w Z Z H T Z Z Z H T Z Z w Z X n h d h d h d Z w Z                         (2-109)
Similarly, with Eq. (2-109) in state relations, we can obtain the following form:

                1/3 3 2 3 0 (1 ( ) ) 1 1 ( ) 3 : 0 2 ( 1 ( ) ) (1 ( ) ) r r H J s T T Z w Z b T R R Q T h d d h d a T Y Y T X X Y C T h d S T h d                       
(2-110) By combining Eq. with Eq.(2-110), we can simplify the equation as follows:

                    1/3 3 3 1 0 0 2 0 2 0 3 0 4 0 1 ( ) 1 ( ) 3 : 2 ( 1 ( ) ) (1 ( ) ) (1 ( ) ) r vp r H J y v y v Part f s T T Part Part Part Z w Z R T T h d h d b T a T Y Y T R XX Y Q T C T h d S T h d h d                                                        0   (2-111)
As shown in Eq.(2-111), all the parts are positive if b(T) has the same sign as Q(T) and a(T) have the same Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 54 sigh as C(T). This means that the constitutive equations defined by the state relations and the evolution equations fulfill the Clausius-Duhem inequality which is the combination of the first and second principle of thermodynamics. The thermodynamic consistency of the proposed model is thus proved.

Contact and friction

The contact and friction between tools and pieces have an important influence on the metal forming process through large inelastic strains. In order to obtain accurate simulation results, the contact-friction phenomenon must be treated carefully in the simulation. In this thesis, we will not propose any new model of contact-friction, however, a general introduction of the common used friction model in FE code Abaqus © is indispensable. Detailed study on modeling of contact-friction phenomenon can be found in literature [START_REF] Shillor | Models and Analysis of Quasistatic Contact: Variational Methods[END_REF][START_REF] Hashiguchi | Elastoplasticity Theory[END_REF][START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]. More deep-going study on this topic is still needed.

Kinematics of contact

Generally, the contact problem is always considered between a pair of solids, as shown in Figure 2 To model the contact between the two solids A and B, the vector related to displacement u  and the vector of force F  are defined as shown in Figure 2.14. The displacement and the force of the contact are given as follows: The unilateral contact is commonly used in the finite element calculation. This system of solids in contact must verify the following additional conditions:

(A) (A) (A) (A) (A) n t u u n u t        (2-112) (A) (A) (A) (A) (A) n t F F n F t        (2-113)
 Impenetrability condition: (A) ( ) 0 n n u u B    Contact surface non-adhesion condition:   (A) ( ) (A) 0 n n n u u B F    Non-contact condition: 0 n u  and 0 n F 

Friction model

Friction models are constitutive equations for contact interfaces between solids under the framework of thermodynamics of irreversible process. Two commonly used friction model are presented below.

Tresca-type friction model

The Tresca-type friction model is a time-independent friction model which only has one friction constant, the sliding limit is not dependent on the normal stress. This model is given as follows:

(A) 0 3 s T g F F v      (2-114) (A) 0, (A) (A) 3 s T g T F F v F            (2-

Coulomb-type friction model

The Coulomb-type friction model is often applied to describe friction between bodies in sliding contact in the case of low contact pressure. It is not dependent on time. The model can be expressed as follows:

Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage

56 (A) (A) 0 T n F F      (2-116) (A) (A) (A) 0 T n g F F v        (2-117) (A) (A) 0, (A) (A) T n g T F F v F             (2-118)
where  is the coulomb friction coefficient, the value of  is often chosen between0.01(with thin films lubrication) and 0.5(without lubrication). It should be noted that Coulomb friction model is appropriate only when the mean contact pressure between the two contact bodies lies below the flow stress of the softer body.

The Coulomb-type friction model is illustrated in Figure 2.16. For the velocity:

1 2 1 2 1 ( ) ( )e v V V V                  (2-119)
For the temperature:

  2 1 0 0 1 / ( ) 1 t t t f T T T T e                       (2-120) where V v   is the slip velocity, T is the current temperature, Tf is the melting temperature 1 2 1 2
, , , , ,

t v t t
      and 0  are material parameters. To take into account the effect of temperature and velocity at the same time, the dependent functions for velocity and temperature can be composed together by multiplication as follows:

1 2 0 ( , ) ( )(1 ) t t T V T V e           (2-121)
In a similar way, some other influencing factors can also be considered, i.e. a friction model considering the contact pressure, velocity and the temperature can be found in [START_REF] Klocke | FE-analysis and in situ visualization of pressure-, slip-rate-, and temperature-dependent coefficients of friction for advanced sheet metal forming: development of a novel coupled user subroutine for shell and continuum discretization[END_REF].

Heat generation from friction

During the frictional processes of two solids, the transformation of mechanical energy into internal energy or heat is often observed, leading to increase of the temperature [START_REF] Kennedy | Frictional heating and contact temperatures[END_REF]. This amount of heat generated at the contact interfaces is distributed in the two solids in contact, by conduction, convection and/or radiation.

Frictional heating and the resulting contact temperatures have an important effect on the contacted solids [START_REF] Martins | Numerical modeling of the thermal contact in metal forming processes[END_REF]. The contact interfaces temperatures can increase high enough to cause changes in the mechanical properties of the sliding materials, oxidation of the surface, and possibly even melting of the contacting solids. The heat generation represented by a general two-body sliding contact is shown in Figure 2.17, the "fine" description of the topology of surfaces showing that the contact is not perfectly flat but made up of numerous "pads" of contact between roughness of contacting surfaces, the friction force and the relative sliding velocity determines the rate of total energy dissipated in the sliding contact. The total heat t q generated from the friction is given by:

t q PV FV    (2-122)
Where  is the coefficient of friction, F is the friction force, P is the contact pressure, Vis the relative sliding velocity. In order to solve the thermal problem (Eq. 2-81) which includes the heat generation (Eq.2-119) at the contact interface, the suitable thermal boundary conditions applied to different operating conditions and geometry of the contacting solid bodies are needed. For example, the Dirichlet condition specifies that a temperature

Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 58 is imposed on one part of the boundary of each of the solids outside the contact interfaces, and the Neumann condition, specifying that a heat flow can be imposed on other parts of the solid boundaries outside the contact interfaces (see [START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF] for further details).

Conclusions

In this chapter, the first two sections (Section 2.2 and 2.3) illustrate the fundamentals of finite strain elastoviscoplasticity theory and some basic elements of mechanics and thermodynamics of the materially simple continua to obtain a better understanding the framework in which the whole model is formulated. The physical aspects of the ductile damage and the different ductile damage models are described in Section 2.4, followed by a summary of the framework of continuum damage mechanics. A review of time-dependent and time-independent plasticity are discussed in Section 2.5. The formulation of the proposed anisotropic thermo-elaso-visco-plasticity model fully coupled with isotropic ductile damage is given in Section 2.6. The proposed fully coupled model, which is under the framework of thermodynamics of irreversible processe with state variables, embraces many phenomena observed in metal inelastic flow (initial and induced anisotropies, SD effect, hardening asymmetry, etc.). Especially, in this section a new formulation of the microcracks closure effect is developed and discussed. Finally, a short discussion of the contact and friction in metal forming are given.

In the next chapter (Chapter 3) the fully coupled constitutive equations developed in the present chapter will be discretized in space and in time and implemented in the general purpose finite element code and used for the simulation of various metal forming processes.

Introduction

The numerical aspects for solving the equilibrium problems of the fully coupled thermo-mechanical equations by large inelastic strains with damage effect is elaborated in this chapter. First, the strong and weak forms of the initial and boundary value problem (IBVP) are defined in section 3.2.1. Then the time and space discretization of the IBVP leading to obtain the highly nonlinear and strongly coupled algebraic system, as given in section 3.2.2. The method to solve the IBVP is the global resolution scheme of the dynamic explicit solver and the local constitutive integration scheme. A fairly detailed presentation of the numerical integration of fully coupled ordinary differential equations (ODEs) associated with thermoelastoviscoplastic constitutive equations with damage are illustrated in section 3.2.3 and 3.2.4. Finally the implementation of a user defined subroutine-VUMAT is given in section 3.3.

Numerical method for thermo-mechanical problems

Thermo-mechanical initial and boundary value problem

The mechanical and thermal problems are posed with the strong and weak forms of the equilibrium equations and the heat equations. As illustrated in Figure 3.1, a deformable and damageable solid occupying The following force fields are imposed:

-Force field  F F  imposed on F t  ; Chapter 3: Numerical aspects 61 -Contact force C F  imposed in C t  ; -Displacement field t U   imposed on u t  . -A temperature field t T  imposed on T t  -A heat flow t q   imposed on q t 
The boundary t  is composed on the following sub-surfaces

u t  , F t  , q t  , T t  and C t 
, for which the following classical relationships are fulfilled at each time t:

u F C q T t t t t t t F u F C u C q T t t t t t t t t                                (3-1)
All the kinematics and state variables must fulfill all of the field equations governing the equilibrium of forces and temperature.

For the mechanical problem, the strong form is defined by the classical equilibrium equations together with mixed Dirichlet and Neumann Boundary conditions:

( ( , )) ( , ) ( , )
,

t t t div x t f x t U x t x t I                (3-2) ( , ) on ( , ) on u t t F F t t t U x t U x t n F                    (3-3)
where  is the Cauchy stress tensor defined by the fully coupled constitutive equations, t  is the solid density and U   is the acceleration vector.

The weak form of the IBVP is deduced from the strong form defined above by Eq. (3-1) to Eq. (3-3), thanks to the well-known weighted residual method together with Galerkin assumption. If the updated Lagrangian formulation is used, the weighted residual method applied to Eq. (3-1) and after the integration by part and the use of the Neumann Boundary condition, the following weak form is obtained:

 ( , , ) : 0 F F C V t J U U T D dV U UdV f UdV F F U ds                                               . . U K A     (3-4)
where U    is the kinematically admissible (K.A.) virtual velocity field and  D is the associated virtual total strain rate tensor.

For the thermal problem on the same volume t  , the heat equation is written in the case of single-surface and isotropic damage theory as follows:

( ( , )) ( , ) ( , ) ( , ) 0 , v t t k Lap T x t C T x t x t Rpl x t x t I                 (3-5)
The appropriated mixed Dirichlet and Neumann boundary conditions should be added after the discretization of Eq. (3-5) by the FEM.

    ( , ) ( , ) .
( , ) .

T t t q t t t t T x t T on q x t n k T x t n q on                      (3-6)
In Eq. (3-5) ( , )

x t  
is the body heat source, while pl R is the internal heat generation from mechanical dissipation which is the sum of the intrinsic volume dissipation and the isentropic term with mixed hardening and isotropic damage, defined by: : : : :

in e pl J X R Y R D X Rr Yd T r d T T T T                                   (3-7)
When considering the contact with friction between different solids, the heat generated at the contact interfaces cannot be underestimated, since this could cause a significant rise of temperature in both solids, the contact interfaces behave as a heat source, which can distribute the heat by conduction, convection and/or radiation. The density of the heat flow comes from the different heat fluxes generated at the contact interface are expressed in the following form:

c r s q q q q        (3-8)
where ( )

B A c q h T T    and 4 4 (( ) ( ) ) B A r r q h T T    
are the heat flux density form the solid by conduction and radiation respectively. The vector s q  is the density of the heat flux generated by friction at the contact interfaces.

The weak form associated with heat equation can be obtained by applying the weighted residuals method to heat equation in Eq. (3-5):

( , , ) ( ) 0 t t t t v p l J T T U T CTdV TR dV T dV Tdiv q dV                       . . T K A   (3-9)
With the Fourier model

( ) q kg T    
, the divergence theorem

( ) t S div Tq dV TqndS         it comes: ( ) ( ) t t S Tdiv q dV TqndS qdiv T dV               (3-10)
Then the weak form of the heat equation is written by: ( , , ) ( ) 0

t t t t v p l S J T T U T C TdV TR dV T dV TqndS qdiv T dV                           . . T K A   (3-11)
Chapter 3: Numerical aspects 63 Taking into account the thermal boundary conditions imposed on T t  and q t  , the final form of the associated weak form is finally expressed as:

( , , ) ( ). ( ) 0 . . t t t t T q t t v p l f T q J T T U T C TdV TR dV T dV kg T g T dV Tq n dS Tqn dS T K A                                     (3-12)

Time and space discretization of IBVP

The two weak forms defining the thermo-mechanical IBVP discussed above are strongly coupled and highly nonlinear [START_REF] Bergheau | Finite element simulation of heat transfer[END_REF]. In order to linearize these variational forms, the discretization of them in time and space are needed, which can obtain the IBVP in an incremental form of a nonlinear algebraic system to be solved over each time increment.

Space discretization of the IBVP

Using the FE method based on the displacement and temperature fields, the discretization of the total volume t  of the solid into e N finite elements with the elementary volume e  leading to:

1 e N e t e     (3-13)
To each finite element e  with e N nodes defined in the working space   1 2 3 , ,

x x x is associated a reference element r  defined in the reference space  

1 2 3
, ,    , so that:

        e e i i x N x             (3-14)
where    

x   is the coordinates vector of any material point of r  ,

 

e i N       



are the shape polynomial functions and   e i

x are the coordinates of the e N nodes of the element.

By applying this type of nodal approximation, together with the Galerkin assumption, the real and the virtual displacement vector and temperature are approximated by:

        ( ) ( ) and ( ) ( ) e e e e e e u i u i U U N U U N U                     (3-15)     ( ) ( ) and
    ( ) ( ) e e e I I U N U            and     ( ) ( ) e e e I I U N U             (3-17)   ( ) ( ) e e e T i T N T         (3-18)
where   e I U  is the accelerations vector of the element nodes.

By using the equations above, the weak form associated to the mechanical and thermal problem written for a single iso-parametric (i.e.

( ) ( )

e e I I N N             
) reference element is easily expressed with the following matrix form:

                    e e e e ext J U U T U M U F F J T T U T C T G G                               ( 3 -1 9 )
For the mechanical problem, the internal element force vector and the external forces vectors are defined for a typical element under the following form:

            int e e F C T e e e V
F B J d V F N F J dV N F J dA N F J dA                                   ( 3 -2 0 )
For the thermal problem, the internal and external thermal force vectors are given for a typical element by:

       
V G N R d V k B B T d V G N d V N qd N q d                                        (3-21)
where e M     is the consistent mass matrix given by:

r T e e e M N N d V                 (3-22)
The capacitance matrix is defined by:

  r e e e T v T C N C N d V         (3-23)
                                     ( 3 -2 4 ) e e e T T T e e N N B x x                                      (3-25)
If  indicates the finite element assembly operator, the discretized weak form of Eq. (3-19) can be written as:

                          int 1 int 1 , , 0 . . , , 0 .
.

J J U U T U M U F F U K A J J T T U T C T G G T K A                                    (3-26)
Leading to:

                  int int 0 0 ext ext M U F F C T G G              ( 3 -2 7 )
where   U  is the global (for all nodes of all of the structure elements) accelerations vector, and   1

Time discretization of the IBVP

The discretization in time of the IBVP consist of dividing the total time into subintervals with empty intersections, which leads to the approximation

  0 1 0 , , t N f n n n n t t t t t t         
. The unknowns of the IBVP are determined by solving the non-linear problem at each time increment. When all the unknowns are assumed to be known at time n t , this problem is to obtain the unknowns at the end of each time increment at time 1 n t  .

Global resolution scheme: computation of displacement and temperature

When dealing with the numerical simulations in metal forming processes, the strong multi-physical coupling between the various partial differential equations (PDEs) is always taken place. Regarding the mathematic Chapter 3: Numerical aspects 66 nature of the problems, there are two frequently used global resolution schemes: dynamic explicit (DE) schemes and static implicit (SI) schemes. They are applied to solve the highly nonlinear algebraic system over each typical time increment

1 n n t t t     .
The static implicit schemes are usually used for solving the static (or quasi-static) problems for which an incremental tangent matrix can be relatively easily computed.

When concerning the damage-induced softening as well as a very large number of contact nodes, it is difficult to obtain an incremental tangent matrix due to strong nonlinearities and physical instabilities induced by the damage effects. In that case, the dynamic explicit scheme is preferred to solve the IBVP using the lumped mass and capacitance matrices [Saanouni 12]. Due to the fact that the damage is fully coupled in our constitutive model, the DE resolution scheme is applied in this work to perform the simulation of sheet metal forming process.

For the coupled thermo-mechanical problem, the DE method is used based on the mass matrix   M and capacitance matrix   C .This method has been implemented into the FE code ABAQUS/Explicit to calculate the unknowns at 1 n t  which are the functions of the values known at n t . For the mechanical problem, the acceleration vector can be obtained when the internal and external force vectors as well as the mass matrix are known as follows:

        1 int L ext n n n n U M F F         (3-28)
Once the acceleration vector is known at the beginning of the time increment, the displacement velocity at the middle of the time increment is approximated by the following Taylor development truncated at the first term:

        1 1/ 2 1/ 2 ... 2 n n n n n t t U U U             ( 3 -2 9 )
Then the displacement vector at the end of the increment 1 n t  is obtained using the same method:

      1 1 1/ 2 ... n n n n U U t U         (3-30)
The same method is applied to the heat equation by computing first the temperature rate at n t :

        1 int L ext n n n n T M G G         (3-31)
The temperature at time 1 n t  is then obtained by: );

      1 n n n n T T t T     (3-32)
(2) Estimate the stable time increment u t  using the Eq. (3-34);

(3) Calculate the displacement field (1) Calculate the lumped capacitance matrix [C ] e and the thermal force vectors at n t with knowing the results of mechanical variable at 1 n t  ;

(2) Estimate the stable time increment (6) Select the critical time step for the thermal problem.

 Numerical treatment of contact with friction

In the simulation of metal forming processes, attention should be paid to the numerical treatment of contact with friction, since the friction has an important influence on the numerical results. A standard way in Abaqus (see Abaqus Analysis User's Guide) to calculate the contact forces is to use the penalty method which calculates the interpenetration gap between solids at a given contact interface node, and uses this displacement and a penalty parameter to compute the contact forces.

Here the numerical treatment of contacting nodes based on master/slave surfaces algorithm is to ensure the proper contact conditions. The contact forces, function of the penetration distance, are applied to the slave nodes to oppose the penetration, while equal and opposite forces act on the master surface distributed at the penetration point. The Coulomb friction model describes the linear relationship between the transmitted shear and the normal forces through the friction parameter, it is used to represent the friction between the contacting bodies.

As illustrated above, the explicit resolution procedure for the fully coupled thermo-mechanical problem is accelerated, since the diagonal character of the lumped mass and capacitance matrices. Furthermore, the mass scaling procedure will be applied in the simulation with ABAQUS/Explicit, it can increase the mass of the elements with small size while preserving the real time scale, in such a manner that the size of the time step is increased, finally, the goal of reducing the CPU time is reached.

Local integration scheme: computation of state variables

In order to calculate the internal force vectors which characterizing the IBVP, the Cauchy stress tensor T  at each integration point at each element at time step 1 n t  need to be computed. The calculation of these Chapter 3: Numerical aspects 69 state variables is performed using the well-known elastic prediction and plastic correction scheme throughout the numerical integration of the ODEs used in constitutive equations. The Newton-Raphson implicit iterative integration method will be applied in this work to solve all the first-order ODEs in constitutive equations which are highly nonlinear and strongly coupled.

Usually, when the time domain is discretized into many small time intervals, as  

1 , n n n t t t t
   , the firstorder ODE integration problem can be expressed in the following form:

      1 1 1 , , n n n n n n y y t t t t y t t y                (3-35)
Suppose that the variable n y at the beginning of the time increment n t is known, the problem is to obtain its value 1 n y  at the end of the time increment 1 n t  . With the  -method, the solution can be written in the following form:

    1 1 ( 0 1 ) n n n n y y t y y for               (3-36)
The forward Euler purely implicit scheme can be set when 1

  , we will confine ourselves to use this scheme in the thesis. For the specific type of Euler first-order differential equations which need to be solved are given below:

        1 1 1 1 1 , ,
                       (3-37)
The asymptotic solution for Eq. (3-37) is given as:

          exp 1 exp (0 1) n n n n n y y y t y t y for                          (3-38)
In the purely implicit case when 1   , this solution becomes:

          1 1 1 1 exp 1 exp n n n n n y y y t y t y                  
(3-39)

Pose of the problem

The anisotropic thermo-elasto-viscoplasticity model fully coupled with isotropic ductile damage constitutive equations developed in Chapter 2 are summarized below:

Yield function   1/3 3 3 ( ) 0 1 ( ) 1 ( ) r c c H J y Z w Z R f T h d h d            Chapter
i i i d i d H i i i i d J Z S X H T S X i c p Z Z H T S X                   0 0 0 1 2 : : 2(1 ( ) ) ( ) / 1 ( ) ( ) 2(1 ( ) ) ( ) / 1 ( ) ( ) r r c d c l y l y S S X X S S X S h d X T R h d T h d X T R h d T                    0 0 1 : 2(1 ( ) ) ( ) / 1 ( ) ( ) r p d p l y S S S S X h d X T R h d T          0 : : x x S X S S S and S X X X    

State equations Evolution equations

Cauchy stress

  0 (1 ( ) ) ( ) 2(1 ( ) ) ( ) 3 1 ( ) ( ) 1 H
h d T T T                            1 f p p p f f F D n n d               Kinematic hardening 2 (1 ( ) ) ( ) 3 X h dC T          1 ( ( ) ) 1 p f f F n a T X d               Isotropic hardening (1 ( ) ) ( ) r R h d QTr r           ( ) r f f F r n b Tr R           Damage governing force e r Y Y Y Y d          0 ( ) ( , ) (1 ( ) ) s f f Y Y T F d Y S T h d                      
where:

  0 1 1 : 1 ( ) (1 ( ) ) ( ) / 1 ( ) ( ) r p D A p l y X S n n I h d h d X T R h d T                                     0 0 1 0 1 ( : ) 2(1 ( ) ) ( ) / 1 ( ) ( ) 1 ( : )( : ) 1 ( ) (1 ( ) ) ( ) / 1 ( ) ( ) ( : ) r r A A p l y x A p l y S S n n h d X T R h d T n S X X n h d S h d X T R h d T X X                                      3 2 0 0 2 1 ( : )( : ) 1 1 1 ( ) 2(1 ( ) ) ( ) / 1 ( ) ( ) r r A r p l y S S n X n h d h d X T R h d T                       Chapter 3: Numerical aspects 71   3 2/3 3 3 ( ): : ( ) 3 T H J A H J w Z H T Z Z Z H T n Z w Z          
These fully coupled constitutive equations presented above have been implemented into ABAQUS/Explicit FE code using the user defined subroutine VUMAT. Since the constitutive equations above have the general form of first-order ordinary differential equations except the kinematic hardening and isotropic hardening having the Euler-type differential equations, their numerical integration can be performed by combining the backward Euler fully implicit ( 1   ) scheme and the fully implicit asymptotic scheme (see [START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]).

The discretized equations can be easily obtained as follows:

  ( ) 1/3 3 3 1 1 1 1 1 1 ( 1 ) 1 ( 1) 1 1 ( , , , , , ) ( )
0 1 ( ) 1 ( ) T r n c c H J vp n n n n n n vn n n y n vn n n Z w Z R f X R d T T h d h d                           (3-40) with       ( 1) 1 ( 1) 1 ( 1) ( 1) 1 3( 1) ( ) : ( ) : ( ) ; ( , ) det det ( ) : ( )
i i i d n n i n d n n H n i i i i n d n n J n Z S X H T S X i c p Z Z HT S X                       (3-41)     0( 1) 0( 1) 1 ( 1) 1 1 1 1 1 1 1 0( 1) 1 2 1 1 : 2(1 ( ) ) ( ) / 1 ( ) ( ) : 2(1 ( ) ) ( ) / 1 ( ) ( ) r r n n c n d n n c n l n n n y n n n n n l n n n y n S S S S X h d X T R h d T X X S h d X T R h d T                                (3-42)   0( 1) 0( 1) 1 ( 1) 1 1 1 1 1 : 2(1 ( ) ) ( ) / 1 ( ) ( ) r n n p n d n n p n l n n n y n S S S S X h d X T R h d T                  (3-43) 1 1 1 0( 1) 1 1 1 : : n n n n n n n S X S S X X X          and 1 1 1 1 ( )1 3 n n n S t r        (3-44)
The stress-like variables are given by:

  1 1 1 1 1 1 1 1 0 (1 ( ) ) ( ) 2(1 ( ) ) ( ) 3

h d T h d T h d T T T I

                                 (3-45) 1 1 1 2 (1 ( ) ) ( ) 3 n n n n X h d CT        (3-46) 1 1 1 (1 ( ) ) ( ) r n n n n R h d QT r        (3-47) 1 1 1 1 e r n n n n Y Y Y Y         (3-48)   1 1 1 1 1 1 1 1 0 1 1 1 ( ) ( )( : ) ( ) ( ) : 2 3 ( ) ( ) ( T ) ( : ) 2 1 ( )
Y h d T h T h d T T I h d                                  (3-49) 1 1 1 1 ( ) ( ) : 3 n n n n Y h CT         (3-50) 1 2 1 1 1 1 ( ) ( ) 2 r r n r n n n Y h d QT r         (3-51)
The damage stress-like variable can be transferred to the stress space by Legendre transformation. The detailed description of this transformation can be found in [START_REF] Issa | Modélisation et simulation numérique des procédés de fabrication sous conditions extrêmes[END_REF].

1 1 1 1 X R n n n n Y Y Y Y         (3-52)   1 1 1 11 11 2 2 1 1 1 0 1 3/ 2 1 ( ) ( ) ( : ) : 4 ( )(1 ( ) ) 2 ( )(1 ( ) ) 3 ( ) ( ) ( T ) ( : ) 2 ( )(
                                       (3-53) 1 1 1 2 1 3 ( ) : 4 ( )(1 ( ) ) X n n n n h Y X X C T h d         (3-54) 1 2 1 1 1 2 1 ( ) 2 ( )(1 ( ) ) r r R r n n n n n h d Y R Q T h d             (3-55)
The related strain-like variables are expressed as below:

1 1 vp vp p n n n n        (3-56) with   1 0( 1) 1 1 1 1 1 1 1 1 : 1 ( ) (1 ( ) ) ( ) / 1 ( ) ( ) r n n p D A n n p n n l n n n y n X S n n I h d h d X T R h d T                                         (3-57) Chapter 3: Numerical aspects 73       1 1 ( ) 1 1 Forward Euler 1 ( ) 1 A s y m p t o t i c ( ) n n x n n n n x a a T n n n n a T n e e a T                             (3-58) with     0( 1) 0( 1) ( : ) 2(1 ( ) ) ( ) / 1 ( ) ( ) 1 1 ( ) ( : )( : ) (1 ( ) ) ( ) / 1 ( ) ( ) ( : ) r r n n A A n n p n l n n y x n A n n n n n n p n l n n n y n n S S n n h d X T R h d T n h d S X X n S h d X T R h d T X X                                                      (3-59)       1 1 ( ) ( ) 1 1 1 Forward Euler 1 ( ) 1 Asymptotic ( ) 1 n n r n n n n r b T b T n n n n n r b T r n r e e b T d                             (3-60) with   3 2 1 0( 1) ( 1) 1 1 2 1 1 1 1 1 ( : )( : ) 1 1 1 ( ) 2(1 ( ) ) ( ) / 1 ( ) ( ) r r A n n n n r n p n n l n n n y n S S n X n h d h d X T R h d T                                (3-61)   2/3
( 1) ( 1) 3( 1) ( 1)

1 3 3 ( 1) ( ): ( ) : ( ) 3 T n n n H n J n n A n H J n w Z H T Z Z Z H T n Z w Z                 (3-62)
And finally the damage equation is given by: ( )

1 0 1 ( ) 1 ( ) ( ) (1 ) n n s T n n n n T n n Y Y T d d S T d           (3-63)

Elastic prediction

The elastic strain at time 1 n t  is given by:

                                                     (3-64)
where

( ) eH H vpH vp tr            and eD D vpD      
In this trial step, we suppose that the total strain increment   is purely elastic under the constant temperature

1 n n T T  
over the current time step, which means no dissipative phenomena(i.e. 0

   ),

Chapter 3: Numerical aspects 74 leading to:

( 0; 0; 0; 0 vp r d           )  ( 1 1 1 1 ; ; ; vp vp n n n n n n n n r r d d             ) (3-65)
The elastic trial strain at time 1 n t  is given by:

* * * 1 1 1 e n eH eD eH eD e H D n n n n n                         
* 1 2(1 ( ) ) ( )( ) 2(1 ( ) ) ( ) eD D D n ne n n n ne n S h d T S h d T                (3-67)   * 1 1 (1 ( ) ) ( ) tr( ) (1 ( ) ) ( ) tr( ) 3 1 ( ) ( ) e e e e H n n en n H n n en n en n n h d T I h d T I h d T T T I                             (3-68)
then the trial stress * 1 n   can be expressed by:

  * * 1 1 * * * 1 1 1 1 2(1 ( ) ) ( ) (1 ( ) ) ( ) 3 1 ( ) ( ) 2(1 ( ) ) ( ) ( 1
h d T h d T h d T T T I S h d T h d T I h d                                                         1 ( ) e n e n n n T T T I      (3-69) * * 0( 1) 0( 1) * * ( 1) 1 ( ) 1 * 0( 1) ( ) 2 : 2(1 ) ( )( / 1 ( )) : 2(1 ) ( )( / 1 ( )) n n n n c n d n n T c n l n n n y n n n n T n l n n y n S S S S X d X T R d T X X S d X T R d T                   (3-70) * * * * * 1 0( 1) 1 ( 1) ( 1)
: : ,,, )

n n n n n x n x n n n S X S S S and S X X X         (3-71)     c* c* c* 1 ( 1) 1 ( 1) 1 c* c* c* 1 ( 1 ) 1 3 ( ) : ( ) : ( ) det det ( ) : ( ) n d n n i n d n n H n i n d n n J Z S X H T S X Z Z H T S X                      (3-72) Noting that * * 1 ( 1 ) c n dn n Z S X     , the trial viscoplastic "yield function" * * 1 1 ( ,
vp n n n n n n f X R d T   
is the particular case of Eq. (3-40) in which ( 1) 0

v n    .   ( ) 1/3 3 * * 1 1 3 * * 1 1 ( , , , , ) ( ) 1 1 T n n n H J vp n n n n n n n y n n n Z w Z R f X R d T T d d              (3-73)
vp n n n n n n f X R d T   
 , then the assumption of elastic behavior in this isothermal trial state is true, and the solution of the present loading step is given by:

* 1 1 1 1 1 1 1 1 1 , ( ) , ( ) , ( ) , (
)

vp vp n n n n n n n n n n n n n n n n a X X b R R r r c Y Y d d d                              
(3-74)

Inelastic correction

If * * 1 1 ( , , , , ) 0 vp n n n n n n f X R d T   
 , the state variables at trial state will be corrected at time 1 n t  by computing the best values of

1 1 1 ( 1 ) 1 , , , , n n n vn n X R d        fulfilling the viscoplastic criterion i.e. 1 1 1 1 ( 1 ) 1 ( , , , , , ) 0 vp n n n n vn n n f X R d T        
 , the corrected stress-like state variables can be obtained. The hydrostatic and deviatoric part of the stress with respect to the trial stress can be written in the following form:

* 1 1 1 1 1 ( ) 2(1 ( ) ) ( )( ) 2(1 ( ) ) ( ) 1 ( ) eD D vp p n n nen n n nen n n h d S h d T S h d T n h d                         (3-75)   1 * 1 1 1 1 (1 ( ) ) ( ) tr( ) 1 ( ) 3 1 ( ) ( ) 1
h d T I h d h d T T T I h d                                 (3-76)
The Cauchy stress tensor is then given by:

  * * 1 1 1 1 1 1 1 1 1 * * 1 1 1 1 1 1 (1 ( ) ) ( )( ( )) 2(1 ( ) ) ( )( ) 3 1 ( ) ( ) 1 ( ) 1 ( ) 2(1 ( ) ) ( )( ) 1 ( ) 1 ( )

d T n h d T T T I h d h d S h d T n h d h d h

                                                               1 1 * * * 1 1 1 1 1 1 1 ( ) ( ) 1 ( ) 1 ( ) ( 1 ) ( 1 ) 2 ( 1 () ) ( ) ( ) 1 ( ) 1 ( )

d T T T I h d h d S h d T n h d h d

                                      (3-77)
The kinematic hardening:

( ) ( ) 1 1 1 1 1 ( ) ( ) 2 1 ( ) (1 ) 1 ( ) 3 ( ) n n x a T a T n n n n n n n n h d C T X Xe h d e n h d a T                    (3-78)
The isotropic hardening:

Chapter 3: Numerical aspects 76 ( ) ( ) ( ) ( ) 1 1 1 1 ( ) 1 ( ) ( ) 1 ( ) (1 ) ( ) 1 ( ) n n n n n T b T T b T r n n n n n n T n n h d Q T R R e h d e n b T h d                       (3-79)
The viscous stress:

    1 1 vp n m T vp vp n n K T t             (3-80)
With the discretized nonlinear equations developed above. Considering the calculation efficiency, the discretization equation system can be reduced to four equations with four unknowns   ,

1 n d  , 1 n X  and 1 p n n  as given below:   ( ) 1/3 3 3 1 1 1 1 1 1 ( 1) 1 1 ( ) * 1 1 1 0 1 1 1 1 1 ( ) 1 1 ( , , , , ) ( ) 0 1 ( ) 1 ( ) ( , , , ) ( ) ( , , , , ) 0 ( )
(1 ( ) ) ( , T n n n c c H J p vp n n n n n n n y n v n n n s T p p n n n n n n n n n n n n T n n n n Z w Z R f d n X T T h d h d Y n d T Y T g d n X T d d S T h d h d                                                     1 ( 1) 0( 1) 1 1 1 1 1 1 1 1 ( ) 1 1 1 1 1 1 , , , ) : 0 1 ( ) ( )( / 1 ( ) ( )) 2 ( , , , , ) ( )(1 ( )
) ( ) 1 0 3 ( ) n n n n A n n p p d n n n n n l n n n y n x p a a T n n n n n n n n n n n n n X S n X T n I h d X T R h d T n d n X T X C T h d T e e a T                                                                              (3-81)
This system is solved thanks to the classical linearization by Newton-Raphson method according to:

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ( ) ( ) ( ) : ( ) : 0 ( ) ( ) ( ) : ( ) : 0 ( ) ( ) ( )
p k k k k k n n n n n n n n p p k k k k k n n n n n n n n p k k k n n n n n p f f f f f d n X d X n g g g g g d n X d X n h h h h d d n                                                                                  1 1 1 1 1 1 1 1 1 1 1 1 : ( ) : 0 ( ) ( ) ( ) : ( ) : 0 p k k n n n p k k k k k n n n n n n n n p h n X X d n X d X n                                                          (3-82)
or under the matrix form:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) k k k k n n n n p k n k k k k n n n n k p n k k k k k n n n n n p k n k k k n n n n p f f f f d X n f g g g g g d X n h h h h h d X n d X n                                                                            1 1 1 1 n p n n k d n X                                                                  (3-83)
This system can be solved with respect to    ,

1 n d   , 1 n X   and 1 p n n   : 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
k k k k n n n n p k k k k n n n n p n p k k k k n n n n n p n k k k n n n n p f f f f d X n g g g g d d X n h h h h n d X n X d X n                                                                                 1 1 1 1 1 k n k n k n k n k f g h                                                 (3-84)
Once the nonlinear system is solved by the Newton-Raphson iterative scheme and the iterative process converges, the state variables should be updated at time 1 n t  using Eqs. (3-71) to (3-74) for the current time increment.

The problem is solved with the material Jacobian matrix which is based on the derivatives below:

(1)

vp f    , 1 vp p n f n    , 1 vp n f X    and 1 vp n f d    ( 1) 1 1 1 ( 1) 1 1 ( 1) ( 1) 1 1 1 1 1 ( 1) ( 1) 1 1 1 1 1 1 ( ) 1 ( ) : 1 ( ) : 1 ( ) r c vp vp d n c n n n A n n n c c vp A n d n n p p p n n n n c c vp A n d n n n n S R X f n h d h d n S X f n n n h d n S f I X X h d                                                                                ( 1) ( 1) 1 3/ 2 ( ) 2(1 ( ) ) 1 ( ) ( ) 1 2(1 ( ) ) 1 ( ) r r r c c vp eq A n d n n n n n n n r n n n n n h n S X f d d d h d h d R h d R d h d h d                                                                              (3-85) Chapter 3: Numerical aspects 78 (2) g    , 1 p n g n    , 1 p n g X    and 1 n g d      1 1 0 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 (1 ( ) ) ( 1 ( ) )
(1 ( ) )

(1 ( ) )

1 s s n n n n n n n n s n n p p n n n n n s n n p n n n n n n n Y Y Y Y Y g s S S h d S h d Y Y Y g s S S h d n n Y Y Y g s S S h d X X s T g d                                                                                                 1 1 0 1 1 1 1 0 1 1 1 ( ) ( ) 1 ( ) n n n n s T n n n T n n n n s T n n T n n Y Y T Y S T d S T h d Y Y T h S T h d                                          (3-86) (3) h    , 1 p n h n    , 1 n h X    and 1 n h d    1 0( 1) A( 1) 1 1 1 1 1 1 A( 1) ( 1) ( 1) 1 ( 1) 1 0( 1) 1 1 1 : 1 ( ) ( )( / 1 ( ) ( )) ( ) ( ) :
:

: : ( ) 1 1 ( ) p n n n d p n l n n n y n p p p n n n d n d n p n d n n n d d p n n l X S n h I h d X T R h d T n S X S X S X X S h I I n h d X                                                                  A( 1) ( 1) 1 1 1 1 1 1 1 ( 1) ( 1) ( 1) 1 ( 1) 1 1 1 0( 1) 1 1 1 : : ( )( / 1 ( ) ( )) ( ) ) (
n n I n n T R h d T n S X S X I I S X n n X S h I X h d X                                                            ( 1) 1 1 1 1 1 ( 1) ( 1) ( 1) 1 1 1 1 ( 1 ) 1 0( 1) 1 1 1 1 : ( )( / 1 ( ) ( )) ( ) ) ( ) ) : : : : : ( ) ) 1 1 ( ) ( )( / 1 ( p A n p n n n n y n p p p p n n A n d n d n A p p n n n n d n n n d p n n l n n n X T R h d T n S X S X n X X n S X X S h I d h d X T R h                                                             ( 1) ( 1) 1/ 2 1 1 1 1 ( 1) ( 1) ( 1) ( 1) 1/ 2 1/ 2 1 1 1 1 ( ) : 2(1 ( ) ) ) ( ) ) ( ) ( ) ( ) : : : 2(1 ( ) ) 2(1 ( ) ) p p A n A n n n n y n p p p p n A n A n A n d n n n n n h n n d h d d T h n n h n S X d d h d h d                                                                                                        (3-87) (4)     , 1 p n n     , 1 p n X     and 1 n d     Chapter 3: Numerical aspects 79       ( ) 1 1 1 ( ) 1 1 1 1 ( ) 1 1 1 1 2 ( )(1 ( ) ) ( ) 3 2 ( )(1 ( ) ) 1 3 () 2 1 ( )(1 ( ) ) 1 3 ()
I C T h d e a T X X                                                                          ( ) 1 1 1 1 2 ( ) ( ) ( ) 1 3 ( ) n n n x a a T n n n n n n n n n X n C T h T e e d d a T                                               (3-88)
The partial derivatives of state variables with respect to   are given by:

1 ( 1) 1 0( 1) 1 2 ( ) ( ) c n d n n n c l n l n A A S S X S X T X T            and 1 ( 1) 1 1 ( ) p n d n n p l n A S S X X T        (3-89)   0( 1) 0( 1) 1 1 : 2 / 1 ( ) ( ) r n n n n yn S S A R h d T            and   ( 1) ( 1) 1 1 : 2 / 1 ( ) ( ) r n n n n yn X X A R h d T            (3-90) 1 1 ( 1) 1 1 1 0 0 1 1 2 2 0 1 2 ( ) ( ) ( ) ( ) ( ) ( )
n c n d n n n n c c c c l n l n l n l n c c n l l S X X T X S A A S A A A A S X S S X X T X T X T X T                                                (3-91) 1 1 ( 1) 1 1 1 1 1 1 ( ) ( ) ( ) ( ) n p n d n n n n p c l n l n p n l X X T S A S A A S X X X T X T                                 (3-92) 0( 1) 0( 1) 0( 1) 1 1 1 1 1 1 1 : 1 : 2( ) ( ) ( ) 1 ( ) r n n n n n n n y n y n y n S S S A A S R X A R R R h d                                                           (3-93) 1 1 1 1 1 1 1 1 1 : 1 : 2( ) ( ) ( ) 1 ( ) r n n n n n n y n y n y n A A R X X X X R R R h d                                               (3-94) 0( 1) 1 1 n n n S S X A                (3-95)   1 1 1 1 ( ) 2 ( ) p n n e n n S h d T n            ( 3 -9 6 ) ( ) ( ) ( ) 1 1 1 2 ( )(1 ( ) ) ( ( )) 3 1 ( ) n n x a T n a T n n n n n n n e X C T h d e a T h d                             (3-97) Chapter 3: Numerical aspects 80 ( ) ( ) ( ) 1 1 1 ( )(1 ( ) ) ( ( )) 1 ( ) n n r r b T r n b T n n n n n n n e R Q T h d r e b T h d                              (3-98) 1 1 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1 1 1 ( ) : 2 ( )(1 ( ) ) 3 ( ) : : 4 ( )(1 ( ) ) ( ) ( )(1 ( ) r X R n n n n n n n e n n X n n n n n n n R n r n n n Y Y Y Y Y S h S T h d Y X X h X X C T h d Y h d R Q T h                                                                                    1 2 1 ) r n n R d                    (3-99)
The partial derivatives of state variables with respect to the 1 n d  are given by:

( 1)

1 1 c d n n n n n n n Z c c n n n n n n n l l S A A S X S X X A S A A d d d d d d d X X                                                           (3-100) ( 1) 1 1 1 1 1 1 1 1 1 p d n n n n c n n n n l S A S X X A d d d d X                            (3-101)     1 0( 1) 0( 1) 0( 1) 0( 1) 1 1 1 2 1 1 1 1 1 1 1 : ( ) 1 : 2(1 ( ) ) 1 ( ) 2 r r r n n n n n n r n n n n n n y n n y S S S S A R R h d d d d h d R h d R                                                 (3-102)     1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 : ( ) 1 : 2(1 ( ) ) 1 ( ) 2 r r r n n n n n n rn n n n n n y n n y A X X X X R R h d d d d h d R h d R                                                 (3-103)       1 * 1 1 1 1 0 1 1 ( ) ( ) : 3 ( ) 1 2 1 ( ) e e p n e n n n n e n n n S h d h T n T T T d h d                            (3-104)                1 ( ) 1 1 1 ( ) 3/ 2 1 2 ( ) 1 3 1 ( ) (1 ( ) ) ( )(1 ) 3 1 ( ) n n n a T a T n n n x n n n n a T n n x n n n h C T e X e n d a T h d h d C T h e n a T h d                                                   (3-105)              1 1 1 1 ( ) 3/ 2 1 1 1 1 ( ) ( 1 ) ( ) (1 ) 1 ( ) 2 1 ( ) n n n r r r b T b T b T r n r n n n n n n n n n n R e h e h Q T d e r d n d b T h d b T h d                                                        (3-106) 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 0 ( ) ( ) ( : ) : 2 ( )(1 ( ) ) ( )(1 ( ) ) 3 ( ) ( ) ( T ) 2 ( )(1 ( ) e e e X
h d T T T h d                                                                          1 3/ 2 1 1 1 1 1 1 1 2 3 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 ( : ) ) 3 3 ( ) : : 2 ( )(1 ( ) ) 8 ( )(1 ( ) ) ( ) ( )(1 ) ( )(1 ( ) ) e r r r r H n n n X n n n n n n n n n n n R n rn n rn n n n n n n n I d Y X h X X X d d C T h d C T h d Y d R R h d R d d Q T d Q T h d                                                                 1 1 2 2 ( 1 ) 2 1 1 1 1 3 1 ( ) ( 1) (1 ( ) ) ( ) 2 ( )(1 ( ) ) r r r r n n r n r n n n r n n R d h R d h d d h Q T h d                                                (3-107)
The partial derivatives of state variables with respect to the kinematic hardening are given by:

1 1 ( 1) 1 1 1 1 0 0 1 1 1 1 1 1 2 2 0 1 2 ( ) ( ) ( ) ( ) ( ) ( ) n c n d n n n d n n c c c c n n n n l n l n l n l n c c n l l S X X T X S A A S X X A A A A S X S S I X X X X X T X T X T X T                                              (3-108) 1 1 ( 1) 1 1 1 1 1 1 1 1 1 ( ) ( ) ( ) ( ) n p n d n d n n p c n n n n l n l n p n l X X T S A S A A S X I X X X X X T X T                              (3-109) 0( 1) 0( 1) 0( 1) 0( 1) 1 1 1 1 1 : : 2( ) ( ) n n n n p n n n n y n y S S S S A X X n R R                                 (3-110) 1 1 1 1 1 1 1 : : 2( ) ( ) d n n n n n n y n y A X X X I X X R R                             (3-111) 0( 1) 1 1 1 n d n n n S S A I X X            (3-112)   1 1 1 1 1 1 ( ) 2 ( ) p n n n e n n n S n h d T X X                 (3-113) 1 1 2 1 1 3 : 2 ( )(1 ( ) ) d n n n n n Y X I X C T h d          (3-114)
Finally, the partial derivatives of state variables with respect to 1 p n n  are given by: Chapter 3: Numerical aspects

82 1 1 ( 1) 1 1 1 1 1 0 0 1 1 1 1 1 1 1 2 2 0 1 2 ( ) ( ) ( ) ( ) ( ) ( ) n c n d n p p n n n n n p p p p p c c c c n n n n n l n l n l n l n c c n l l S X X T X S A A S n n A A A A S X S S X n n n n n X T X T X T X T                                                  (3-115) 1 1 ( 1) 1 1 1 1 1 1 1 1 1 1 1 ( ) ( ) ( ) ( ) n p n d n n n n p p p p p p c n n n n n l n l n p n l X X T S A S A A S X X n n n n n X T X T                                  (3-116) 0( 1) 0( 1) 0( 1) 1 1 1 1 1 1 1 1 : : 2( ) ( ) n n n n n p p p p n n n n n y n y S S S A S X A n n n n R R                                             (3-117) 1 1 1 1 1 1 1 1 1 : : 2( ) ( ) n n n n p p p n n n n y n y A X X X X n n n R R                                 (3-118) 0( 1) 1 1 1 1 1 n n n p p p n n n S S X A n n n                (3-119)   1 1 1 1 ( ) 2 ( ) d n n e n p n S h d T I n             (3-120)   ( ) 1 1 1 1 ( ) 2 (1 ( ) ) 1 3 ( ) n x a T n n n p p n n n C T X n h d e a T n n                (3-121) 1 1 1 1 1 2 2 1 1 1 1 1 ( ) 3 : : 2 ( )(1 ( ) ) 2 ( )(1 ( ) ) n n n n n p p p e n n n n n n n Y S X h S X T h d C T h d n n n                                
(3-122)

Calculate the heat generation

In order to calculate the internal flux vector   int G which is described in the resolution scheme, we need to calculate the internal heat generation. The expression of pl R  is given below:

1 1 1 1 1 1 1 1 : : : : : : p pl n n n n e n n n n n n n n n R X R r Y d X R Y T r d T T T T                                              (3-123)
Where:

  * * 1 1 1 1 1 1 1 1 1 1 : [(1 ( ) ) ( )( ( )) 2(1 ( ) ) ( )( ) 3 1 ( ) ( ) 1
d T T T n                                       1 1 1 1 1 1 1 2 2 : ( 1 ( ) ) ( ) : ( 1 ( ) ) ( ) ( : : ) 3 3 n n n n n n n n n n X h d C T h d C T                        2 1 1 1 1 1 1 (1 ( ) ) ( ) (1 ( ) ) ( )( ) r r n n n n n nn n n R r h d Q T r r h d Q T r r r                  ( ) 1 1 0 1 ( ) 1 ( ) ( ) (1 ( ) ) n n s T n n n n T n n Y Y Y T Y d S T h d            
The isentropic term needs to be calculated with the derivatives of the state variables

1 1 1 1 ( , , , ) n n n n X R Y     
with respect to temperature. This process concerns the calculation the derivatives of the temperature dependent material parameters, see the work of [START_REF] Lestriez | Modélisation numérique du couplage thermo-mécaniqueendommagement en transformations finies[END_REF][START_REF] Issa | Modélisation et simulation numérique des procédés de fabrication sous conditions extrêmes[END_REF] for more details.

Practical implementation of user defined subroutine-VUMAT

The user defined subroutine provides us an approach to perform simulations using proposed constitutive models which are not existed in commercial finite element (FE) code to accurately represent the material behavior. The user defined material model can be implemented in user subroutine UMAT and VUMAT in ABAQUS/Standard and ABAQUS/Explicit respectively. In this thesis, the VUMAT is chosen to make the simulation in ABAQUS/Explicit, since it is relatively easy to reach convergence with a reasonable computational time. The significant features of the VUMAT as well as the differences between the UMAT and VUMAT are listed below:

 The stresses and strains are stored as vectors. For plane stress elements, plane strain and axisymmetric elements, the storage scheme is the same in UMAT and VUMAT. The difference is in storage of three-dimensional elements, as shown below:

VUMAT:    In VUMAT, the initial values at the beginning of each increment are put in the 'old' arrays, the updated values at the end of the each increment must be assigned in the 'new' arrays.

 There is no need to define the Jacobian matrix.

 The time increment cannot be redefined.

 The vectorized interface is applied. This structure allows vectorization of the routine, but branching inside loops has to be avoided.

Conclusions

In this chapter, the numerical aspects of the fully coupled thermo-mechanical model were discussed. The following aspects were presented: the strong and weak forms of the initial and boundary value problem, the time and space discretization of the IBVP, the global resolution scheme of the dynamic explicit solver and the local constitutive integration scheme, the implementation of user defined subroutine-VUMAT.

In the next chapter, we will show the high predictive capabilities of using the constitutive equations with FEM to account for the complex phenomena in hot sheet metal forming.

Introduction

In order to show and analyze the capabilities as well as the limitations of the proposed model(summarized in Table 4.1), in this chapter an exhaustive parametric study was conducted to explore its newly developed features. The parametric study of the local response was done with Representative Volume Element (RVE)

(which in fact is nothing but a material point representing an integration point), to analyze the effects of anisotropy and tension-compression asymmetry, the distortion of the yield surface, the temperature and strain rate, the hardening asymmetry, the microcracks closure, the stress triaxiality and Lode angle.

Table 4.1 Summary of the constitutive equations

 

1/3 3 3 ( ) 0 1 ( ) 1 ( ) r c c H J y Z w Z R f T h d h d                3 ( ): ( ):( ) ;( , ) det det ( ) : ( ) i i i d i d H i i i i d J Z S X H T S X i c p Z Z H T S X               (a)     0 0 0 1 2 : : 2(1 ( ) ) ( ) / 1 ( ) ( ) 2(1 ( ) ) ( ) / 1 ( ) ( ) r r c d c l y l y S S X X S S X S h d X T R h d T h d X T R h d T                  (b)   0 0 1 : 2(1 ( ) ) ( ) / 1 ( ) ( ) r p d p l y S S S S X h d X T R h d T         
with 0 : :

x x S X S S S and S X X X     (c)   0 (1 ( ) ) ( ) 3 1 ( ) ( ) 1 2(1 ( ) ) ( ) H e e

S eH eD e e e e h d T h d T T T h d T

                            (d) 2 (1 ( ) ) ( ) 3 X h dC T          (1 ( ) ) ( ) r R h d QTr r         e r Y Y Y Y d          (e) vp p f D n    ( ( , ) ) x f n a p           ( , ) r f r n b pr       ( ) 0 ( ) ( ) ( , ) (1 ( ) ) s T f T Y Y T d S T h d                  (f) 1 ( ) 1 ( ) ( ) Tanh( ) 2 2 c c h h T h T h        (g) 1 1 2 2 3 2 1 ( , ) ( ) [1 ( )][ ( ) ( ( ) ( )) ( ( ( )))] 2 L L s a p a T Tanh a T a T a T Tanh p p T           (h) 1 1 2 2 3 2 1 ( , ) ( ) [1 ( )][ ( ) ( ( ) ( )) ( ( ( )))] 2 L L s b p b T Tanh b T b T b T Tanh p p T           (i)   ( , ) ( ) ( ) ( ) Tanh( ) s L s h t e n s h L S T S T S T S T       with 3 3/ 2 2 3 3 6 2 1 1 arccos( ) 2 L L J J         (j)

Initial anisotropy and asymmetry

Recall that in the case of plane stress, the yield function in Table 4.1 (Eq. (a)) contains 9 parameters which determine the shape of the yield surface. The case of initial plastic isotropy (or von Mises yield function) corresponds to the generic values given by: F=G=H=F'=G'=H'=0.5 and N=N'=1.5. We will focus on the study of parameters F, G, H, F', G', H' and w effects, by varying one by one every parameter while the others have fixed values. Note that, as shown in Figure 4-1, the case of w=0 reduces the yield surface to the classical von Mises yield stress, while when considering w=2 or w=-2, the yield asymmetry is activated to get respectively larger yield stresses in tension (w=2) and in compression (w=-2). F, G and H (involved in J 2 part of equivalent stress) are more significant than those of parameters F', G' and H' (involved in J 3 part of equivalent stress). It is also observed that the effect of these anisotropic parameters are highly coupled and requires an appropriate optimization method to be accurately determined.

Distortion of the yield surfaces

In order to examine the distortional parameters effects, new simulations are made with the proposed model using the material parameters given in Table 4.2, the initial plastic isotropy is considered with w=0 to avoid the initial plastic flow anisotropy effect on the yield surface.

Table 4.2 Assumed material parameters 2 0 0 0 1 0 0 0 1 0 0 0 0 0 0

E(GPa) υ σ y (MPa) C(MPa) a Q(MPa) b F=G=H L=M=N F'=G'=H' L'=M'=N
1 0 0 1 1 1 1 1 1 0 1 0 , 1 0 0 , 0 0 0 , 0 0 1 , 0 1 0 , 0 1 0 6 2 2 2 2 3 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 e e e e e e                                                                     
From Figure 4.3, it can be found that, when the effect of 2 c l X is neglected, the distortional parameter 1 c l X affects the distortional ratio of the yield surface, but the size of the yield surface remains the same. Smaller value of 1 c l X leads to a high distortion of the yield surface.

Combined SD effect and distortion of yield surface with thermal effects

The influence of temperature on the initial yield surface is treated in this section, the initial anisotropy and tension-compression asymmetry are affected by the temperature through changing the values of anisotropic parameters at different temperatures. We use a real case of Mg alloy AZ31 to show this phenomenon by using the identified parameters in Table 5.5 of Chapter 5. From Figure 4.4b, the tension-compression asymmetry is evolving with the temperature. At room temperature, the yield surface shows significant tension-compression asymmetry of the yield stresses, while at high temperature, this phenomenon becomes very weak. The comparison of yield surfaces obtained from these two models shows the interest of the proposed model to capture the evolving tension-compression asymmetry for different temperatures.

The temperature effect on the induced anisotropy is included by using temperature dependent distortional parameters. The assumed material parameters are given in Table 4.3. In order to study the combined distortion hardening effect and tension compression asymmetry effect under complex loading paths, a parametric study based on the identified material parameters [START_REF] Zhang | Thermomechanical modeling of distortional hardening fully coupled with ductile damage under non-proportional loading paths[END_REF] of Mg alloys, given in Table 4.4, was conducted. During this loading path, the changes of size, location and shape of the yield surfaces are revealed. In these evolving surfaces on different directions, the yield surfaces show a high curvature in loading direction for all the applied loading paths. The evolution of the yield surfaces at these two loading paths have the same trend with the numerical predictions of other distortional hardening models which was reported in [START_REF] Shi | On the macroscopic description of yield surface evolution by means of distortional hardening models: application to magnesium[END_REF].

Table 4.3 Assumed material parameters

No. T σ y (MPa) C(MPa) a Q(MPa) b 1 c l X 1 c l X 2 c l X T1
By comparing the yield surfaces at two temperatures, the mentioned effects at different loading paths are more pronounced at temperature of 65°C compared to that at 150°C. This is caused by the decreasing of isotropic, kinematic and distortional hardening with the increasing of temperature i.e. temperature-induced softening.

Effect of the asymmetric hardening parameters

The hardening asymmetry is introduced through the kinematic and isotropic hardening parameters a and b by defining dependence functions to both normalized Lode angle and equivalent plastic strain. The value of a and b are reduced to a1 and b1 respectively under tension (

1 L   )
, so a1 has the same effect as the classical model in which a higher value of a1 leads to a quicker saturation of kinematic hardening. The same remark can be made for isotropic hardening for which a higher value of b 1 leads to a faster saturation (see [START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. The effect of the parameters 1  , 2  , a 2 , a 3 , b 2 , b 3 and p s on the hardening asymmetry without the damage effect is investigated using assumed values of the material parameters given in Table 4.5. First, we study the influence of different choices of 1  for the stress-strain evolution, as shown in Figure 4.8.

This parameter only affects the shear response, without changing the response in tension and compression.

Indeed, from Table 4.1(h) the parameter 1  is designed to adjust the effect of L  only in shear. When the value of 1  decreases, the shear stress is decreased at the beginning and increased at high plastic strain level. The effect of varying parameter a2 on the stress-strain curves in tension and compression is now investigated while keeping the kinematic hardening parameter a always positive, to ensure the thermodynamic admissibility of the model which requires that the ratio C/a (as well as the ratio Q/b for isotropic hardening) should be positive (see Chapter 2). Figure 4.9(a) indicates that a2: (i) has no effect in tension curves, (ii)

while it has an important effect in compression. Indeed, the smaller values of a 2 lead to increase the stress level in compression from the early stage of hardening. This effect is much depreciated for high values of plastic strain (p>0.2). As shown in Figure 4.9(b), the value of a2 affects deeply the value of parameter a at the beginning of hardening (lower values of plastic strain), which decreases significantly from its initial value going to its asymptotic value which is a = a 1 -a 3 = 50.0 -35.0 = 15.0 (see Table 4.4 and Table 4.1 (h)). Accordingly, when the value of a is high, the kinematic hardening saturates very fast giving the nonlinear (compression) hardening curves at low plastic strain (p<0.05). As the value of a decreases approaching its asymptotic value, the kinematic hardening nonlinearity decreases leading to a quasi-linear kinematic hardening and the stress increases linearly again. Finally, when a = a1 -a3= 15.0 becomes constant a slightly nonlinear kinematic hardening is observed again and the kinematic hardening tends to slightly saturate (Figure 4.9a). It is worth noting that a is still constant in tension with a = a 1 =50.0 (see Table 4.4

and Table 4.1 (h)). The impact of this evolution of the kinematic hardening parameter with respect to the parameter a2 for the shear loading is shown in Figure 4.9c. Similarly, the larger value of a2 results in a decrease of the shear stress for plastic strain lower than 0.1. in tension, a strong effect is observed in compression and a much less effect is noticed in shear. Indeed, from these figures we can see that parameter a 3 has no effect on the first nonlinear hardening in compression (for p<0.05), but has a significant effect on the second hardening stage (for p>0.05). Note that when considering a 3 = a 1 = 50, the hardening will have linear evolution (no saturation will occur since a = 0). The increase of a 3 leads to the increase in shear stress at high plastic strain, but the effect is very weak. The isotropic hardening parameter b3 has the same influence with the kinematic hardening parameter a3 on the stress strain evolution for tension and shear, but the effect is delayed for large plastic strain, as illustrated in Figure 4.12. Through this parametric study, the different sigmoidal hardening of Mg alloys in compression can be achieved by adjusting the parameters of the proposed model. We can finally note that a 2 controls the first hardening saturation stage and a3 controls stress saturation at high plastic strain in compression. For shear loading path, a 2 affects the stress evolution at low plastic strain and a 3 has a weak effect on stress curve at high plastic strain level. The influence of isotropic parameters b 2 and b 3 are similar with kinematic hardening parameters respectively a2 and a3. However, depending on modules C and Q the effect is important form the start of the plastic flow (as the case of the kinematic hardening) or for large plastic strain (as the case of isotropic hardening). The parameter 1  has influence only under shear loading paths without effect on uniaxial tension and compression. The parameter 2  has influences on the curvature of the stress-strain curve concerning loading cases defined by normalized Lode angle 1 L   (i.e. different from uniaxial tension). The parameter p s defines the size of the initial low strain plateau in compression, and higher value of p s leads to stress decrease in shear.

Hardening asymmetry fully coupled with ductile damage

We focus in this section on the effect of the full coupling of the ductile damage with the strength-differential (SD) effect and hardening asymmetry (HA) effect. The generic set of damage parameters are given in Table 4.6. If only the SD effect is taken into account, the stress-strain curves for tension and compression loadings will have the same trend in evolution but at different stress levels due to the difference in terms of initial yield stresses between these two loading paths. So the final fracture is delayed for the loading having the lower yield stress which is the case of the compression as shown in In order to obtain a more predictive ability related to the physical phenomena, the microcracks closure effect is considered to enhance the damage coupling for loading cases altering tension and compression stress states. As described in detail in chapter 2, damage evolution shall be different in tension and compression due to the microcracks closure effect described by introducing parameter h, as given in Eq.(2-60). When this parameter varies in the range from 0 (no damage effect in compression) to 1 (damage effect in compression is similar to tension) the response in compression will be varied without any change in tension response, as illustrated in Through this parametric study we can highlight, for compression loading, competition ruling damage evolution between hardening (due to HA effect allowing increasing of damage evolution) and microcracks closure (allowing reduction of damage evolution).

Effects of temperature and strain rate fully coupled with ductile damage

In this section, parametric study on the effects of temperature and strain rate (viscoplasticity) is performed in a material point to examine the response of the proposed anisotropic thermo-elasto-visco-plastic anisotropic model fully coupled with isotropic ductile damage. For each of the temperature dependent

parameters                 , , , , , v y

P T E T C T Q T S T K T T

  , the dependence to the absolute temperature is postulated as [START_REF] Lestriez | Modélisation numérique du couplage thermo-mécaniqueendommagement en transformations finies[END_REF]:

  0 0 0 1 [( ) / ( )] f P T P T T T T          (4-1)
where 0 P is its value at the reference temperature 0 T , f T is the melting temperature of the material and  is a temperature independent material parameter.

The simulations are made with the proposed model using the material parameters given in Table 4.7, assuming the initial plastic isotropy with w=0; F=G=H=0.5 and L=M=N=1.5; Here f T =650 °C; 0 T =20 °C the parameters  =4 for Young's modulus;  =1.08 for capital S of damage equation, and  =1.03 for the remaining parameters.

Table 4.7 Assumed material parameters To study the damage effect on the complex loading path with temperature effects. A non-proportional loading path shown in figure 4.18 is applied to RVE at three temperatures and constant strain rate of 10 -4 s -1 using the material parameters reported in [START_REF] Zhang | Thermomechanical modeling of distortional hardening fully coupled with ductile damage under non-proportional loading paths[END_REF]. Note that all the simulations performed under this complex loading path are conducted until the final fracture of the RVE. With the temperature increase, both of these two stresses are decreased (temperature-induced softening) at all the stages of this loading path. The damage is obviously delayed as the temperature is high. In order to compare the damage and distortional hardening effects, a case uncoupled with damage and a case without distortion are computed as shown in Figure 4.19b. The damage effect results in stress decreasing and it becomes more obvious at higher plastic strain level. Accordingly, at the final stage of the loading path III, the stress vs strain curve evolves to fracture very quickly. The distortional hardening has a very significant impact on the stress strain evolution after comparing with the case without distortion. The notable changes can be observed when the loading direction is suddenly changed. -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 0,0 0,1 0,2 0,3 0,4 0,5

E(GPa) υ σy(MPa) C(MPa) a Q(MPa) b K v m v S  s γ β Y 0
h c =0.2 h c =0.1 h c =0.3 h c =0.4 h c =0.5 h c =0.6 h c =0.7 h c =0.8 h c =0.9 Effect of h c  eq  h c =1.0 -0,2 0,0 0,2 0,4 0,6 0,8 0,0 0,1 0,2 0,3 0,4 0,5 h c =0.2 h c =0.1 h c =0.3 h c =0.4 h c =0.5 h c =0.6 h c =0.7 h c =0.8 h c =0.9 Effect of h c  eq  h c =1.0
Chapter 4: Parametric study of the proposed model 105 In Figure 4.24a is displayed the effect of parameter hc on the curves of accumulated plastic strain vs triaxiality (for critical damage value d=0.9, near the final fracture), when the condition S ten =S sh is considered, for which no Lode angle effect is considered for the damage evolution. Note that the case when hc=1.0 corresponds to an h parameter independent of triaxiality (h=1.0). With the decrease of hc , the fracture strain at low triaxiality is increased. When h c reaches a small value (i.e. h c =0.1), the difference between final fracture strain at low and high triaxiality is significantly important. Since the negative value and positive value of triaxiality represent the compressive and tensile loading respectively, the parameter hc controls the different damage evolutions from tension to compression.

The Lode angle effect with Sten=0.5 and Ssh=0.1 is shown in Figure 4.24b, the equivalent plastic strain at fracture is no longer a monotonic evolution with stress triaxiality, it increases with the increase of stress triaxiality within the range ( 0 0.33

  
) and ( 0.566 0.66

  
), meanwhile, it decreases with the increase of stress triaxiality within the remaining ranges ( The difference between the fracture plastic strain increases when the difference between S ten and S sh is -0,2 0,0 0,2 0,4 0,6 0,8 0,0 0,1 0,2 0,3 0,4 0,5 significant. When Sten=Ssh, Lode angle effect disappears from the accumulated plastic strain at fracture. With the decrease of the S sh , the Lode angle effect is more obvious.

Recall that the parameter h

 controls the shape of the evolution curve of parameter h with respect to stress triaxiality. In Figure 4.26 are displayed the curves in terms of fracture plastic strain with respect to stress triaxiality for different values of parameter h  . Note that a small value of parameter h  (i.e. h  = 0.5) gives quasi-linear evolution with insignificant difference between fracture plastic strain for negative and positive stress triaxiality. However, the difference becomes more significant when h  is increased. For example, for the case h  =6, the maximum value of fracture plastic strain reaches 0.42 for =-0.33, while the fracture strain is only about 0.12 for the stress triaxiality higher than 0.33   . -0,4 -0,2 0,0 0,2 0,4 0,6 0,8 0,0 0,1 0,2 0,3 0,4

 h =6.0  h =5.0  h =5.5  h =4.5  h =3.5  h =3.0  h =2.5  h =2.0  h =1.5  h =1.0 Effect of  h  eq Triaxiality  h =0.5
-0,2 0,0 0,2 0,4 0,6 0,8 0,0 0,1 0,2 0,3 0,4

 S =2.0 Effect of  S  eq Triaxiality  S =6.0
The effect of parameter S on the evolution of accumulated plastic strain at fracture vs stress triaxiality considering d=0.9, S ten =0.5, S sh =0.1 and h c =0.1 is shown in Figure 4.27. The parameter  S affects the shape of the curve giving, for high values of this parameter, singularities for The iso-damage surfaces including the Lode angle effect are plotted in Figure 4.28. The circular surface can be obtained when S ten =S sh (no Lode angle effect). When S sh <S ten , the shape of fracture locus at constant triaxiality varies from a right hexagon to a six-point star. When S sh >S ten , the fracture locus at constant stress triaxiality is transferred to a shape of flower with six petals. Note that for this case of isotropic plasticity, the iso-damage surfaces has a / 3  symmetry so that fracture locus at six vertices with

1 L    are the same.
It should be bear in mind that the ductile fracture tends to occur at the region between the two vertices (shear or plane strain mode), which is proved experimentally by [START_REF] Lou | Prediction of ductile fracture for advanced high strength steel with a new criterion: Experiments and simulation[END_REF].

The effect of parameter  S on the isotropic damage surface displayed in the strain deviatoric plane (considering S ten =0.5, S sh =0.1 and h c =1.0) is given in Figure 4.29. When  S = 2, the shape of fracture locus looks like a six-point star with maximum fracture strain at the axis of - The effect of stress triaxiality on the iso-damage surface without Lode angle effect (S ten = S sh =0.5) are shown in Figure 4.30. The fracture plastic strain are the same at different normalized Lode angles for a given stress triaxiality. The size of the iso-damage surface becomes smaller with the increase of triaxiality, which makes the fracture plastic strain larger for negative value of stress triaxiality than for positive stress triaxiality. The combined effect of stress triaxiality and Lode angle are clearly presented. The Lode angle and triaxiality dependence of the damage surfaces is included in our proposed fully coupled CDM model, which could make more accurate predictions of damage in sheet metal forming simulations.

)  p la n e  1    0   1   3 ' e  2 ' e  1 ' e  0,3 -0,2 -0,1 0,0 0,1 0,2 0,3 -0,3 -0,2 -0,1 0,0 0,1 0,2 0,3 (
9  s  4  s  3  s  2  s  ( )  p la n e  1    0   1   3 ' e  2 ' e  1 ' e   D II  D I Effect of parameter  s -1,0 -0,5 0,0 0,5 1,0 -1,0 -0,5 0,0 0,5 1,0 0.33    0   0.26    0.29    0.33    ( )  p la n e  0   1    1    D II

Conclusions

In this chapter, a systematic parametric study was conducted with the proposed model applied to a typical RVE subject to various loading paths. From section 4.2 to 4.5, the isotropic ductile damage was not considered, we have mainly presented the effects of initial anisotropy and tension-compression asymmetry, hardening asymmetry parameters, distortion of the yield surface, as well as the interaction of these phenomena with temperature. Then these effects were investigated again using the fully coupled model with isotropic ductile damage in section 4.5 and 4.6. A detailed description of temperature and strain rate effects was given in section 4.7. Finally, the stress triaxiality and Lode angle effect were included in the fully -1,0 -0,5 0,0 0,5 1,0 -1,0 -0,5 From all these parametric studies in this chapter, the novelty of the proposed fully coupled model were clearly presented. In the next chapter, the methodology of the identification for material parameters will be introduced. The fully coupled damage model will be applied to the simulation of different hot sheet metal forming processes.

0,0 0,5 1,0 Effect of triaxiality  0.33   0   0.26   0.29   0.33   0   1    0   1   0   1    0   1   0   1    0   1   3 ' e  2 ' e  1 ' e 

Introduction

This chapter consists of two parts: the first part focuses on the identification procedure required to determine the values of the material parameters entering the proposed fully coupled constitutive equations for Mg and Titanium alloys, and the second part concerns the applications of hot sheet metal forming for the same materials.

According to the modular structure of the fully coupled constitutive equations, the material parameters for inelastic flow anisotropy, hardenings and damage parameters are determined separately based on specific identification methodology using Python scripts, ABAQUS FE software and MATLAB-based minimization code. Both AZ31B and Ti-6Al-4V alloys are investigated.

The applications are made to hot sheet metal forming processes using the calibrated material parameters of Mg and Titanium alloys. The proposed constitutive model is firstly validated using simple tensile tests (ST), then three point bending test (TPB), circular cup deep drawing test (CCD) and cross-shaped cup deep drawing test (CSD) are performed. Finally, some comparisons between simulation results and experimental results are performed, analyzed and discussed.

Material parameters identification

The identification of the developed model material parameters is performed in various steps as schematized in Figure 5.1. This figure illustrates the ductile damage effect on a typical stress-strain curve from initial safe state (no damage) to the final fracture state. First of all, the elasticity modulus E as well as the initial yield stress y  are easily obtained from uniaxial tensile test using quasi-linear path O-A (Figure 5.1).

Second, the parameters of initial plastic anisotropy ( , , , , ', ', ', ' H G F N H G F N ) are identified using the yield stresses measured at the early stage of plastic strain (i.e. point A of Figure 5.1 conventionally defined at 0.02% of total strain) of uniaxial tensile and compressive tests in different orientations lying from the rolling to the transverse directions together with the balanced biaxial tests. Third, the hardening parameters  In tension:

( 1 2 3 1 2 3 1 2 , , , , , , , , , , s C
2 11 = cos T T     , 2 22 = sin T T    , 12 = cos sin T T       In compression: 2 11 = cos C C      , 2 22 = sin C C      , 12 = cos sin C C      
Substituting the above equations into the yield function defined by Eq. (2-82) in chapter 2, T   and C   can be obtained as follows:

1/3 3/ 2 1 2 / T y K wK         
(5-1)

1/3 3/ 2 1 2 / C y K wK          (5-2) where     4 2 2 4 2 2 1 cos 2 cos sin sin 2 sin cos K G H H H F N                 2 2 6 2 2 2 2 2 4 2 2 2 2 6 2 2 2 2 2 4 2 cos ( 2 4 ) sin cos sin ( 2 4 ) cos sin K G H GH GH H F HF F G FN H F H F F H G H H G F G G N                                                
The uniaxial tensile yield stresses 0

    1/ 3 3/ 2 0 / T y G H wH G G H              
(5-3)

      1/3 3/ 2 2 2 2 45 1 / 2 4 4 8 T y w F G N F G N G F F G                              
(5-4)

    1/ 3 3/ 2 90 / T y F H wH F F H              
(5-5)

    1/ 3 3/ 2 0 / C y G H wH G G H               (5-6)       1/ 3 3/ 2 2 2 2 45 1 / 2 4 4 8 C y w F G N F G N G F F G                              
(5-7)

    1/ 3 3/ 2 90 / C y F H wH F F H              
(5-8)

The

    1/ 3 3/ 2 / T b y G F wF G G F              
(5-9)

    1/ 3 3/ 2 / C b y G F wF G G F              
(5-10)

Together with all the functions Eq. (5-3) to Eq. (5-10), the set of eight parameters , , , , , , ,

G H F N G H F N    
are expressed using the algebraic system of eight equations expressed for the eight considered tests. In order to identify the anisotropic parameters, an objective function is built based on stress ratio as follows:

_ exp i_ exp 90 2 2 _sim _sim 0 ( ) ( 1) ( 1) b i C C b i i i T i T b Min Min                            (5-11)
where the superscripts 'exp' and 'sim' refer to the experimental and simulation data respectively and the index i represents the loading paths of uniaxial tension (i=T) and compression (i=C). The variables ,

      1/ 3 3/ 2 2 2 2 45(11 33) 1 / 2 4 4 8 T y w F H M F H M H F F H                                (5-12)       1/3 3/ 2 2 2 2 45(11 33) 1 / 2 4 4 8 C y w F H M F H M H F F H                               
(

      1/3 3/ 2 2 2 2 45(22 33) 1 / 2 4 4 8 T y w G H L G H L H G G H                                (5-14)       1/3 3/ 2 2 2 2 45(22 33) 1 / 2 4 4 8 T y w G H L G H L H G G H                               
(5-15)

For the case of non-AFR model, the objective function like Eq.(5-11) can be constructed using different Lankford R-values, giving the following expressions:

f f f f f r                                             
(5-16) The identification of the hardening and damage parameters are performed within an inverse approach, which is based on minimizing deviation value between numerical and experimental results [START_REF] Yue | Ductile damage prediction in sheet metal forming processes[END_REF][START_REF] Souto | Material parameter identification within an integrated methodology considering anisotropy, hardening and rupture[END_REF]. The procedure is realized using Matlab and Python subroutines in connection with Abaqus/Explicit FE software, the detailed process of this identification scheme is illustrated in Figure 5.2.

Anisotropic and hardening parameters calibration

In this work, both Mg alloys and Titanium alloys are investigated to validate the proposed model. The anisotropic and hardening parameters are calibrated using the experimental data in the following published works [Khan11, Gilles11,

. The results of identification are discussed in the following sub-sections.

Titanium alloys

Using the identification method described above, together with the experimental data in [START_REF] Gilles | Experimental characterization and elasto-plastic modeling of the quasistatic mechanical response of TA-6 V at room temperature[END_REF], the anisotropic parameters of Titanium alloy Ti-6Al-4V are obtained, as given in Table 5.1.

Mg alloys (Experimental data from UMD)

In this section, the experimental data are taken form the published works [START_REF] Khan | Mechanical response and texture evolution of AZ31 alloy at large strains for different strain rates and temperatures[END_REF][START_REF] Habib | Anisotropy, tension-compression asymmetry and texture evolution of a rare-earthcontaining magnesium alloy sheet, ZEK100, at different strain rates and temperatures: experiments and modeling[END_REF] of a research group from University of Maryland (UMD). Both Mg alloys AZ31 and ZEK100 are investigated to validate the predictive capability of the proposed model in describing the initial anisotropy and asymmetry at room temperature, since both of them exhibit strong initial anisotropy and tension-compression asymmetry in yield stress. The anisotropic parameters of AZ31 and ZEK100 are given in Table 5.2. ( 2 )/4 ( , )

i i i i b i T C        
).

The yield function correctly describes the experimental yield stresses of AZ31 with few differences except for the assumed balanced biaxial compressive yield stress (Figure . 5.3a). The predicted yield surface of ZEK100 fits well with the experimental measurements (Figure . 5.3b). The initial anisotropy and asymmetry of AZ31 and ZEK100 are well captured by the proposed model. The identified hardening parameters of AZ31 and ZEK100 at room temperature using the experimental data of UMD [Khan11, Habib17] are shown in Table 5.3 and Table 5.4 respectively. They are determined by using tension, shear and compression tests, while only tension and compression experimental results are available for ZEK100. The predicted stress-strain curves for AZ31 (see Figure 5.4a) in tension, shear and compression show a nice fit with the experimental results. For ZEK100, good fitting with experimental data can also be observed in both tension and compression as shown in Figure . 5.4b. The evolution of the compressive stress for both these two Mg alloys follows a typical sigmoidal function, the low plastic strain is governed by the twinning mechanism. While at higher strain level, it is mainly dominated by classical plastic slip mechanism, as we can see the sigmoidal hardening of AZ31 is more obvious than the case of ZEK100. The proposed model can well predict the unusual hardening evolution for both of these Mg alloys.

Mg alloy AZ31 (Experimental results from UW)

Tension and compression tests have been conducted over a wide range of temperatures lying between 23 to 250°C and strain rates varying between 0.001 to 1.0 s -1 taken from [START_REF] Tari | Experimental studies of deep drawing of AZ31B magnesium alloy sheet under various thermal conditions[END_REF][START_REF] Tari | Mechanical response of AZ31B magnesium alloy: experimental characterization and material modeling considering proportional loading at room temperature[END_REF][START_REF] Tari | Elevated temperature constitutive behavior and simulation of warm forming of AZ31B[END_REF]. By using these experimental data, the anisotropic parameters and hardening parameters at different temperatures are calibrated, and given in Table 5.5 and 5.6. Comparison of predicted yield surfaces and experimental measurements of AZ31B at different temperatures are shown in Figure 5.5. The SD effect is significant at room temperature, the predicted yield surface could fit the experimental data at room temperature with a satisfied accuracy. It should be noted that the high temperature weaken the SD effect, when the temperature is elevated (150°C, 200°C, 250°C), the predicted yield surfaces have a good agreement with the experimental data. Stress-strain curves in tension and compression at room temperature are slightly over estimated compared to experimental results at low plastic strain. This is may be due to the limited form of viscous stress which is introduced instantaneously for a given strain rate at the beginning of inelastic flow. The hardening asymmetry is prominent at room temperature and the numerical results describe correctly the sigmoidal shape evolution for compression at room temperature. With the temperature increase, the effect of hardening asymmetry is reduced significantly. The simulated stress-strain curves at elevated temperatures fit well with Chapter 5: Parameter identification and hot sheet metal forming applications 121 the experimental data for both tension and compression at constant strain rate 0.001s -1 . The viscosity parameters are identified using the tensile tests at different strain rates, the simulation results are compared with the experimental data in Figure 5.7. At room temperature, the strain rate sensitivity is weak and the viscous stress is small. The strain rate sensitivity becomes more evident at higher temperatures as shown in Figure 5.7c, and 5.7d. A large deviation can be observed in the prediction of stress at low viscoplastic stain for higher strain rate at 200°C, it is mainly due to the fact that the equation to calculate the viscous stress is directly added in the yield function, giving rise to a significant increase in the stress for low viscoplastic strain range. The stress-strain curves at higher viscoplastic strain can be well fitted by the simulation results. The evolution trends caused by the strain rate sensitivity can be reasonably captured by our proposed model. However, improvement of the proposed viscous stress can be acted by introducing for example a dependency of the viscous module on the equivalent viscoplastic strain.

Ductile damage parameters calibration

The influence of the damage parameters ), so the cracks start at the border of the specimen for Case_1.

Even though this case is not consistent with the experimental measurement for Ti-6Al-4V, it is discussed here to demonstrate the capacity of the proposed model in predicting ductile damage occurrence at different location of the specimen [START_REF] Thuillier | Influence of ductile damage on the bending behaviour of aluminium alloy thin sheets[END_REF]. Using the same methodology as for damage parameters identification of Ti-6Al-4V, the best values of the calibrated damage parameters for AZ31B at different temperatures are given in Table 5.8 The predicted force-displacement responses at different values of the imposed temperature are compared to the experimental results in Figure 5.17 The final fracture of notched AZ31 specimen are all located in the center of the notched specimens, which is similar to titanium specimens presented in the previous sections. The numerically predicted final fracture of the shear specimens are compared with the experiments, as shown in Figure 5.18. The fracture often takes place at one side (upper or lower part) of the sheared regions. At room temperature, the fracture occurs faster at the displacement of 1.25 mm, as presented in Figure 5.17. However, at high temperature, a severe rigid body rotation can be found at the shear region, and the fracture surface is inclined. The predicted crack path in shear specimen has a good consistency with the experimental observations (see Figure 5.18).

Material parameters dependent to the temperature

There are many choices for the thermal softening functions P(T), namely: Linear, Power Law and Polynomial functions, and so on. The Polynomial functions were chosen for the current work for more flexibility by varying the polynomial order. ( 5 -1 7 )

However, the evolution of some parameters with respect to temperature can be well captured by a lower order function, as the Young's module E. The detailed fitting curves and functions for the temperature dependent parameters of AZ31 are given in Figure 5.19, Figure 5.20 and Table 5.9.

Validation with simple test for Ti-6Al-4V

In this section, the proposed model is validated through simple tension tests using the calibrated material parameters in section 5.2 for Ti-6Al-4V. The uniaxial tension test using flat specimen of Ti-6Al-4V was reported in [Giglio12, Allahverdizadeh13]. Simulation using the identified material parameters with mesh size 0.1 mm in the central part of the specimen and C3D8R element type has been done. The mesh of the specimen is shown in Figure 5.21, a constant loading rate is controlled to ensure the quasi-static deformation state which is same with the identification processes. The numerical predicted force-displacement curve compared with the experimental results are shown in Figure 5

Circular cup deep drawing tests (CCD)

The simulations of circular cup deep drawing test for Magnesium alloy AZ31B are performed using FE code Abaqus/Explicit, the objective is to demonstrate the capability of the proposed model for hot sheet metal forming processes. Figure 5.28 shows the schematic representation of die, blank holder and punch [Tari 13]. The dimensions are given in Table 5.10. A near iso-thermal forming of CCD test was conducted in [Tari 13], in which the temperature of the die and holder are 225 °C, the temperature of punch is 215 °C, the blank center at about 211°C. In this work, the heat transfer by radiation and conduction between the blank sheet and tools is ignored for simplicity, only the heat generated from the inelastic deformation is considered. However, to approach the real condition of the experiment, the blank in the numerical simulation was divided into three sections to include the temperature gradient, as illustrated in Figure 5.29. The thermal properties of material used in numerical analysis are given in Table 5.11.

Table 5.11 Thermal properties [START_REF] Tari | Experimental studies of deep drawing of AZ31B magnesium alloy sheet under various thermal conditions[END_REF] Thermal conductivity Heat Capacity Density Inelastic heat fraction 95.8 W/(m K) 1049.3 J/(kg K) 1770.0 kg/m 3 0.9

In the numerical analysis, the blank holder force is 80.0 kN with punch speed of 4 mm/s. The punch, holder and die are rigid body, the blank sheet is meshed with element type C3D8RT. The friction coefficient is fixed to 0.05 for the contacted surfaces. The time step is increased to 5E-5 to reduce the simulation time thanks to the mass scaling method. The force vs displacement curves obtained from the simulation and experiment are compared in Figure 5.30. The simulation result over-predicts the force at large displacement slightly, and the predicted fracture occurs a little earlier compared to the experimental result. Despite all this, the maximum force and the displacement at fracture are relatively well predicted by the proposed model.

Cross-shaped cup deep drawing tests (CSD)

In order to present the capability of our numerical simulation methodology for hot forming processes with complex geometry, the proposed fully coupled damage model was used for the cross-section deep drawing analysis. The geometry of the CSD process is shown in Figure 5.32. No symmetry condition is considered and the simulation was performed with the full part. The element type for the blank is brick element C3D8RT from Abaqus element library and the mesh is shown in Figure 5.32 where different mesh size was set for the blank with the minimum mesh size of 0.35 mm in the region of interest. The punch displacement at fracture under high temperature is highly increased compared to that at room temperature (i.e. 5.5 mm for room temperature and 11.97 mm for temperature of 250°C), as shown in Figure 5.39. The crack path is well predicted by the numerical simulations. Therefore, the proposed fully coupled

Conclusions

In this chapter, the identification methodology for anisotropic parameters, hardening parameters and damage Conclusions and main perspectives 150  Mg alloys often exhibit more complex stress vs. strain response when applying cyclic loading path at room temperature, which is caused by the interaction of slip, twinning, de-twinning deformation modes respect to different loading paths. The current model could not describe this complex behavior which need to be exploited later, the use of multi-surface yield criterion could be feasible.

 The proposed viscous stress could not accurately describe the stress vs. strain response at small plastic strain level. Enhanced viscous stress model by introducing dependency of the viscous module on equivalent viscoplastic strain should be considered.

 The yield criterion used in the current work can describe the SD effect while it is unable to apply for pressure sensitive metals, to obtain a wide application range for different metals of the proposed model, extension by using more sophisticated yield functions to include the pressure sensitivity should be made

      , X    et   , r R  
comme décrit par les équations suivantes :

(

1 e e d      , 1 d      ) (1) ( 1 d      , 1 X X d    ) (2) ( 1 r d r     , 1 R R d     ) (3)
Le paramètre  est un paramètre de couplage de l'endommagement et l'écrouissage isotrope [START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF].

Dans ce travail une nouvelle approche est proposée dans le but d'éviter le problème de continuité induit par la décomposition spectrale pour tenir compte des fermetures de microfissures. L'idée de base consiste à considérer une fonction h de fermeture de microfissures dépendant de la triaxialité des contrintes comme suit :

1 1 ( ) tanh( ) 2 2 c c h h h h        ( 4 )
      et 1 ( ) e H H h d        (5) 1 ( ) eD e h d       et 1 ( ) S S h d     ( 6 ) 1 ( ) h d       et 1 ( ) X X h d     ( 7 ) 1 ( ) r r h d r      et 1 ( ) r R R h d      ( 8 )
Les paramètres e  et r  sont des paramètres de couplage qui contrôlent l'effet de l'endommagement sur la partie hydrostatique de la contrainte et de l'écrouissage isotrope.

                    ( 9 ) 2 0 0 0 2 0 0 0 1 ( )( : ) ( )( : ) 3 ( ) ( )( :1) ( ) 2 2 1 (1 ( ) ) ( )( : ) (1 ( ) ) ( )( : ) 2 3 1 ( ) ( ) ( )( :1) ( ) 2 1 1 ( ) : ( ) 3

T T T T h

d T h d T C h d T T T T T T C T Q T                                                           2 2 1 1 (1 ( ) ) ( ) : (1 ( ) ) ( ) 3 2 r r h dCT h d QTr                        ( 1 0 ) Dans ces équations, ( ) e T  et ( ) e T  sont les constantes de Lamé (   ( ) ( ) / (1 )(1 2 ) e T E T        ) et   ( ) ( ) / 2(1 ) e T E T     , ( ) (3 ( ) 2 ( ))/3 ( )/(
T T T T T T C T Q T r                                       
( 1 1 )

Les relations d'état s'obtiennent par dérivation du potentiel d'état par rapport aux variables d'état comme suit :

  0 (1 ( ) ) ( ) 3 1 ( ) ( ) 1 2(1 ( ) ) ( ) H e e
eH e e e

S eD e h d T h d T T T h d T

                            ( 1 2 ) 0 0 3 1 ( ) ( :1) (T T ) e eH e e C s h d T T                 ( 1 3 ) 2 (1 ( ) ) ( ) 3 X h d C T          ( 1 4 ) (1 ( ) ) ( ) r R h d QTr r         ( 1 5 ) e r Y Y Y Y d          (16)   1 1 0 1 ( ) ( )( : ) ( ) ( ) : 2 3 ( ) ( ) (T T ) ( :1) 2 1 ( )
                         (16a) 1 ( ) ( ) : 3 Y h C T      (16b) 1 2 1 ( ) ( ) 2 r r r Y h d QT r      (16c)

Analyse des dissipations

Il est établi que la dissipation mécanique et la dissipation thermique sont toutes les deux séparément positives, ce qui permet de séparer l'analyse des dissipations en partie mécanique (ou intrinsèque) notée M  et partie thermique notée T  , de sorte que :

0 M T       ( 1 7 ) ( ) 0 T q grad T T        (18) : : vp M D X Rr Yd          
( 1 9 )

Analyse de la dissipation thermique

Le vecteur flux de chaleur q  peut être obtenu à partir du potentiel de Fourrier où k représente le coefficient de conductivité thermique. On arrive alors à :

( ) q kgrad T     ( 2
                                           (21)

Analyse de la dissipation intrinsèque

Pour définir les évolution des variables d'état ( , , , ) ,,,; ) 22)) la contrainte équivalente de type [START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF] peut être retrouvée et en prenant 0 w  la contrainte équivalente classique de Hill peut être obtenue. Pour ce dernier cas si on fixe En suivant les travaux de François [START_REF] François | A plasticity model with yield surface distortion for non-proportional loading[END_REF], une façon simple pour décrire la distorsion de la surface de charge dans le cadre de la thermodynamique consiste à définir l'écrouissage cinématique comme étant l'élément de contrôle de la distorsion. Dans ce travail, le modèle de François a été modifié permettant de mieux décrire le changement de la surface de charge dans la direction orthogonale à la sollicitation [START_REF] Yue | Ductile damage prediction in sheet metal forming processes[END_REF]Badreddine17]. Dans cette approche la partie déviatorique de la contrainte S qui intervient dans les contraintes équivalentes du critère et du potentiel d'écoulement inélastique est remplacée par une contrainte de distorsion d S définie par : Tenseur taux de déformations viscoplastiques :

F X R Y d T  permettant de vérifier à tout moment l'inégalité résiduelle :   1/3 3 3 ( ) 1 ( ) 1 ( ) r c c H J y Z w Z R f T h d h d           ( 2 2 )                  1/3 3 3 ( ) 1 2 0 ( ) 3 ( : ) 4 1 ( ) 1 ( ) 1 ( ) , ( ) , 2 1 ( ) ( ) 1 1 ( ) r r p p H J s T L T L Z w Z R aX X F C T h d h d h d S T Y Y T bR S T Q T h d s T h d                         ( 2
1 2 2 3 2 1 ( , ) ( ) [1 ( )][ ( ) ( ( ) ( )) ( ( ( )))] 2 L L s a p a T Tanh a T a T a T Tanh p p T           ( 2 4 ) 1 1 2 2 3 2 1 ( , ) ( ) [1 ( )][ ( ) ( ( ) ( )) ( ( ( )))] 2 L L s b p b T Tanh b T b T b T Tanh p p T           ( 2
L  ( 1 1 L     ) est défini comme suit : 3 3/ 2 2 3 3 6 2 1 1 a r c c o s ( ) 2 L L J J         ( 2
i i i d i d H Z S X H T S X i c p     ( 2 8 )       3 det det ( ) : ( ) ; ( , ) i i i i d J Z Z H T S X i c p       ( 2 9 ) 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 2 i i i i i i i i i i i i i i i i G H H G H H F F G F F H H N M L                               ' ' ' ' 0 0 0 ' ' ' ' 0 0 0 ' ' ' ' 0 0 0 0 0 0 2 ' 0 0 0 0 0 0 2 ' 0 0 0 0 0 0 2 ' i i i i i i i i i i i i i i G H H G H H F F G F F H H N M L                                (
    0 0 1 0 2 : 2(1 ( ) ) ( ) / 1 ( ) ( ) : 2(1 ( ) ) ( ) / 1 ( ) ( ) r r c d c l y l y S S S S X h d X T R h d T X X S h d X T R h d T                  ( 3 1 )   0 0 1 : 2(1 ( ) ) ( ) / 1 ( ) ( ) r p d p l y S S S S X h d X T R h d T          ( 3
vp p f f F D n          ( )
Tenseur d'évolution de l'écrouissage cinématique:

( ( , ) ) x f f F n a p X              ( ) Evolution de l'écrouissage isotrope :   ( , ) r f f F r n b p r R            ( ) Evolution de l'endommagement ductile : ( ) 0 ( ) ( ) ( , ) (1 ( ) ) s T f f T Y Y T F d Y S T h d                       ( ) avec   0 1 1 : 1 ( ) (1 ( ) ) ( ) / 1 ( ) ( ) r p D A p l y X S n n I h d h d X T R h d T                                 ( )     0 0 1 0 1 ( : ) 2(1 ( ) ) ( ) / 1 ( ) ( ) 1 ( : )( : ) 1 ( ) (1 ( ) ) ( ) / 1 ( ) ( ) ( : ) r r A A p l y x A p l y S S n n h d X T R h d T n S X X n h d S h d X T R h d T X X                                    ( )   3 2 0 0 2 1 ( : )( : ) 1 1 1 ( ) 2(1 ( ) ) ( ) / 1 ( ) ( ) r r A r p l y S S n X n h d h d X T R h d T                       ( 4 0 )   3 2/3 3 3 ( ) : : ( ) 3 T H J A H J w Z H T Z Z Z H T n Z w Z           ( 4 1 )
où l'amplitude de la déformation viscoplastique est défini par f   déterminé par la forme de Norton-Hoff [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF]: A partir de l'équation Eq.( 42) on peut déduire le critère de charge viscoplastique suivant : comme suit :

  ( ) v n T f v f K T    ( 4
  1 ( ) ( , , , , , ) ( , , , , ) ( ) ( , , , , ) 0
v vp v n T v f v f X R d T f X R d T K T f X R d T             ( 4
  ( ) 1/3 3 3 1 1 1 1 1 1 ( 1 ) 1 ( 1) 1 1 ( , , , , , ) ( )
0 1 ( ) 1 ( ) T r n c c H J vp n n n n n n vn n n y n v n n n Z w Z R f X R d T T h d h d                           (44) avec       (
                      1 ( )1 3 n n n S t r        (47)
Les variable force sont données par:

 

                                        (56)                                (60)
                                           
                                                       
                                                                           (74)
Chapitre 4: Etude paramétrique du modèle proposé 4.1 Anisotropie initiale et asymétrie La Figure 4.1 illustre différentes formes de la surface de charge, dues à la variation des paramètres d'anisotropie considérant une asymétrie donnant une plus grande limite élastique en traction (w=1). Il est important de signaler que les paramètres F et F' affectent la surface de charge le long de la direction transverse sans modifier la forme de la surface dans la direction de laminage. Les paramètres G et G' affectent la forme de la surface de charge dans la direction de laminage sans modifier la forme de la surface dans la direction transverse. La variation des paramètres H et H' n'affectent pas la taille de la surface en expansion équibiaxiale et permet la variation de la surface dans les deux directions, de laminage et transverse.

Figure 1 . 1 :

 11 Figure 1.1: Schematics of hot forming process.

Figure 1 .

 1 Figure 1.2 illustrates an example of deep drawing process and its loading conditions. The punch moves towards the metal sheet clamped between the die and the blank holder and deforms it to obtain a desired shape. The external force from the blank holder is to prevent wrinkling and allow the control of the homogeneous material flow. The two principal deformation modes in this process may depend directly on the blank holder action and its corresponding force [Teixeira12]. During this drawing process, the deformation mode changes from shear, compression to tension or from shear, tension to compression (bending/unbending), depending on the location in the sheet thickness.

Figure 1 . 2 :

 12 Figure 1.2: Schematic view of a deep drawing operation and its deformation modes [Kim13].

Figure 1 . 3 :

 13 Figure 1.3: Innovative lightweight construction materials in automotive, aircraft, ship and petroleum industries.

Chapter 1 :

 1 Major physical phenomena exhibited in hot sheet metal forming 8 deformation mode in ideal HCP monocrystal textures.

Figure 1 . 4 :

 14 Figure 1.4: Schematic of the ideal HCP mono-crystal textures with loading direction to activate each deformation mode [Chapuis16].

Figure 1 . 5 :

 15 Figure 1.5: Tension responses of AZ31: (a) Different temperatures with strain rate10 -3 s -1 ; (b) Different

Figure 1 . 6 :

 16 Figure 1.6: Yield stresses and R-values at different orientations of titanium alloy Ti-6Al-4V.

Figure 1 . 7 :

 17 Figure 1.7: Yield loci of Mg alloy AZ31 at different plastic strains [Cazacu06].

(

  ii) ductile shear fracture and (iii) necking, as shown in Figure 1.10. Both diffuse necking and localized necking are two types of instability occurring in sheet metal tensile test before the sheet fracture, as shown in Figure 1.10. Diffuse necking is caused by a reduction of the width over a length. Localized necking occurs inside the diffuse zone because of a reduction of thickness along a narrow band. Sometimes, localized necking can take place without the preceding diffuse necking.

Figure 1 . 10 :

 110 Figure 1.10: Different failure types of sheet metals [Björklund14].

Figure 1 . 11 :

 111 Figure 1.11: Creep strain versus time. As illustrated in Figure 1.11, creep is subdivided into three stages: primary creep where the material hardening is active and is the most important phenomenon, (ii) secondary (or steady state) creep where the hardening is saturated leading to a constant creep strain rate, and (iii) tertiary creep where the creep damage is active and is the most important phenomenon evolving quickly to the final fracture. Creep damage should

  ) Wrinkling Wrinkling generally occurs in the wall or flange of the part, as shown in Figure 1.12a and Figure 1.12b. The flange of the blank undergoes radial drawing stress and tangential compressive stress during the stamping process, which sometimes results in the formation of wrinkles. When the wrinkled flange is drawn into the cup, these ridges appear in the vertical wall. (3) Tearing Tearing usually occurs near the base of the drawn cup leading to a macroscopic open crack in the vertical wall, due to high tensile stresses that cause plastic flow localization, ductile damage occurrence, sheet thinning (where the plane stress condition is no longer valid) and failure of the metal at this location as shown in Figure 1.12c.

Figure 1 . 12 :

 112 Figure 1.12: Sheet Metal defects in drawn parts: (a) Wrinkling in the flange, (b) Wrinkling in the wall, (c) Tearing, (d) Earing.

Figure 2 . 1 :

 21 Figure 2.1: Reference and current configurations of a material body.

Figure 2 . 2 :

 22 Figure 2.2: Schematic illustration of the polar decomposition of F.

Figure 2

 2 Figure 2.3, the elementary force vector t dF is defined at the current configuration t C at point t P , the elementary section t dS  oriented by the normal t n  in t C is transformed from 0 dS  oriented by the normal 0 n 

  2.4. The elastic part of deformation gradient is obtained by unloading from actual configuration t C giving the isoclinic inelastic configuration vp t C (i.e. having the same orientation as the reference configuration 0 C ).

Figure 2 . 4

 24 Figure 2.4 The rotating frame concept and the multiplicative decomposition of the total deformation gradient.

vpF

  is the viscoplastic deformation gradient with respect to the isoclinic configurationvp C . Finally, any tensorial quantity defined in the current configuration is mapped to the locally rotated configurations vp C and C on which all the constitutive equations are expected to be objective. According to the work of[START_REF] Badreddine | On non-associative anisotropic finite plasticity fully coupled with isotropic ductile damage for metal forming[END_REF], by considering Eq.(2-14) as well as the small elastic strain assumption (i.e.

Figure 2 . 5 :

 25 Figure 2.5: Multiscale damage in nature [Bonora17].

Figure 2 . 6 :

 26 Figure 2.6: Different fracture locus: (a) Johnson-Cook model independent of Lode angle; (b) Symmetric fracture locus with Lode angle dependency proposed by Xue and Wierzbicki[Xue08];(c)Asymmetric fracture locus with Lode angle dependency proposed by Bai and Wierzbicki [Bai08].

Figure 2 . 7 :

 27 Figure 2.7: Schematic representation of the damage effect on the force-elongation curve [Saanouni12].

  (effective configuration) is introduced by removing all the voids and cracks. Then we can obtain F A      . Equating the two expressions of F obtained from both configurations, one obtains the following expression for the effective uniaxial stress: bound, d = 0, represents the intact material without any damage, and the upper bound, d= 1, represents complete fracture of the material.

Figure 2 . 8 :

 28 Figure 2.8: Effective stress concept defined in a typical representative volume element (RVE).

Figure 2 . 9 :

 29 Figure 2.9: Yield surface evolution according to different hardening models. -Distortional hardening In order to capture the yield surface distortion (Figure 2.9) during the large plastic flow under nonproportional loading paths, distortional hardening should be taken into account as an induced anisotropy. Attempts to model distortion of the yield surfaces have been made by various approaches in several works [Ortiz83; Voyiadjis90, Kurtyka96, Teodosiu98, François01, Haddadi06, Feigenbaum07, Barlat11]. Teodosiu model [Teodosiu98] can capture the cross hardening effect, but it cannot describe the evolving distortion of

Figure 2 .

 2 Figure 2.10: A schematic view of the two-surface model.(2)Three-surface model proposed by Nguyen et al.[START_REF] Nguyen | A practical constitutive model for AZ31B Mg alloy sheets with unusual stress-strain response[END_REF] Three yield functions are proposed to handle three different deformation mechanism, namely, slip mode, twinning mode and detwinning dominant modes:

  are the yield stresses for different deformation modes. The modeling of the hardening asymmetry is achieved through defining three different isotropic hardening given as below:

e

  ) describing the elasto-plastic flow with the small elastic strain tensor and the Cauchy stress tensor; flux vector and its conjugate force;  ( , r R ) describing the isotropic hardening depicting the change of the yield surface size (radius);  ( , X  ) describing the kinematic hardening i.e. the change of the yield surface center location;  ( , d Y ) describing the isotropic ductile damage. The scalar variable d is an average measure of damage for different directions; with d = 0 for the initial undamaged (or safe) material point, while d = 1 for the fully damaged material point. The strong coupling between the plastic flow with hardening and the ductile damage is performed in the framework of the total energy equivalence assumption [Saanouni94, Saanouni12], leading to the definition of the effective state variables   ,

Figure 2 .Figure 2 . 11 :

 2211 Figure 2.11a, h is set to 2, c h is the critical value of h when  is smaller than a certain negative value, if

(



  are the temperature-dependent Young's modulus and yield stress, ( ) C T and ( ) Q T are the temperature dependent kinematic and isotropic hardening moduli. Also the microcracks closure parameter ( ) c h T can be defined as a function of temperature. The thermal expansion coefficient   , the specific heat at constant volume v C , the density and the Poisson's

Figure 2 .

 2 Figure 2.12 shows how the effective Lame's constant e   evolves with the triaxiality for a given value of damage d=0.3. The values of T e   , C e   and S e   represent its value at tension, compression and shear. With the effect of the microcracks closure effect, the value of e   is increased form tension to compression, which leads to T S C e e e         . This case is also suitable for the other effective damage material properties( , , e C Q    

Figure 2 . 12 :

 212 Figure 2.12: The value of e   varying with the triaxiality.

  , a and b define the nonlinear evolution of mixed isotropic and kinematic hardening. To capture the asymmetric behavior in hardening exhibited by Mg alloys with accuracy, the two hardening parameters a (for kinematic hardening) and b (for isotropic hardening) are

Chapter 2 :L

 2 Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage 49 high values at the beginning of the plastic flow in order to achieve a very fast saturation of hardening (see the first plateau which appears in the stress-strain curve in Figure1.9). This stage is twinning dominant hardening, then their values decrease when the equivalent plastic strain increases and an important hardening rate can be found in this stage. After passing the inflection point, the slip dominant mode is activated and the hardening parameters a and b reduce gradually to a small constant value which leads to another value of the hardening stress saturation.The typical variation of the kinematic hardening parameter a and isotropic hardening parameter b with the normalized Lode angle and plastic strain is shown in Figure2.13. The value of parameter a and b reaches a high value at uniaxial compression ( 1   ) and at low plastic strain, while it stays at a small value at uniaxial tension ( 1 L   ) and high plastic strain level. It should be noted that the calculation of the normalized Lode angle does not depend on the rotation frame as shown below: det( ) det( ) det( ) det( ) det( ) det( ) det( ) det( )

Figure 2 . 13 :

 213 Figure 2.13: Typical variation of parameters a and b versus Lode angle and equivalent plastic strain. It has been proved that the ductility of some materials is dependent on the stress state which is represented by the stress invariants, especially the Lode angle [Bao04]. Different damage and fracture models have included the Lode angle effect [Xue08, Bai08, Cao13, Chbihi17]. In this thesis, the Lode angle effect on the damage is considered. Since the ductile damage parameter S which concerns the material ductility controls directly the value of the equivalent plastic strain at fracture, a new form of dependent function of damage parameter S on the Lode angle is proposed:

  and s  are used to adjust the value of S according to different stress states. The detailed description can be found in Chapter 4.The equivalent relative stress is characterized by the accounting for both second and third invariant of the Hill equivalent stress in order to take into account the Strength Differential effect exhibit by HCP materials (e.g., Mg alloys)[START_REF] Cazacu | A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals[END_REF][START_REF] Yoon | Asymmetric yield function based on the stress invariants for pressure sensitive metals[END_REF] :

  .14, two solids occupying configurations with the volume A  and B  with boundaries A  and B  .

Figure 2 . 14 :

 214 Figure 2.14: Schematization of the contact problem between two solids A and B.

F

  measurement of normal and tangent displacement, (A) are the normal contact force and tangent contact force respectively.

  115) Where  is the Tresca coefficient (0 1)    , this model shows that the friction stress (A) T F  is proportional to the sliding force s F , the description of this model at space ( , ) g T v F is given in Figure 2.15.

Figure 2 . 15 :

 215 Figure 2.15: Schematic representation of Tresca-type friction model.

Figure 2 . 16 :

 216 Figure 2.16: Schematic representation of Coulomb friction model.The tribological tests show that the friction coefficient is dependent on temperature and sliding velocity[START_REF] Grüebler | Temperature dependent friction modeling for sheet metal forming[END_REF]. Moreover, using friction model considering the temperature and velocity dependence of friction coefficient can obtain more accurate results than only using constant friction coefficient in metal forming simulation[START_REF] Klocke | FE-analysis and in situ visualization of pressure-, slip-rate-, and temperature-dependent coefficients of friction for advanced sheet metal forming: development of a novel coupled user subroutine for shell and continuum discretization[END_REF]. Some proposed dependent functions [Juanicotena06, Nayebi08, and Grüebler09] of the friction coefficient are given below:

Figure 2 . 17 :

 217 Figure 2.17: Heat generation during friction.

  volume t  with boundary t  .

Figure 3 .

 3 Figure 3.1 A deformable solid and boundary condition at a given time t .

   are the interpolation functions of the displacement and temperature fields respectively. The velocity and acceleration vectors, as well as the rate of temperature can be obtained by the Chapter 3: Numerical aspects 64 derivative with respect to time:

  two interpolation matrix derived from the derivatives of the interpolation functions with Chapter 3: Numerical aspects 65 respect to the coordinates and, defined by:



  is the global lumped mass matrix,  is the lumped overall thermal capacitance matrix, the overall internal and external thermal force vectors.

F

  Update the state variables, e.g. 1 n   , (n 1) p   , using the local integration scheme(see next section); Chapter 3: Numerical aspects 68 (5) Compute the internal forces   int e and external   e ext F forces using the Eq. (3-20); (6) Select the critical time step for the mechanical problem.  Resolution of the equation of local temperature:

  the hydrostatic and deviatoric part of the trial strain, so the corresponding hydrostatic



  The shear strain is stored as tensor shear strains in VUMAT(

Chapter 3: Numerical aspects 84 

 84 The data are passed in and out in large blocks('nblock'), it has a value of 64 or128. All material points in the same block have the same material name and belong to the same element type.

Figure 3 .

 3 Figure 3.2 shows how the fully damaged Gauss points and fully damaged elements are treated in the implementation in Abaqus/Explicit. Once the unknowns at time 1 n t  are obtained, the value of damage will be compared with the critical damage value c d . If it has not reached the critical damage value, all the statevariables will be updated. If damage reaches its critical value, the values of the state variables will be set to zero. The element in FE code will be deleted, which leads to the micro-crack, with accumulating the microcracks, finally the macro-crack can appear in the structure. This is how the damage is used to describe the cracking inside the structure.

Figure 3 . 2

 32 Figure 3.2 Numerical aspects of the proposed model implementation in FE code.

Figure 4 . 1 :

 41 Figure 4.1: Yield surfaces with different value of parameter w.

Figure 4 . 2 :

 42 Figure 4.2: Effect of the anisotropic parameters on the yield surface with w=1 (plane stress condition): (a) parameter F; (b) parameter G; (c) parameter H; (d) parameter F'; (e) parameter G'; (f) parameter H'.

Figure 4 .

 4 Figure 4.2 illustrates different shapes of the yield surfaces with varying the anisotropic parameters and considering yield asymmetry with larger yield stress in tension (w=1). It is worth noting that the parameters F and F' affect the yield surface shape in the transverse direction as shown in Figures 4.2aand 4.2d without changing the yield stresses in rolling direction. The parameters G and G' affect the yield surface in rolling direction without changing the yield stresses in the transverse direction (Figures 4.2b and 4.2e). Varying parameters H and H' conserve the equi-biaxial stresses and allow changes in both rolling and transverse directions as shown in Figure 4.2c and Figure 4.2f. It is to be highlighted that the sensitivity of parameters

,

  The Figure4.3 shows the evolution of the yield surface with distortion effect after 5% of equivalent plastic strain on different stress planes, namely, the deviatoric plane (  -plane), the principal stresses plane ( 11 22 Each symmetric stress tensor can be expressed by the following six-dimensional unit tensor [Yue14, Badreddine17]:

X

  Figure 4.3b, 4.3d, 4.3f show that with the value of 2 c l X decreasing, the yield surface is expanded. It means that 2 c l X controls the cross size of the yield surface orthogonal to the loading direction. It should be noted that negative value of 2 c l X leads the yield surface be shrink in the orthogonal direction as shown in Figure 4.3d and 4.3f. Another noteworthy feature is the loss Chapter 4: Parametric study of the proposed model 89 of yield surface convexity happening for the case the value of 1 =40.0. In order to keep the convexity of the yield surface, the distortional parameters should fulfil to the condition c deviatoric plane (π-plane) Yield surfaces in principal stress plane Yield surfaces in tension-shear plane

Figure 4 . 3 :

 43 Figure 4.3: Effect of the distortional parameters on the yield surface of the proposed model.

Figure 4 . 4 :

 44 Figure 4.4: Effect of the temperature on the predicted initial yield surfaces with or without considering the tension-compression asymmetry: (a) Hill yield function; (b) proposed model. As shown in Figure 4.4, the size of the yield surface is shrunken by the increase of temperature. Since the Hill criterion can only describe initial anisotropy in tension similar to compression, the yield surfaces are always symmetric between tension and compression and keep the same shape at different temperatures.

Figure 4 .X

 4 Figure 4.5 shows, in the -plane, the different yield surfaces with distortion at different temperatures and 5% plastic strain. The surfaces have a big change when the distortion (lines with hallow symbol in Figure 4.5)

Figure 4 . 5 :

 45 Figure 4.5: Effect of temperature on the distorted yield surfaces.

Figure 4 . 6 :

 46 Figure 4.6: Different loading paths.

Figure 4 . 7 :

 47 Figure 4.7: Evolution of the yield surfaces predicted by the proposed model for AZ31 under two loading paths at two temperatures: (a) loading path I at 65°C, (b) loading path II at 65°C, (c) loading path I at 150°C, (d) loading path II at 150°C. The evolution of the yield surfaces during path I and path II at temperature of 65°C and 150°C with strain rate of 10 -4 s -1 are computed. As shown in Figure 4.7, for each loading path, four yield surfaces at different strain are illustrated to show the evolution of yield surface. During the loading path I in Figure 4.7a and Figure 4.7c, the yield surface has a distinct change in the size, location and shape, which is caused by the coupling of isotropic, kinematic and distortional hardening. A slight rotation of the main diameter of the yield surface to the orthogonal direction is observed in the loading path II (Figure 4.7b and Figure 4.7d).

Figure 4 . 8 :

 48 Figure 4.8: Effect of parameter 1 on the shear stress-strain response.

Figure 4 . 9 :

 49 Figure 4.9: Effect of parameter a2: (a) Stress-plastic strain curves in compression and tension; (b) Evolution of the parameter a in compression; (c) Stress-plastic strain curves in shear.The same observation can be made concerning the isotropic hardening, but with small effects on stressplastic strain curves in compression and in shear as shown in Figure4.10. Due to the lower value of isotropic hardening modulus Q compared to kinematic hardening one C, the effect is less obvious compared to the parameter a2, especially at the starting stage of the plastic strain. For most important class of metallic materials, the values of the parameters we chose to do the parametric study determine that the kinematic

Figure 4 . 10 :

 410 Figure 4.10: Effect of parameter b2: (a) Stress-plastic strain curves in compression and tension; (b) Evolution of the parameter b in compression; (c) Stress-plastic strain curves in shear.

Figure 4 .

 4 Figure 4.11 represents the effects of the parameter a3 for the same loading paths. No effect of this parameter

Figure 4 . 11 :Figure 4 . 12 :

 411412 Figure 4.11: Effect of parameter a 3 : (a) Stress-plastic strain curves in compression and tension; (b) Evolution of the parameter a in compression; (c) Stress-plastic strain curves in shear

Figure 4 . 13 :

 413 Figure 4.13: Effect of parameter p s : (a) Stress-plastic strain curves in compression and tension; (b) Evolution of the parameter a in compression; (c) Stress-plastic strain curves in shear.

Figure 4 .

 4 Figure 4.14 depicts the stress-strain evolution with different value of 2  . A lower curvature is found in compression when 2  decreases, in other words, the hardening rate is smaller with the lower values of 2  . The parameter 2  controls the curvature of the hardening curve in compression. From Figure 4.14(b), the parameter a starts with different values and evolves by decreasing to reach the same value around 150 for p=ps=4% after which this evolution is inverted going to a clear independence of a from 2  after 12% plastic strain. A very limited rise in shear stress can be found with the increase of 2  as shown in Figure 4.14(c).

Figure 4 . 14 :

 414 Figure 4.14: Effect of parameter 2  : (a) Stress-plastic strain curves in compression and tension; (b) Evolution of the parameter a in compression; (c) Stress-plastic strain curves in shear.

Figure 4 .

 4 figure are obtained without microcracks closure effect i.e. h=1). In this figure, as expected, it is observed that the damage evolution is more important in compression than that in tension, so that final fracture occurs at lower plastic strain in compression compared to tension.

Figure 4 . 15 :

 415 Figure 4.15: Effect of SD and HA fully coupled with ductile damage in stress vs plastic strain curves in uniaxial tension and compression: (a) Stress-plastic strain curves; (b) Ductile damage evolution.

Figure 4 . 16 :

 416 Figure 4.16: Effect of combined SD and HA fully coupled with ductile damage considering the microcracks closure effect through parameter h in stress vs plastic strain curves in uniaxial tension and compression:(a) stress-strain evolution; (b) damage energy density release rate.

  Figure 4.16a. When small values of this parameter (h<1) are taken, the damage energy release rate (Y * ) is delayed obviously in compression (see Figure 4.16b) giving rise to a clear delay of the final fracture in compression with lower values of h as indicated in Figure 4.16a. Indeed, it is observedthat the plastic strain at the maximum point of Y * in compression is decreased by increasing the value of h, namely, the final failure is delayed with the decreasing of h in compression due to the fact that damage evolution (or rate) in compression is much lower than in tension.

Figure 4 .

 4 Figure4.17a shows the stress vs strain response and damage evolution of the uniaxial tension test at different temperature with constant strain rate 10 -4 s -1 , with the fully uncoupled and fully coupled models results

Figure 4 .Figure 4 . 17 :

 4417 Figure4.17b shows the strain rate effects on the uniaxial tension test at T=200 °C, with the increase of the strain rate the yield stress increase drastically (lower viscosity effect) while the final accumulated plastic strain at fracture becomes slightly smaller due to the higher damage rate. The two figures show clearly the viscous effects. From the damage evolution at various temperature and strain rate, it is clear that the evolution of damage is delayed by the raise of the temperature and the total strain rate decelerated. So the proposed model could describe the obvious change in ductility with temperature change and the strain rate sensitivity at a constant temperature.

Figure 4 . 18 :

 418 Figure 4.18: Non-proportional loading path. The results of equivalent stress and equivalent kinematic hardening stress are illustrated in Figure 4.19a.

Figure 4 . 19 :

 419 Figure 4.19: Evolution of the equivalent stress and kinematic hardening stress invariant vs the plastic strain for non-proportional loading path III (A-B-C-D-A) for AZ31: (a) comparison of different temperatures fully coupled with damage, (b) four cases at temperature of 65°C. Lines with solid symbols are equivalent stress, lines with hollow symbols are kinematic hardening stress. The yield surfaces at points A, B, C, D (where the loading path changes) of this non-proportional loading path plotted in tension-shear plane are depicted in Figure 4.20 together with the yield surface without

Figure 4 . 20 :

 420 Figure 4.20: Yield surfaces at the end of each loading stage of path III, B: end of stage AB, C: end of stage BC, D: end of stage CD, A: end of stage DA. Grey lines are the yield surfaces without distortion.

Figure 4 . 21 :

 421 Figure 4.21: Detailed evolution of the yield surfaces from point B to point C. A specific evolution of the yield surface from point B to point C is displayed in Figure 4.21, from the red line to the purple line, the clear change is the rotation, once the rotation reach a certain level, more obvious change of the yield surface is the extension in the shear direction, which can be seen from green line to blue

Figure 4 . 22 :

 422 Figure 4.22: Stress state in the space of stress triaxiality() versus normalized Lode angle ( L  ). In this work, the stress triaxiality effect is embedded in the microcracks closure effect with the parameter h depending on triaxiality as shown in Table 4.1 (g). From this dependent function, the desired values of h in tension, shear and compression can be achieved by adjusting the two parameters c h and h  . The evolution of microcracks closure effect parameter versus the stress triaxiality is shown in Figure 2.11 (Chapter 2).

Figure 4 . 23 :Figure 4 . 24 :

 423424 Figure 4.23: Influence of normalized Lode angle parameter L  on the ductility parameter S: (a) for different values of  s (S ten = 2.0, S sh = 1.0); (b) for different values of S ten and S sh ( s = 2.0). The Lode angle effect is taken into account in ductility parameter S which becomes function of the normalized Lode angle as given in Table 4.1 (j). The parameter sh S , ten S and s  are used to adjust the value of S according to different stress states. Figure 4.23 shows the Lode angle dependency of the ductility parameter S, when Sten = 2.0 and Ssh = 1.0, showing that the value of S is symmetric with respect to the vertical axis defined L  = 0, the smallest value of S is obtained for L  = 0, then it increases symmetrically at both the negative and positive sides of L  . The parameter s controls how the value of S evolve with L  , when s increases, the curve becomes more flat at the middle part (corresponding to the small absolute values of L  ). As shown in Figure 4.23b, the evolution curve is convex when Sten > Ssh , otherwise, it is concave when S ten < S sh , and the minor and major values of S can be achieved by varying the values of S ten and S sh .

  parameter h c does not affect the value of fracture strain when 0.33   , a smaller value of hc leads to a high equivalent plastic strain at fracture for the lower stress triaxiality range. The parameters S ten and S sh in the CDM damage model make it possible to describe the Lode angle dependency of fracture strain as reported in [Bao04, Xue07, Bai07, Lou13a].

Figure 4 . 25 :

 425 Figure 4.25: Effect of parameter Ssh on the curves of accumulated plastic strain at fracture vs triaxiality given for d=0.9 considering S ten =0.5 and h c =0.1.

Figure 4 .

 4 Figure 4.25 shows the effect of varying the parameter Ssh (under the condition Ssh <Sten) on the evolution of accumulated plastic strain at fracture vs stress triaxiality. This figure shows a local maximum for positive triaxiality corresponding to both simple tension (

Figure 4 . 26 :

 426 Figure 4.26: Effect of parameter h on the curves accumulated plastic strain vs stress triaxiality for d=0.9 and considering hc=0.1 and Sten=Ssh=0.5.

Figure 4 . 27 :

 427 Figure 4.27: Effect of parameter  S on the curves accumulated plastic strain vs stress triaxiality given for d=0.9 considering Sten=0.5, Ssh=0.1 and hc=0.1.

Figure 4 . 28 :

 428 Figure 4.28: Effect of (S ten -S sh ) on the iso-damage surface displayed in the deviatoric strain plane.

1 L   and minimum fracture strain at the axis of 0 LLL

 10  , with the increase of S, the fracture locus becomes more approaching at the two sides of 1   (as shown in the case  S =9), and a plateau with symmetric axis of 0  appears.

Figure 4 . 29 :

 429 Figure 4.29: Effect of parameter S on the iso-damage surface displayed in the strain deviatoric plane given for d=0.9 and considering Sten=0.5, Ssh=0.1 and hc=1.0.

Figure 4 . 30 :

 430 Figure 4.30: Effect of stress triaxiality on the iso-damage surface displayed in the deviatoric strain plane for d=0.9 and considering S ten = S sh = 0.5 and h c = 0.2.

L

  Parametric study of the proposed model 109 When considering Lode angle effect with Sten=0.5, Ssh=0.1, the shape of fracture locus, like a six-point star, expanded with the decrease of the value of triaxiality, as illustrated in Figure 4.31. However, the expansion of the iso-damage surfaces is no longer isotropic, compared to the surfaces shown in Figure 4.30. The points on the axis with 1 L    move further than the points on the axis 0  with the decrease of the stress triaxiality value. It also shows that the fracture is easy to occur for shear loading ( 0 L   ) under high values of stress triaxiality.

Figure 4 . 31 :

 431 Figure 4.31: Effect of stress triaxiality on the iso-damage surfaces displayed in the deviatoric strain plane for d=0.9 and considering S ten =0.5, S sh =0.1 and h c =0.2.

  Parametric study of the proposed model 110 coupled model and their effects on the ductile fracture locus were carefully illustrated in section 4.8.

  Q a a a b b b p   ) are calibrated based on uniaxial tension and compression tests to fit part A-B of the experimental curve and can be extrapolated to point F (Figure 5.1). Fourth, the damage parameters ) are fitted from the path B-C-D (Figure 5.1) where the damage-induced softening occurs. Finally, the values of hardening parameters and damage parameters can be readjusted in order to get a best fit with the complete nonlinear curve A-B-C-D (Figure 5.1). In summary, we assume that at the Chapter 5: Parameter identification and hot sheet metal forming applications 113 beginning of the plastic flow, the damage effect is small enough so that we can neglect it to only identify the undamaged behavior parameters based on the curve A-B-F. Once the hardening parameters are obtained and fixed, the damage parameters can be identified according to B-C and C-D parts of the stress-strain curve. For the temperature dependent material parameters, first, their values are calibrated with respect to different isothermal test conducted with various constant temperature values. A polynomial function is then used to calibrate these material parameters over the range of temperature under concern. In this work, both analytical and inverse methodologies are applied to obtain the best sets of adequate values of material parameters.

Figure 5 . 1 :

 51 Figure 5.1: Schematic representation of the material parameters identification procedure The identification of anisotropic parameters needs the definition of analytical expressions for yield stress function at different material orientations. In 11 22 ( )   

  using Eq.(5-1) and Eq.(5-2). These yield stresses are expressed in three classical directions 0

  in-plane balanced biaxial tensile T b  and compressive C b stresses are also obtained by substituting the stress components of( , ,0)

  to balance the different experimental results. The objective function is optimized using down-hill simplex method (Nelder and Mead, 1965) , also recently adopted by[START_REF] Li | Anisotropic and asymmetrical yielding and its distorted evolution: modeling and applications[END_REF][START_REF] Yoon | Asymmetric yield function based on the stress invariants for pressure sensitive metals[END_REF]. Under Chapter 5: Parameter identification and hot sheet metal forming applications 115 3D loading paths, four additional anisotropy parameters (namely , , , M L M L   ) need to be calibrated. To achieve this goal, uniaxial tensile and compressive yield stresses in the 11

Figure 5 . 2 :

 52 Figure 5.2: Identification methodology of the material parameters by combining different software.

Figure 5 .

 5 Figure 5.3 displays the predicted yield surface and experimental data points. The proposed yield function is compared to Hill yield function, which only describes the initial anisotropy in tension. It is observed that the proposed model can capture the anisotropy and tension compression asymmetry of titanium alloy Ti-6Al-4V accurately.

Figure 5 . 3 :

 53 Figure 5.3: Yield surfaces predicted by Hill criterion and proposed yield function.

  The predicted yield surfaces and experimental data of two Mg alloys are shown in Figure.5.3. It should be noted that the experimental data of balanced biaxial tests are not available in the work of[START_REF] Khan | Mechanical response and texture evolution of AZ31 alloy at large strains for different strain rates and temperatures[END_REF][START_REF] Habib | Anisotropy, tension-compression asymmetry and texture evolution of a rare-earthcontaining magnesium alloy sheet, ZEK100, at different strain rates and temperatures: experiments and modeling[END_REF], they are assumed to the mean values of three uniaxial yield stresses (

Figure 5 . 3 :

 53 Figure 5.3: Comparison between numerically predicted yield surfaces and experimental data: (a) AZ31; (b) ZEK100 (Solid symbols are experimental data; hollow symbols are mean values).

Figure. 5. 4 .

 4 Figure. 5.4. Comparison of experimental data and simulation results of stress-strain evolution under uniaxial tension, shear and compression: (a) AZ31B; (b) ZEK100.

Figure 5 . 5 Figure 5 . 6 :

 5556 Figure 5.5 Comparison of predicted yield surfaces and experimental measurements of AZ31B at different temperatures. (Solid symbols are experimental data; hollow symbols are mean values) Comparison of stress-strain evolution between simulation and experimental results of AZ31B at four different temperature values at constant strain rate 0.001s -1 are presented in Figure 5.6.

Figure 5 . 7 :

 57 Figure 5.7: Tensile stress-strain evolution of simulation and experimental results of AZ31B along the rolling direction at different temperatures and strain rates.

Figure 5 . 8 :Table 5 . 7

 5857 Figure 5.8: FEM mesh of the notched and shear specimen (Element type C3D8R).

Figure 5 . 9 :Figure 5 . 10 :Figure 5 . 11 :Figure 5 . 12 :Figure 5 . 13 :

 59510511512513 Figure 5.9: Comparison of numerical and experimental response of tensile tests: (a) shear specimen; (b) notch radius 6.67 mm; (c) notch radius 10 mm; (d) notch radius 20 mm. Improvement of the different models is presented by comparison of the results of each one with experimental data. As illustrated in Figure 5.9, the different numerical simulation results and the experimental response of tensile tests are compared. All the predictions made by the three models have discrepancy with the experimental data, our fully coupled damage model (Model3) has a better fitting with experimental results compared to other models, and especially the displacement of the specimen at final fracture is well predicted.

Figure 5 .L

 5 Figure 5.15 shows the normalized Lode angle along the horizontal axis of the notched part, at the two edge sides the value of L is near 1 which is higher than that in the center not exceeding 0.5. Recall the Eq. (2-88) in Chapter 2 as well as the Figure4.28 in Chapter 4, for Case_1 (

Figure 5 . 14 :

 514 Figure 5.14: Comparison of von Mises stress, damage (SDV14) and plastic strain (SDV15) for different cases of damage parameters (Case_1 with

Figure 5 . 16 :

 516 Figure 5.16: Geometry of the notched and shear specimens with dimensions.

  . The predicted results fit well with the experimental data for both notched tensile and shear tests. Small discrepancies are observed for 250°C for these two tests at large displacement level. This seemed caused by the necking which occurs in experiments. With the temperature increase, the fracture occurrence is delayed. The onset of fracture for notched and shear specimen are predicted with a satisfied accuracy.

Figure 5 . 17 :

 517 Figure 5.17: Comparison of numerical and experimental response of tensile tests: (a) notched specimen with radius 4.0 mm; (b) shear specimen.

Figure 5 . 18 :

 518 Figure 5.18: Comparison of the crack path of the shear specimen obtained from simulation and experiments at different temperatures.

Figure 5 . 19 :

 519 Figure 5.19: Plot of elasticity, hardening and damage parameters as function of temperature for AZ31.In this work all the temperature dependent material parameters are expressed by a third-order polynomial function, as given in Eq. (5-17):

Figure 5 . 20 :

 520 Figure 5.20: Variation according to temperature of anisotropic parameters.

  .22. The FE simulation results reported in [Giglio12] and [Allahverdizadeh13] are also presented with the name of 'FEM_01' and 'FEM_02'. The proposed model ('FEM_03') has a better prediction of fracture compared to the result of 'FEM_01', and the hardening evolution are fitted better compared to the result of 'FEM_02'.

Figure 5 . 21 :

 521 Figure 5.21: Mesh of the flat tensile test.

Figure 5 . 22 :Figure 5 . 23 :

 522523 Figure 5.22: Comparison of numerical and experimental force-displacement curves for tensile test on flat specimen for titanium alloy Ti-6Al-4V.

Figure 5 . 24 :

 524 Figure 5.24: Three-point bending test: (a) Test apparatus [Tari 2015]; (b) FE model with mesh Figure 5.25 shows the force-displacement curves obtained by three different simulations compared to experimental results. These simulations are conducted using von Mises model (SIM_01), current model without hardening asymmetry (SIM_02), and current model with hardening asymmetry (SIM_03). A higher level of loading force is observed for the von Mises model (SIM_01), both the current model without hardening asymmetry (SIM_02) and with hardening asymmetry (SIM_03) show a satisfactory agreement with the experimental curve. It is worth mentioning that the Model taking into account the hardening asymmetry (SIM_03) predicts more accurately the force displacement experimental response.

Figure 5 . 25 :

 525 Figure 5.25: Comparison of predicted force vs. displacement responses with experimental data Two elements located at the top (compression state) and bottom (tension state) of the sheet, as shown in Figure 5.26a, are selected to compare their response obtained from different models. As illustrated in Figure 5.26b, the stress vs. plastic strain obtained from von Mises model (SIM_01) is symmetric in yield stress and hardening. The stress vs plastic strain response of SIM_02 and SIM_03 at bottom point (tension state) is almost the same, meanwhile, both simulation results of SIM_02 and SIM_03 can capture the tension

Figure 5 . 27 :

 527 Figure 5.27: Distribution of von Mises stress, damage (SDV14) and plastic strain (SDV15) for TPBT at displacement of 30.0 mm.

Figure 5 . 28 :

 528 Figure 5.28: Schematic representation of die, blank holder and punch.

Figure 5 . 29 :

 529 Figure 5.29: Different temperatures used for the blank sheet.

Figure 5 . 30 :Figure 5 . 31 :

 530531 Figure 5.30: Comparison of force vs. displacement curves for simulation and experimental results for CCD test. The distribution of von Mises stress, damage (SDV14) and plastic strain (SDV15) for CCD test at different veritical punch displacements are given in Figure 5.31. At the displacement of 25.0 mm (Figure 5.31(a)),the damage is relatively low (D max = 2.4%), but is starts to localize at the top edge of the deformed cup. As the displacement increases, more significant damage localization at this place is observed (Dmax = 6.3% at U=40.0 mm; Dmax = 13.6% for U = 48.0 mm; Dmax=96% for U=52.0 mm). It should be noted that a localization at the bottom edge of sheet close to the punch is appeared with a lower magnitude compared to

  Figure 5.32: Numerical model for CSD test. The punch velocity is 1.5mm/s, blank-holding force and pad force are 2.5 kN and 1.5 kN respectively, the temperature of the blank sheet is assumed to isothermal and the contact with the tools will not change the temperature of the blank. The force vs displacement curves of AZ31B simulation results at different temperatures are shown in Figure 5.33. A temperature-induced softening leading to a decrease of the force and increase of displacement is clearly observed as the temperature increases. At room temperature, final fracture appears at small displacement (5.1 mm), which indicates a low formability of CSD test for AZ31B. However, the displacement at final fracture increases with the temperature. Accordingly and as expected the formability is improved at high temperature.

Figure 5 . 33 :Figure 5 . 34 :

 533534 Figure 5.33: Force vs. displacement curves of CSD test at different temperatures.

Figure 5 .Figure 5 . 35 :Figure 5 . 36 :

 5535536 Figure 5.35 shows the distribution of von Mises stress, damage (SDV14) and plastic strain (SDV15) for AZ31B under different temperatures at displacement U=5.1 mm. At room temperature, the cracks appear at the corners of the punch for this displacement, and propagated along the horizontal direction of the blank.Due to the full coupling effects, severe drop of the equivalent stress (as well as all the internal stresses) is

Chapter 5 :Figure 5 . 37 :

 5537 Figure 5.37: Temperature (SDV25) distribution at different initial blank temperatures. For the purpose of comparison, an example of fractured cross-shaped blank sheet of AZ31B [Lee15] formed under non-isothermal condition is given in Figure 5.38, in which the experimentally determined forming height at fracture is 11.5 mm. The numerically predicted final fracture of the CSD test at different values of the initial blank temperature is shown in Figure 5.39. The failure of the blank sheet at different temperatures locates at the corners of the cross-shaped die for all the cases, which is consistent with the experimentally observed crack path as shown in figure 5.38.

Figure 5 . 38 :

 538 Figure 5.38: Example of fractured cross-shaped blank sheet of AZ31B [Lee15].

Chapter 5 :

 5 Parameter identification and hot sheet metal forming applications 146 damage model could predict the location of failure with good accuracy. The simulated forming height at fracture is increased with the temperature, which shows a significant improvement of the formability of AZ31B sheet at elevated temperature. The current model presents well the formability and fracture behavior of AZ31B.

Figure 5 . 39 :

 539 Figure 5.39: Final fracture and punch displacement of AZ31B for CSD test at different temperatures.



  parameters of the proposed model is presented in details. Titanium alloy Ti-6Al-4V are investigated at room temperature, and validated by simulation of flat tensile test, the simulation results show a good agreement with the experimental measurements. For Mg alloy AZ31B, the material parameters are calibrated under different temperatures. Simulations of three-point bending (TPB) tests at room temperature and circular cup deep drawing (CCD) tests at high temperature are performed to validate the capability the proposed constitutive model in sheet metal forming process for Mg alloy. Finally, the simulations of cross-shaped cup Chapter 5: Parameter identification and hot sheet metal forming applications 147 deep drawing (CSD) tests are conducted at different temperatures, which shows a satisfied damage and fracture prediction ability of the proposed model in hot sheet metal forming processes.Conclusions and main perspectives 148Conclusions and main perspectivesSince the constitutive equations are vital in FE simulations of sheet metal forming processes, this work aims to develop the thermodynamically-consistent non-associative constitutive equations fully coupled with isotropic ductile damage for hot sheet metal forming. It should be addressed that the proposed theory embraces many phenomena observed in metal inelasticity. For example, the initial anisotropy, tensioncompression asymmetry, hardening asymmetry are successfully included in the proposed model. Furthermore, the CDM damage model is enhanced by accounting for the Lode angle and stress triaxiality effects. These achievements constitute the main contribution of this work.The proposed anisotropic elasto-visco-plastic constitutive model accounting for mixed nonlinear isotropic and kinematic hardenings strongly coupled with isotropic ductile damage was formulated in the framework of thermodynamics of irreversible processes. Then numerical simulations of sheet metal forming processes were applied. All the contents are presented sequentially in each Chapter, as summarized below: The major phenomena observed in light weight metallic materials during hot sheet metal forming processes are presented and discussed. Anisotropic thermo-elasto-visco-plastic constitutive equations fully coupled with isotropic ductile damage for metallic materials under large strains, are formulated based on continuum damage mechanics framework. The third stress invariant is added to Hill equivalent stress to describe the tension-compression asymmetry exhibited by HCP materials. Distortion of the yield surface is taken into account. The kinematic and isotropic hardening parameters are assumed to be function of the normalized Lode angle and equivalent plastic strain to capture the hardening asymmetry. The introduction of tension-compression asymmetry in yield stress and hardening helps to capture accurately the complex behavior exhibited by Mg alloys. A new approach to include the microcracks closure effect to influence the damage evolution under a wide range of stress triaxiality is proposed.Lode angle-dependent enhanced fully coupled damage model is developed to treat the Lode angle effect on ductile damage evolution. The associated numerical aspects of the proposed model are presented. The constitutive equations have been discretized in time domain using appropriate finite difference scheme and in space domain using FEM. The global resolution scheme as well as the fully implicit local integration scheme of the viscoplastic model are elaborated. The capability of the proposed fully coupled model is presented through a systematic parametric study.The initial anisotropy, tension-compression asymmetry, hardening asymmetry and distortion of theConclusions and main perspectives 149 yield surface are investigated first without coupling with isotropic ductile damage, while the interaction with temperature effect is carefully studied. When fully coupled with ductile damage is considered, significant effects on damage evolution induced by these phenomena are observed. Through the introduction in the fully coupled damage model of stress triaxiality and Lode angle effects, the capability of the proposed model is enhanced tremendously in terms of damage prediction under different loading paths (shear, tension, compression) within a wide range of stress triaxiality. The sensitivity of the newly proposed parameters (hardening parameters, damage parameters) on the macroscopic behavior shows clearly the strong ability of the current model in capturing the complex behaviors of light weight metallic materials. The material parameters identification strategy is presented in detail. The anisotropic parameters, hardening parameters, damage parameters are calibrated successively at certain temperatures, then temperature dependent functions are used to describe the evolution of these parameters with temperature. The procedure is realized by MATLAB-based minimization software interfaced with ABAQUS FE code through the Python script. The desired values of the material parameters are determined through comparing the simulation results to the experimental data by minimizing the objective functions. The identification is applied for two materials (Titanium alloy Ti-6Al-4V at room temperature and Mg alloy AZ31B at elevated temperatures). Three sheet metal forming processes namely three point bending test (TPB), circular cup deep drawing test (CCD) and cross-shaped cup deep drawing test (CSD), are simulated for Mg alloys AZ31B under room and high temperatures using the proposed fully coupled damage model. The high capability of the fully coupled CDM model is validated by the comparisons of the damage initiation and propagation.It bears to emphasize that the pursuit of more accurate prediction in material behaviors with simple numerical implementation is always desired. The current proposed model seeks to include the main physical phenomena under large strains of light weight metallic materials, such as various types of hardening (isotropic and kinematic hardening, distortional hardening), the various anisotropies (Initial anisotropy, tension-compression asymmetry) and fracture by ductile damage under complex loading paths. Besides the merits of the proposed model, it also has limitations: (1) the predicted yield surface shows large equi-biaxial compression stress, which indicated that more advanced yield function is still needed; (2) in order to capture the pronounced anisotropic and asymmetric behaviors of Mg alloy at room temperature, new material parameters are added, when applied for elevated temperatures, the identification work to obtain the material parameters are heavy. Considering the limitation of the model as well as the limited time for this work, further developments are needed to be done in the near future to enhance the proposed model, as summarized below:

Chapitre 1 :

 1 Les principaux phénomènes physiques exhibés en mise en forme à haute températureLes nouveaux alliages légers à hautes propriétés mécaniques, comme les alliages d'Aluminium, de Magnésium et de Titane, sont des matériaux très prometteurs pour la construction mécaniques et sont de plus en plus étudiés durant les dernières années. En effet, la demande industrielle de ces matériaux ne cesse de s'accroitre dans plusieurs secteurs industriels (voir Figure1.1) sous la pression des nouvelles normes de réduction des oxydes de carbone en faveur de la protection de l'environnement. Les alliages d'Aluminium ont déjà été très utilisés dans les secteurs de construction automobile (cadres, panneaux de carrosserie et roues), de construction navale et de l'industrie aéronautique et spatiale. Les alliages de Titane, à cause de leurs modules spécifiques important, ils sont spécialement dédiés aux applications aéronautiques et spatiales pour la fabrication par exemple d'aubes de turbes de moteurs à réaction. Egalement, à cause de leurs excellentes résistances à la corrosion et à la fatigue, ils sont utilisés dans l'industrie pétrolière pour la fabrication des outils de forage.

Figure 1 . 1 :

 11 Figure 1.1: Matériaux légers de construction innovante en industries automobile, navale, aéronautique et pétrolière. Par contre leur utilisation dans le secteur de l'automobile est fortement limitée par leur coût onéreux. Les aciers à haute et à très haute résistance ayant des modules spécifiques très importants comparés aux aciers conventionnels, ont également fait leur preuve dans la construction automobile. Les alliages de Magnésium représentent un candidat intéressant pour le secteur de l'automobile et une alternative aux alliages d'Aluminium et aux aciers à haute résistance pour le rapport important entre leur contrainte à la rupture et leur densité permettant d'avoir des structures plus légères réduisant la consommation de l'énergie. Par contre, leur utilisation industrielle a été limitée par leurs fortes anisotropies initiale et induite [Agnew05, Khan11, Shi13] incluant une asymétrie importante entre traction et compression (effet-SD) [Kelley68, Khan11, Steglich11]. En effet, leur comportement mécanique exhibe un effet important d'asymétrie de l'écrouissage entre traction et compression [Khan11, Kabirian16, Lee17] ainsi qu'une très faible formabilité à température ambiante [Chen03, Kim13, Rodriguez16]. Ces caractéristiques sont surtout reliées à leur structure Hexagonale-compacte (HC). Par conséquent, la modélisation du comportement mécanique complexe de ce type de matériau demeure un défi à relever.

  0 ) L'équation de la chaleur peut être obtenue en combinant cette équation avec la première loi de thermodynamique [Lestriez02,[START_REF] Saanouni | Computational damage mechanics: Application to metal forming simulation[END_REF][START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]. Dans un milieux élastoviscoplastique endommageable elle s'écrit :

  le cadre d'une théorie nonassociée un critère d'écoulement ( , , , ; ) f X R d T  et un potentiel inélastique ( ,



  sont des paramètres caractérisant l'évolution de l'endommagement, a et b sont des paramètres qui définissent la non linéarité de l'écrouissage cinématique et isotrope respectivement. Pour rendre compte de l'asymétrie de l'écrouissage exhibé par les alliages de Magnésium, ces deux paramètres sont définis comme fonction de l'angle de Lode normalisé L  ( 1 1 L     ) ainsi que de la déformation plastique équivalente p comme suit :

  1

  matériaux. L'angle de Lode normalisé

  l'amplitude de la distorsion du critère de plasticité et du potentiel plastique respectivement. Le paramètre 2 ( ) c l X T affecte la taille de la surface de charge dans la direction orthogonale de chargement. Pour assurer la continuité et la convexité de la surface de charge les paramètres de l'écrouissage cinématique doivent vérifier la condition ( la viscoplasticité, les variables flux qui définissent l'évolution des phénomènes dissipatifs peuvent être postulés comme suit[START_REF] Lemaitre | Mécanique des matériaux solides[END_REF][START_REF]Mechanics of anisotropic materials[END_REF]:

  2 ) Le scalaire positif f   peut être admis comme un multiplicateur viscoplastique en comparaison avec la plasticité indépendante du temps. Par contre ce multiplicateur n'est pas déterminé par une condition de consistance comme pour le cas de la plasticité indépendante du temps.

  3 ) où la contrainte visqueuse est définie par la forme scalaire caractérisant l'effet visqueux et sont dépendant de la température.Chapitre 3: Aspects numériquesLe modèle développé au chapitre 2, décrit par un système d'équations fortement non linéaires, a été implémenté dans le code de calcul par E.F. ABAQUS/Explicit à travers la routine utilisateur VUMAT. L'ensemble des équations ont une forme classique d'équation différentielle ordinaire d'ordre 1 sauf les équations d'écrouissage cinématique et isotrope qui ont une forme spécifique qui se prête bien à une intégration numérique par schéma d'Euler purement implicite et un schéma asymptotique[START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]

  expression de l'endommagement en fin du pas de temps:

Figure 4 . 1 :

 41 Figure 4.1: Effet des paramètres d'anisotropie sur la forme de la surface de charge en prenant w=1 (cas des contraintes planes): (a) paramètre F; (b) paramètre G; (c) paramètre H; (d) paramètre F'; (e) paramètre G'; (f) paramètre H'.

Figure 4 . 2 : 1 L

 421 Figure 4.2: Effet de la température sur la prédiction de la surface initiale de plasticité avec et sans prise en compte de l'asymétrie: (a) Critère de Hill classique; (b) Modèle proposée. Comme illustré dans la Figure 4.2, la taille de surface est réduite en augmentant la température (sorte d'adoucissement thermique). Le critère de Hill ne permet de décrire que l'anisotropie initiale en traction ou en compression tout en gardant une forme symétrique en traction et en compression. La Figure 4.2b montre l'évolution de l'asymétrie traction-compression en fonction de la température. A la température ambiante l'asymétrie est significativement importante qui s'estampe lorsque la température croît. La comparaison des résultats des deux modèles montre un intérêt important du modèle proposé vis-à-vis de la description de l'évolution de l'asymétrie avec la température.

Figure 4 . 3 :

 43 Figure 4.3: Effet du paramètre a 2 sur les courbes contrainte-déformation viscoplastique: (a) Traction et compression (b) Evolution du paramètre a en compression; (c) Cisaillement.

Figure 4 . 4

 44 Figure 4.4 Effet du paramètre a 3 sur les courbes contrainte-déformation viscoplastique (a) Traction et compression; (b) Evolution du paramètre a en compression; (c) Cisaillement.

Figure 4 . 5 :

 45 Figure 4.5: Effet du paramètre p s sur les courbes contrainte-déformation viscoplastique (a) Traction et compression; (b) Evolution du paramètre a en compression; (c) Cisaillement.

Figure 4 . 6 :  et minimale en cisaillement 0   et maximale en cisaillement 0 

 4600 Figure 4.6: Effet du paramètre 2  sur les courbes contrainte-déformation viscoplastique (a) Traction et compression; (b) Evolution du paramètre a en compression; (c) Cisaillement. 4.3 Effet de la triaxialité des contraintes et de l'angle de Lode sur l'évolution de l'endommagement Dans ce travail, l'effet de triaxialité des contraintes est introduit dans le paramètre de fermeture de microfissures ( ) h  . Les valeurs souhaitées du paramètre ( ) h  en traction et en compression peuvent être ajustées en pilotant les paramètres c h et h  . L'effet de l'angle de Lode étant introduit dans le paramètre de ductilité ( ) S  , les paramètres sh S , ten S et s  sont utilisés pour ajuster les valeurs de S pour différents état des

Figure 4 . 7 :

 47 Figure 4.7: Effet de (S ten -S sh ) sur la surface d'iso-endommagement exprimée dans le plan déviateur des déformations. L'effet de la triaxialité sur la définition de la surface d'iso-endommagement dans le cas (S ten = S sh =0.5) est montré sur la Figure 4.8a. Dans ce cas la forme circulaire est obtenue donnant la même déformation plastique à rupture quel que soit l'angle de Lode dû à l'état des contraintes. L'augmentation de la triaxialité des contraintes induit une réduction de la taille de la surface, ce qui conduit à une augmentation de la déformation viscoplastique à rupture pour des valeurs négatives de la triaxialité des contraintes. Lorsque l'effet de l'angle de Lode est considéré avec Sten=0.5> Ssh=0.1, la forme en étoile à six branches de la surface d'iso-endommagement gonfle homothétiquement en décroissant la triaxialité des contraintes comme montré dans la Figure 4.8b.

  La comparaison des surfaces initiales de plasticité prédites par le modèle les points expérimentaux de AZ31B à différentes température est réalisée dans la Figure5.1. Cette figure montre clairement que l'effet de l'asymétrie est important à température ambiante et le modèle arrive à le prédire correctement. Il est à noter également que plus la température augmente plus l'effet de l'asymétrie est estampé et cette variation est correctement prédite par le modèle proposé.

Figure 5 . 1

 51 Figure 5.1 Comparaison des surfaces initiales de viscoplasticité prédites par le modèle et les mesures expérimentales de l'alliage AZ31B à différentes températures.

Figure 5 . 2 :

 52 Figure 5.2: Comparaison des courbes issues du modèle et obtenues expérimentalement en termes de contrainte en fonction de la déformation viscoplastique équivalente -de l'alliage AZ31B la longue direction de laminage à différentes températures et à vitesse de déformation 0.001s -1 . Dans la Figure 5.2 sont comparées des courbes en termes de contrainte en fonction de la déformation viscoplastique équivalente entre prédiction du modèle et résultats expérimentaux de l'alliage AZ31B à différentes températures et vitesse de déformation constante 0.001s -1 . Cette figure montre que pour des déformations faibles la prédiction du modèle surestime le résultat expérimental. Ceci est due essentiellement à la forme de la contrainte visqueuse utilisée dans le modèle qui est constante pour une vitesse de chargement donnée quel que soit le niveau de la déformation plastique. Notons également que l'asymétrie de

Figure 5 . 4 :

 54 Figure 5.4: Schématisation du procédé d'emboutissage d'un godet cylindrique. Les conditions de chargement du procédé sont définies par une force serre-flan de 80.0 kN et une vitesse de déplacement du poinçon de 4.0 mm/s. Le coefficient de frottement est fixé à 0.05 pour toutes les surfaces en contact. Dans la Figure 5.5 sont comparées les courbes force-déplacements issues de la simulation et de l'expérience. Nous pouvons constater que la courbe numérique surestime légèrement l'effort maximal d'emboutissage et par conséquent la rupture intervient légèrement en avance par rapport aux résultats expérimentaux.

Figure 5 . 5 :Figure 5 . 6 :Figure 5 . 7 :

 555657 Figure 5.5: Comparaison entre simulation et expérience en termes de courbe force -déplacement pour l'essai d'emboutissage du godet cylindrique. Les distributions des champs de contrainte de von Mises, d'endommagement et de déformation plastique

Figure 5 . 8 :Figure 5 . 9 :Figure 5 . 10 :Figure 5 . 11 :procédés de mise en forme de tôles minces à hautes températures sont réalisées et leurs comparaisons aux résultats expérimentaux sont présentées

 5859510511 Figure 5.8: Courbes force -déplacement de l'essai d'emboutissage en croix à différentes températures. Les courbes en termes de force-déplacement obtenues pour l'essai d'emboutissage en croix sur l'alliage AZ31B à différentes températures sont données dans la Figures 5.8. Cette figure montre que le niveau d'effort décroit considérablement avec l'élévation de la température. A température ambiante, la rupture finale de l'embouti intervient avec un faible déplacement du poinçon autour de 5.1 mm (faible profondeur de l'embouti) qui caractérise une faible formabilité du matériau à cette température. Ce déplacement à rupture croît en fonction de la température ce qui indique que la formabilité est améliorée par l'élévation de la température.

Chapter 2 Theoretical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage

	Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage Chapter 2: Theorectical modeling of thermo-elasto-visco-plasticity coupled with isotropic ductile damage
	2.6.4 Thermodynamical admissibility
	2.7 Contact and friction
	2.7.1 Kinematics of contact
	2.7.2 Friction model
	2.7.3 Heat generation from friction
	2.8 Conclusions
	Content
	2.1 Introduction
	2.2 Fundamentals of finite viscoplasticity
	2.2.1 Homogeneous transformation and deformation gradient
	2.2.2 Polar decomposition of deformation gradient
	2.2.3 Strain and strain rate tensors
	2.2.4 Stress tensors
	2.2.5 Kinematics of elasto-viscoplastic decomposition
	2.3 Fundamental laws of thermodynamics
	2.3.1 Mass conservation law
	2.3.2 Momentum conservation law
	2.3.3 The first law of thermodynamics
	2.3.4 The second law of thermodynamics
	2.3.5 The Clausius-Duhem inequality
	2.4 Modeling of ductile damage
	2.4.1 Physical aspects of ductile damage
	2.4.2 Classification of ductile damage models
	2.4.3 The concept of effective state variables in CDM
	2.5 Time-independent plasticity and viscoplasticity theories
	2.5.1 Time-independent plasticity
	2.5.2 Thermo-viscoplasticity (time-dependent plasticity)
	2.6 Formulation of the proposed fully coupled model
	2.6.1 State variables and effective state variables
	2.6.2 State potential and state equations
	2.6.3 Dissipation analysis
	16

Table 2 .1

 2 Relations between different stress measures

	Cauchy

  The continuum damage mechanics (CDM) is based on the thermodynamics framework leading to describe the ductile damage by a scalar or tensorial variable[START_REF] Lemaitre | Mécanique des matériaux solides[END_REF][START_REF] Lemaitre | A Course on Damage Mechanics[END_REF][START_REF] Voyiadjis | A plasticity-damage theory for large deformations of solids. Part I: theoretical formulation[END_REF]

	Saanouni03, Lemaitre05, Lemaitre09, Besson10, Saanouni12, Murakami12]. Due to the strong
	coupling between the ductile damage and the material behavior, the material properties are deeply
	affected by the damage occurrence. The concept of effective stress was first introduced by Kachanov
	[Kachanov58] to define the concept of continuum damage mechanics. Since that, this kind of CDM
	approaches have been used for damage prediction in various situations assuming various types of
	coupling in many published works [Lemaitre85, Lemaitre92, Voyiadjis92, Saanouni03, Lemaitre05,
	Lemaitre09, Besson10, Saanouni12, Murakami12, Baddridine10, Bouchard11, Soylansan16,
	Wulfinghoff17]. Within the CDM approach, isotropic damage assumes that spherical micro-void
	cluster are homogeneously distributed without any privileged evolution direction. While the

description of the anisotropic damage assumes that the microcracks with various shapes and orientations are highly influenced by the initial microstructure of the material and its evolution (texture) as well as by the direction of the applied loading paths. Many extensions have been made based on this framework in order to avoid some limitations, such as nonlocal formulation to avoid the mesh dependencies of the local models

[START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF][START_REF] Saanouni | Micromorphic approach for finite gradientelastoplasticity fully coupled with ductile damage: Formulation and computational aspects[END_REF]

  In this DE resolution scheme neither iterations nor tangent stiffness matrix are required, however the accuracy and the stability of the solutions are highly conditioned by the time step amplitude t  is the highest eigenvalue, and  is the damping parameter. In practice, the estimation of the largest eigenvalue are determined by the highest value of volume expansion of all elements. In ABAQUS/Explicit, the stable time increment is set to be the minimum value of the mechanical step and thermal step which is expressed below: First, the mechanical problem is solved over the time interval under isothermal condition (i.e.  . The resolution of the mechanical step and thermal step are in the following manner:

													Chapter 3: Numerical aspects
													. Accordingly,
	attention should be paid to the control of time increment step t  to ensure the stability and accuracy of the
	solution [Hughes87]. The stable time increment is defined:
													max 2 ( 1 t w  	2    	)	(3-33)
	where max w
													t  	min(	u   t	e d L C	,	T t  	2 e L 2 C 	)	( 3 -3 4 )
	where e
	element, approximated by	d C		(  e		2 )/ e  	and the material diffusion can be expressed by:
						/ (				)
	increment 	t	n	,	t	n	1	t	n	t	
	n T	1	n T
													n T  at 1 n t 
	terminating the resolution of the thermomechanical problems over the same time increment
		t	n	,	t	n	1	t	n		t		final
	time 1 n t	t	f
				 Resolution of the equation of displacement field:
					(1) Calculate the lumped mass matrix	 	M	e	  at time n t ;
					(2) Obtain the values of the variables at time n t (	....
													67

L is the smallest mesh element dimension. d

C

is the dilatational wave speed travelling in the

v C k C   

, where  is the material density, e

 and e  are the Lame's constants.

In this work, the mechanical and thermal problems are solved sequentially over the same typical time    assuming the weak thermomechanical coupling (See

[START_REF] Saanouni | Damage Mechanics in Metal Forming: Advanced Modeling and Numerical Simulation[END_REF]

for more details).

  ) in order to obtain the mechanical variables at 1 n t  . Once the stress, hardening and damage variables are known at 1 n t  , the thermal problem is then solved in order to compute the temperature 1    . This incremental procedure will not be stopped unless the total time exceeds the  n u u 

  ,

	n y t t y	n	n y t n y	n y t y	n	n n t t t

  ( )

	H n	n	e	en	e n		v p
		e e n n	H n	n	e	e	n	n	n

  ]:

	n	p	h d	n	e	e	eH nn T	tr	p n n	h d	n	e	eD nn T	p n n
		h	n	e	e		n		n	n		p n

Table 4 .4

 4 Model parameters for the AZ31with experimental results from[START_REF] Khan | Mechanical response and texture evolution of AZ31 alloy at large strains for different strain rates and temperatures[END_REF]

	T(°C)	H	G	F	N	H'	G'	F'	N'	L=M	L'=M'
	65	1.203704 0.67226	0.539721 3.335803 2.286705 0.219813 0.208985 2.968164	1.5	1.5
	150	0.905765 0.627792 0.423295 2.378798 2.470310 0.138736 0.113971 2.544293	1.5	1.5
		E (GPa)	υ	σ y (MPa)	C(MPa)	a	Q(MPa)	b	X	1 c l	X	l	1 p	X	2 c l
	65	42.0	0.35	145.0	3117.7	50.5	376.9	6.4	60	60	100
	150	32.18	0.35	110.0	1118	50.5	176.9	6.4	30	30	30

Table 4

 4

		.5 Assumed material parameters
	H	G	F	N	H'	G'	F'	N'	L=L' M=M'	σ y
	1.2255 0.48639 0.42001 3.44122 2.70885 0.15525 0.13115 3.49908 C Q a1 a 2 a 3 b 1 b2 b 3	1.5 ps	1.5 1 	161.6 2 
	3000	500	50.0	100.0	35.0	6.0	10.0	3.0	0.04	3.0	40.0

Table 4

 4

		.6 Assumed damage parameters
	S	s	γ	β	Y 0 (MPa)	h
	2.2	3.0	3.0	2.0	0.0	0.2-1.0

Table 5

 5

	.1 Calibrated anisotropic parameters of Ti-6Al-4V
	F	G	H	N	F	G	H 	N 
	0.51268	0.57484	0.37089	1.53660	-0.42908	-0.19191	0.72270	0.63820

Table 5

 5

	.2 Anisotropic parameters of AZ31 and ZEK100
	L=M= L'=M'=1.5
	Material	H	G	F	N	H'	G'	F'	N'
	AZ31	1.2255	0.48639	0.42001	3.44122	2.70885	0.15525	0.13115 3.49908
	ZEK100 1.030511 0.420433 1.607213 3.44436 0.287824 1.473031 -0.26309 0.39605

Table 5 .

 5

3 Material parameters of AZ31

Table 5 .5

 5 Anisotropic parameters of magnesium alloy AZ31

	L=M=L'=M'=1.5
	T (°C)	H	G		F	N		H'		G'		F'	N'
	RT	1.205501 0.496391 0.418510	3.441221	2.608850 0.155251 0.160252 3.499081
	150	0.477407 0.773161 0.621849	1.652781	1.332151 0.193945 0.175041 1.460198
	200	0.474259 0.555571 0.467317	1.463941	0.427075 0.230738 0.326213 0.608006
	250	0.480181 0.557180 0.444049	1.483791	0.390420 0.232581 0.373659 0.437981
	Table 5.6 Hardening parameters of magnesium alloy AZ31
	T (°C) E(GPa)	y  (MPa) Kvp m Q	C	b1	b2	b3 a1	a2	a3	1 	2 	ps
	RT	43.05	161.6	5	3 800 3500 3.0 30 2.5 35 200 32	10	100 0.035
	150	33.68	112.6	15 3 580 1950 3.0 20 1.2 35 90	25	10	50	0.011
	200	30.50	73.49	24 3 300 1150 3.0 10 1.0 35 40	23	10	30	0.0
	250	27.40	58.39	25 3 200 800	3.0 5.0 0.9 35 10	22	10	10	0.0

Table 5 .8

 5 Damage parameters for AZ31B

	T(°C)	Ssh	S ten	ξs	s	γ	β	Y0(MPa)	h
	RT	1.6	3.0	1.0	2.0	4.0	1.0	0.0	0.2
	150	3.0	5.7	1.0	1.0	4.0	1.0	0.0	0.2
	200	5.0	8.9	1.0	1.0	4.0	1.0	0.0	0.2
	250	8.0	12.5	1.0	1.0	4.0	1.0	0.0	0.2

Table 5 .9

 5 Fit functions for anisotropic parameters Fit functions for anisotropic parameters

	7.262e-08*  N ?	T	3		8.531e-05*	T	2		0.02692*	T		4.01 6
	H			1.331e-07 * 	T	3		8.169e-05*	T	2		0.01634*	T		.54 1
	G		3 2.877e-07 * -0.0001401* T	T	2		0.01866*	T		0.1 363
	F		3 1.876e-07 * -9.343e-05* T	T	2		0.01271*	T		0.1 724
	N	'		3 6.257e-07* -0.000239* T	T	2		0.008719*	T		. 7 3 4 1
	H	3 8.018e-07* -0.0003343* ' ? T 	T	2		0.02636*	T		2.15
	G	'		3 -2.761e-08* +1.217e-05* T	T	2		0.001058	*	T		0. 42 1 7
	F	'		3 -8.17e-08* +4.046e-05* -0.004371 2 T T	*0 . 2 6 8 2 T 

Table 5 .

 5 10 Dimensions of the component for CCD test

extensif en Français Introduction Ce

  travail de thèse est dédié au développement d'une modélisation phénoménologique thermodynamique pour une meilleure description du comportement thermomécanique complexe avec endommagement ductile en mise en forme par grandes déformations à hautes températures de tôles métalliques. Ce modèle permet de prendre en compte plusieurs couplages forts entre les déformations thermoélastiques, les grandes déformations inélastiques, les écrouissages (isotrope et cinématique), les échanges thermiques et l'endommagement ductile. En se basant sur les compétences en modélisation théoriques et numériques acquises au sein de l'équipe LASMIS, un modèle anisotrope thermo-élasto-viscoplastique couplé à l'endommagement ductile isotrope a été proposé dans ce travail. L'effet de la fermeture des microfissures, de la triaxialité des contraintes ainsi que l'angle de Lode, ont été pris en compte dans le modèle. L'asymétrie de la surface de charge a été prise en compte en introduisant la dépendance au 3 ième invariant des contraintes Cette méthode a nécessité le développement d'un programme Matlab couplé au code de calcul ABAQUS/Explicit via des scripts en Python[START_REF] Yue | Identification of fully coupled anisotropic plasticity and damage constitutive equations using a hybrid experimental-numerical methodology with various triaxialities[END_REF].

	[Chaboche06, Saanouni12, Lou13b, Yoon14].
	[Wu15]	Wu, Z., & Curtin, W. A. (2015). The origins of high hardening and low ductility in
	magnesium. Nature, 526(7571), 62. Une étude paramétrique exhaustive a été conduite sous divers trajets de chargement simples et complexes
	[Wulfinghoff17] pour vérifier les capacités prédictives du modèle proposé. Enfin, des applications à des procédés de mise en Wulfinghoff, S., Fassin, M., & Reese, S. (2017). A damage growth criterion for anisotropic
	damage models motivated from micromechanics. International Journal of Solids and forme à haute température de tôles minces sont réalisées et les résultats des simulations numériques sont Structures, 121, 21-32. comparés aux résultats expérimentaux s'ils sont disponibles.
	La thèse est organisée en cinq chapitres.

 Only anisotropic plastic behaviors are taken into account in this work, while ductile damage can be highly anisotropic in sheet metal. The proposed model could be extended to couple the anisotropic damage using existing anisotropic CDM damage model developed by our team

[START_REF] Nguen | Anisotropie de l'endommagement et simulation numérique en mise en forme par grandes déformations plastiques[END_REF][START_REF] Rajhi | Anisotropic ductile damage fully coupled with anisotropic plastic flow: Modeling, experimental validation, and application to metal forming simulation[END_REF][START_REF] Badreddine | Damage anisotropy and its effect on the plastic anisotropy evolution under finite strains[END_REF]

.

 To obtain more accurate damage prediction in various metal forming processes, adaptive meshing methodology should be taken into account

[START_REF] Bouchard | Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria[END_REF][START_REF] Labergere | Numerical simulation of continuous damage and fracture in metal-forming processes with 2D mesh adaptive methodology[END_REF]

. To avoid the dependence of the numerical solution to the discretization aspects, the current fully coupled local model should be extended to the non-local formulation using the generalized framework of micromorphic theory [Saanouni12, Diamantopoulou17, Liu17].  Systematic experimental works for Mg alloys or Titanium alloys are required to complete the identification of the proposed model correctly. Meanwhile, this model need to be applied for much more materials (i.e. ZEK100, pure Ti) to extend the applicability of the present model. W X [Xue08] Xue, L., & Wierzbicki, T. (2008). Ductile fracture initiation and propagation modeling using Résumé J3 dans le critère d'écoulement de Hill. Egalement la dissymétrie de l'écrouissage a été représentée par l'introduction, dans les paramètres de vitesse de saturation de l'écrouissage, d'une dépendance à l'angle de Lode. La formulation du modèle a été réalisée dans le cadre de la thermodynamique des processus irréversibles en utilisant une théorie générale de plasticité non associée à normalité non-associée en grandes déformations inélastiques. Les aspects numériques, associés à l'implémentation du modèle dans le code de calcul par E.F. ABAQUS/Explicit via la routine utilisateur VUMAT, ont été développés. La détermination des paramètres d'anisotropie est assurée par la minimisation d'une fonction objective qui représente l'écart au sens des moindres carrés entre la prédiction du modèle et la mesure expérimentale [Yoon14]. Pour identifier les paramètres d'écrouissage et d'endommagement, une méthode inverse a été mise en oeuvre basée sur la minimisation de l'écart entre les réponses du modèle et les résultats expérimentaux en termes de courbes force-déplacement.

Chapitre 2: Formulation théorique du modèle thermo-élasto-viscoplastique couplé à l'endommagement ductile isotrope 2.1 Formulation du couplage du modèle proposé

	2.1.1 Variables d'état et variable effectives
	Dans ce travail, le model est formulé dans le cadre de la thermodynamique des processus irréversibles
	avec variables internes [Lemaitre90, Saanouni12, Skrzypek13]. Les couples de variables observables et
	internes sont les suivantes :
	 ( ,   ) pour la déformation totale et le tenseur des contraintes de Cauchy,
	 ( , e T s ) pour la température absolue et l'entropie spécifique,
	 ( ,
	 ( / q T 	,	( ) g grad T   	) représentant le vecteur flux de chaleur et sa force conjuguée,

e   ) représentant les déformations élastoplastiques avec l'hypothèse des petites déformations élastiques et le tenseur des contraintes de Cauchy,  ( , r R ) représentant l'écrouissage isotrope qui représente la variation du rayon de la surface de charge,  ( , X  ) représentant l'écrouissage cinématique qui représente la variation de la position du centre de la surface de charge  ( , d Y ) représentant l'endommagement ductile isotrope. La variable scalaire d est une moyenne de tous les endommagements dans différentes directions, avec d=0 pour un élément de volume initialement saint et d=1 pour un élément de volume complétement endommagé. Le couplage fort entre l'écoulement plastique, l'écrouissage et l'endommagement ductile est réalisé dans le cade de l'hypothèse d'équivalence en énergie totale [Saanouni94, Saanouni12], permettant de définir les couples de variables effectives   , e

    sont les parties hydrostatique et déviatrice de la déformation d'essai élastique.

					(	vp   	0;			  	0;	r  	0;		0   )  ( 1 d vp n  			vp n	;		n	1 			n	;	n r	1 		; r d n n	1 	n d  )	(64)
	La déformation d'essai élastique à 1 n t  est donnée par :
																			* 1  e n			* 1 eH n 			eD n 	* 1		e n       eH eD HD n n          	( 6 5 )
	dans laquelle	* 1 eH n   et		* eD 1 n
	Les parties hydrostatique	* 1   et déviatrice * H n n S  de la contrainte sont données par : 1
													* n S    1 2(1 ( ) ) ( )( n en h d T   	eD n	  	D	)	2(1 ( ) ) ( ) n en h d T      n S  	D	( 6 6 )
										* 1  H n	(1 ( ) ) ( ) tr( (1 ( ) ) ( ) tr( ) 3 1 ( ) ) e e e n en n H n n en h d T I h d T I h d                    	n 	e		 T T T ( ) en n n 	1 				I	( 6 7 )
	D'où la contrainte d'essai élastique aura la forme suivante:
									S	* n	1														H n	* 1
		* n	1	2(1 ( ) ) ( ) n e n h d T	* 1 eD n	(1 ( ) ) ( ) e n e n h d T	* 1 eH n	3 1 ( ) h d	n	e	e	 T T T ( ) n n n	1		I
							n
					S	n	H n	2(1 ( ) ) ( ) n e n h d T				D	(1 ( ) ) ( )tr( ) 3 1 ( ) e n e n h d T I h d
																							  	  	 	( )
	avec			eH	H   		vpH	(    tr 		vp	)	et			eD	D   		vpD
	Durant la prédiction élastique nous supposons que l'incrément de déformation imposé est purement
	élastique et pour une condition isotherme de température 1 n T   ne produisant donc aucune dissipation ce n T
	qui permet d'avoir :

Table 5 .1

 5 Paramètres d'anisotropies identifiées pour différentes températures de l'AZ31

	L=M=L'=M'=1.5
	T (°C)	H	G		F	N		H'		G'		F'	N'
	RT	1.205501 0.496391 0.418510	3.441221	2.608850 0.155251 0.160252 3.499081
	150	0.477407 0.773161 0.621849	1.652781	1.332151 0.193945 0.175041 1.460198
	200	0.474259 0.555571 0.467317	1.463941	0.427075 0.230738 0.326213 0.608006
	250	0.480181 0.557180 0.444049	1.483791	0.390420 0.232581 0.373659 0.437981
	Table 5.2 Paramètres d'écrouissage identifiés à différentes températures de l'AZ31
	T (°C) E(GPa)	y  (MPa) Kvp m Q	C	b1	b2	b3 a1	a2	a3	1 		2	ps
	RT	43.05	161.6	5	3 800 3500 3.0 30 2.5 35 200 32	10	100 0.035
	150	33.68	112.6	15 3 580 1950 3.0 20 1.2 35 90	25	10	50	0.011
	200	30.50	73.49	24 3 300 1150 3.0 10 1.0 35 40	23	10	30	0.0
	250	27.40	58.

1

)

Correction viscoplastique avec endommagement:

 les variables d'état d'essai élastique doivent être corrigées à l'instant 1

 . Les équations discrétisées du système non linéaire peuvent être réduites à uniquement 4 équations fonctions des 4 variables L'accroissement de la température induit un retardement de la rupture du à l'endommagement ductile.

L'amorçage de l'endommagement pour les essais de traction entaillés et de cisaillement sont prédits correctement par le modèle

Applications à des procédés de mise en forme