Abderrahman Clément

Axel Lili

Zacharie Nathan

Nicolas

Deyo

Jean-Baptiste David

Manon Paul

Daria Un

Arrivé à la fin de mes trois riches et passionnantes années de doctorat au CEA Valduc et au Centre des Matériaux de Mines Paris, je tiens à remercier ici l'ensemble des personnes sans qui ce travail n'aurait pu avoir eu lieu.

Je remercie tout d'abord l'ensemble des membres du jury pour avoir accepté d'examiner mon travail. Merci au Pr. Kondo de m'avoir fait l'honneur de présider le jury. Merci aux Prs Fritzen et Abed-Meraim pour leurs rapports détaillés et inspirants sur le manuscrit, et pour les perspectives ouvertes par leurs recherches qui m'ont guidé tout au long de cette étude. Merci au Pr Faleskog d'être venu depuis la Suède pour discuter de rupture ductile. Je remercie également Laurianne Pillon, avec qui j'aurai le plaisir de continuer à travailler sur l'endommagement ductile, pour la confiance qu'elle m'a accordée.

Pour m'accompagner et m'aider dans ce projet ambitieux, j'ai pu compter sur une riche équipe d'encadrement. Merci à Samuel pour son enthousiasme communicatif, sa pédagogie, pour m'avoir initié à la beauté et les nombreuses facettes de la mécanique, depuis ma première année aux Mines. Merci à Pierre de m'avoir fait découvrir les méthodes bayésiennes, et d'avoir toujours été disponible pour m'aider à interpréter leurs résultats. Merci à Sylvain et Victor pour avoir su se rendre disponibles afin de suivre mon travail, et d'avoir été mes relais à Valduc. Merci à Jacques, le maître des arcanes et des formulations exotiques de Z-set. Merci à Laurent pour sa débrouillardise avec Zset et Python, mais aussi pour nos passionantes discussions autour de sujets variés. Merci aussi à Franck pour l'analyse des images de tomographie.

Ce projet a été effectué pour le CEA Valduc, et je remercie l'ensemble de ma hiérarchie pour m'avoir permis de réaliser cette thèse dans d'excellentes conditions logistiques et expérimentales et de m'avoir toujours soutenu pour valoriser et communiquer mes recherches. Merci à tous les membres de mon laboratoire pour leur excellent accueil et leur aide expérimentale précieuse lors de mes (trop rares) visites à Valduc.

Au Centre des Matériaux, j'ai pu bénéficier d'un environnement de travail idéal. Merci en particulier à Basile, pour son aide avec le cluster, que j'ai bien fait chauffer pendant ces trois années. Merci à Henry pour son aide indispensable concernant les expériences de tomographie.

Merci à tous mes camarades futurs docteurs avec qui j'ai vécu l'aventure de la thèse. Je garde un souvenir ému de nos nombreuses discussions autour d'un café ou d'un déjeuner, ou pendant les conférences. Merci donc à

CONTENTS

Notations

Le soudage est un procédé courant permettant d'assembler de façon permanente des pièces métalliques par fusion locale et resolidification du matériau. Du fait du chargement thermique, les soudures ont des caractéristiques microstructurales complexes. En particulier, des pores (ou soufflures) peuvent être créés par piégeage de bulles de gaz lors de la solidification. Ces pores sont des défauts, pouvant réduire les performances mécaniques des pièces lors de leur utilisation. Nous nous intéresserons en particulier à leur résistance à la rupture ductile, mode de ruine fréquent associé à une forte déformation plastique au niveau local. Les pores peuvent être observés par des moyens de contrôle non destructif comme la tomographie X. Cependant l'utilisation directe des défauts observés dans une simulation mènerait à des coûts de calcul prohibitifs. L'objectif de l'étude est d'analyser l'influence des pores sur la rupture ductile de pièces soudées, et de proposer une méthode permettant de prendre en compte de façon simplifiée les images de tomographies afin de prévoir efficacement et précisément la résistance en rupture ductile.

L'influence des pores sur la rupture ductile peut être étudiée grâce à des simulations micromécaniques sur des cellules élastoplastiques contenant une distribution aléatoire de pores. De telles simulations permettent de prendre en compte l'interaction entre défauts. La détection de la rupture à l'échelle de la cellule, l'influence de l'orientation du chargement et l'effet de différentes populations seront en particulier étudiés.

A partir des simulations sur microstructures aléatoires, la dépendance de la déformation à localisation aux conditions de chargement peut être représentée par un métamodèle obtenu par krigeage. Une approche multifidélité est utilisée pour combiner les résultats de nombreuses simulations sur cellules à un pore, simulations rapides mais peu précises et un nombre plus limité de simulations plus coûteuses sur cellules aléatoires. Un critère d'endommagement ductile peut alors être déduit du métamodèle .

Enfin, la rupture ductile est étudiée expérimentatement par des essais in-situ observés par tomographie synchrotron. L'évolution des pores peut être suivie lors de leur déformation. Les résultats expérimentaux sont également comparés à ceux de simulations d'éprouvettes avec défauts. La comparaison est menée à la fois au niveau du comportement macroscopique, et à 10 CHAPTER 1. INTRODUCTION celui de la déformation locale des défauts.

Welding and its consequences on materials

Welding is a widespread assembly process used in numerous industries to permanently join metallic parts while ensuring material continuity [START_REF] Paillard | Procédés de soudage[END_REF]. Welding generally consists in locally and suddenly heating the base metal past its melting point so as to join the two parts when the metal solidifies back. Many variants of welding are available, differing by the origin of the energy needed to melt the metal.

For instance, in laser welding, the energy is provided by the absorption of a high power density laser beam [START_REF] Bruyère | A Phase Field Approach to Model Laser Power Control in Spot Laser Welding[END_REF]. The welding can be carried out with a series laser pulses of high power but low duration. This causes a sudden heating of the material, which melts and may even vaporize. A weld pool is created, but also a capillary, that is a zone left empty by the sublimated material. When the pulse ends, the material solidifies back and creates a bond between the two parts to be joined. To continue the weld joint, a new pulse must then be emitted on an area connected to the one already targeted.

The thermal cycling leads to a modification of the material characteristics. The metallographic section of welded TA6V plates in Figure 1.1 shows the extent of the transformation. Far from the weld, the metal remains in its initial state. Closer to the joint is the heat-affected zone, which has been heated during the pulse but has remained in a solid state. Finally at the heart of the joint is the fusion zone, where the material has been transformed into a liquid or even gaseous state before solidifying. Each zone exhibits a different microstructure. Thermal cycling may also induce residual stress in the materials. These metallurgical transformations were studied by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF]. Moreover, the fusion zone may contain weld pores, which are created if the weld capillary collapses and traps gas bubbles. [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF] characterized the population of such pores in a weld joint.

Ductile fracture of welded parts

Pores constitute defects that weaken the weld joint, due to the local lack of matter, and to stress concentration. In order to efficiently design parts, it is of relevance to engineers to determine to what extent the defects affect the mechanical resistance of welded parts.

The mechanical resistance depends on the applied loading conditions. For cyclic loading and therefore for fatigue failure mode, defects are known to reduce lifetime, especially if they are close to the surface. The influence of the weld pores on the lifetime was studied extensively by [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF].

Another loading condition of interest consists in applying a monotonic loading on the part until it fails. At room temperature, numerous metals (including TA6V alloy, and even in weld joints [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF]) fail by ductile fracture. This type of fracture is characterized by significant plasticity, at least at small scales. The typical failure pattern shows several stages (see Figure 1.2). Voids are first nucleated within the metal, for instance by decohesion near inclusions or between phases. Increasing the applied displacement then leads to the growth of the voids. The final stage of failure corresponds to voids coalescing, i.e. merging together either by internal necking, or by the nucleation of many new small secondary voids between existing voids. Depending on the loading conditions and the materials, the extent of each of the previous phases may change. Voids facilitate ductile fracture as stress concentrators favoring plasticity and nucleation in their vicinity, and by reducing effective cross-section. However, unlike fatigue, nucleation would naturally produce voids, even if pores were not initially present.

Non destructive controls and justification of exemptions

As gas pores are defects, process parameters should be chosen to limit their number and size. However, even with an optimized process, pores may still arise in weld joints. Under these conditions, it is important to verify that the produced parts meet the mechanical strength requirements despite the presence of defects. In the case of critical or high value-added parts manufactured in small batches, it is neither viable nor sufficient to sacrifice a part of the production for destructive mechanical tests so as to deduce the average properties of the product. Instead, each part produced must be non-destructively tested.

Non destructive controls allow to detect pores. Surface defects can be easily detected by visual inspection, dye penetrant inspection, magnetoscopy or even analysis of eddy currents [START_REF] Cherfaoui | Essais non destructifs[END_REF]. However defects deeper within the weld joint require volume methods such as ultrasonic testing or X-ray radioscopy. This study will focus on X-ray methods. In X-ray radioscopy, an image is produced by illuminating the tested specimen with an X-ray beam, producing images where voids appear as zones of lower X-ray absorption. By itself, radioscopy only yields two-dimensional information and not the precise location of voids in space.

By combining several images of the same specimen taken under different angles (for instance by rotating the specimen), a three-dimensional image can be reconstructed by Radon transform. This process is known as X-ray Computed Tomography (CT). Pores are then precisely located in three dimensions within the weld joint. The resolution and thus ability to detect defects naturally depends on the apparatus used. Whereas laboratory tomographs, which are increasingly available, allow to detect defects of the order of 50 µm, higher resolution values, better than 1 µm, can be reached with synchrotron beam lines. [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF] identified with synchrotron CT the pore population in a weld joint (Figure 1.3). Figure 1.3: Pore population in a weld joint, observed with synchrotron tomography by [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF] The defect population obtained with tomography does not directly indicate the mechanical resistance of the specimen. A first approach is to only use very simple geometric information like the total porosity or the diameter of the largest pore, as in the standard NF EN ISO 13919-1 (2019). If such quantities are larger than some threshold, the weld joint is deemed insufficiently resistant. This approach is simple to use, but makes limited use of the physical understanding of failure. Therefore, it is generally too conservative and may lead to discard parts which would satisfy otherwise all the requirements of the industrial application. Another approach is to use numerical simulations to predict the part's mechanical behavior, and include as much information as possible from the tomography images. For instance, in finite element simulations, all the visible defects could be meshed and the specimen could then be simulated until its failure. Provided that a suitable model is used for the matrix, this approach would lead to a precise estimation of mechanical resistance to monotonic loading. With this precise estimation, parts not satisfying the standards on geometric information, but with enough predicted resistance, could be sensibly accepted. This justification of exemption allows reducing process costs and delays associated with discarding functional parts.

However the simulations precisely accounting for all information would require a high number of degrees of freedom and be untractable. At least they would require too long computational power to be used in a production process and to respect production schedule. How should therefore tomography information be taken into account for a precise estimation of resistance to ductile loading without prohibitive computational times?

Objectives and methodology of the study

The main objective of the present study is therefore to investigate the role of weld pores on the ductile fracture behavior, and to propose a method to account for the presence of these defects in a simplified way. The study will focus on a welding process of industrial interest for the CEA. More precisely, plates of titanium alloy TA6V are joined with pulsed laser welding. This process was already studied by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF] (metallurgy and ductile fracture from a macroscopic point of view) and [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF] (pore population and their behavior in fatigue). However this study could be of interest for the ductile fracture of any porous material, including additively manufactured materials.

More generally, the key research problem is to understand how a population of preexisting pores influences ductile fracture. This issue can be investigated in a simplified and model setting: a small volume of material containing a population of pores is subjected to monotonic loading. This micromechanical approach captures elementary ductile failure mechanisms and may be representative of the local behavior of the material in a structure.

The micromechanical approach used initially purely mathematical methods and allowed creating analytic ductile failure criteria such as [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] well-known model. However, such analytic methods require simplifying hypotheses and may not represent the whole complexity of the ductile fracture mechanisms. In a computational homogenization approach, the behavior of the volume element may also be simulated, by finite elements for instance. However most existing studies only considered cells with a single void. Such unit cells cannot represent correctly the interaction between different voids. This is especially an issue for high porosity microstructures. Simulations of random cells, i.e volume elements with randomly dispersed voids, may therefore help assess the interaction of the defects. [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] obtained yield surfaces from simulations on such cells, but this only corresponded to the early stage of ductility. In order to understand ductile fracture, simulations should reach the coalescence or final failure of the random cells.

The present study therefore deals with random microstructures composed of an elastoplastic matrix surrounding a population of randomly dispersed voids. They are simulated until ductile failure for various loading conditions. The effect of the loading and simulation conditions are first assessed to verify the soundness of the micromechanical approach. Proportional loading conditions are chosen. Ductile failure is known to depend, for isotropic materials, mainly on stress invariants, such as the stress triaxiality and the Lode parameter. However the loading orientation with respect to the (cubic) cell may also play a role, and this model-induced anisotropy is carefully studied. Moreover, detecting the onset of failure in a cell simulation is not straightforward: There is no damage in the matrix material, and the ductile fracture is then a purely cell-level phenomenon. Two failure criteria are used and compared: a geometric coalescence criterion and a more mathematically-justified localization criterion (the formation of strain localization bands is a precursor to ductile fracture). The influence of void interaction is studied by comparing several void populations. For simplicity, populations of identical spherical defects are used. However the dispersion of failure-related quantities, and the role of the number of simulated voids are assessed.

Cell simulations allow to understand the influence of voids on ductile fracture. In order to simulate ductile fracture in real structures and weld joints, a hybrid approach can be used. Whereas the largest visible voids may be explicitly represented (i.e meshed) in a finite element simulation, the smaller voids have to be considered implicitly, for instance with a ductile failure model. This model can be constructed from the simulations on random microstructures.

So as to construct the model, the evolution of strain localization with respect to loading conditions is required. A metamodelling strategy is used to construct a response surface from a set of simulations at various loading conditions, so as to allow non parametric models. However, simulations on random microstructures are computationally expensive, much more than simulations on unit cells. A multifidelity approach is therefore used to combine precise, but expensive results on random cells, and results on unit cells, which only give the overall trends but require less computational budget. This metamodelling approach uses kriging [START_REF] Matheron | Traité de géostatistique appliquée[END_REF], also known as Gaussian process regression. A key requirement is to compensate the effect of the loading orientation on random cell simulations, for it is not related to physical mechanisms but due to biased model simulation conditions.

Finally the weld joint may be studied from a more experimental point of view. The macroscopic mechanical behavior was already studied by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF] who proposed a ductile failure model. However, how do the defects evolve during deformation and influence ductile failure? In situ tensile tests are thus carried out and imaged by synchrotron tomography. Unlike [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF], the focus is not on the initial population, but on observing the deformation of individual defects throughout the tensile experiment. The tensile tests are also simulated, so as to compare with experimental results. The largest defects are explicitly taken into account. Moreover a ductile damage criterion is derived from the macroscopic response surface on random microstructures and implemented in the FE software. With the full simulation of the specimens with defects, simulations can be compared to experiments both for the macroscopic behavior and the microscopic deformation of defects.

Organization of the dissertation

The dissertation is organized as follows:

• Chapter 2 describes the methodology for creating random microstructures made of an elastoplastic matrix with a population of random dispersed voids, and then simulating its behavior until failure. Simple loading conditions are applied and cell failure is detected with a geometric coalescence criterion. The dispersion of results with regard to the realization of random microstructure is quantified, and the role of the number of voids is assessed.

• Chapter 3 extends the preceding methodology by considering arbitrary loading orientations.

A more powerful and mathematically-justified localization criterion (Rice criterion) is introduced and implemented. The anisotropy due to loading orientation is carefully studied.

• Chapter 4 uses the results of cell simulations to construct a (non-parametric) response surface for strain localization. Results from unit and random cells are combined with a multifidelity Gaussian Process Regression approach.

• Chapter 5 describes in situ tensile tests experiments on TA6V weld joints, observed by synchrotron tomography. The tests are also simulated, allowing to compare experimental and numerical results both for the overall mechanical behavior and the local evolution of defects.

Chapter 2

Ductile Fracture of Materials with Randomly Distributed Voids

A reliable determination of the onset of void coalescence is critical to the modelling of ductile fracture. Numerical models have been developed but rely mostly on analyses on single defect cells, thus underestimating the interaction between voids. This study aims to provide the first extensive analysis of the response of microstructures with random distributions of voids to various loading conditions and to characterize the dispersion of the results as a consequence of the randomness of the void distribution. Cells embedding a random distribution of identical spherical voids are generated within an elastoplastic matrix and subjected to a macroscopic loading with constant stress triaxiality and Lode parameter under periodic boundary conditions in finite element simulations. The failure of the cell is determined by a new indicator based on the loss of full rankedness on the average deformation gradient rate. It is shown that the strain field developing in random microstructures and the one in unit cells feature different dependencies on the Lode parameter L owing to different failure modes. Depending on L, the cell may fail in extension (coalescence) or in shear. Moreover the random void populations lead to a significant dispersion of failure strain, which is present even in simulations with high numbers of voids.

This chapter is mainly based on the homonymous article, published in the International Journal of Fracture [START_REF] Cadet | Ductile fracture of materials with randomly distributed voids[END_REF]. A new appendix dedicated to the position of localization bands was added.

Résumé détaillé

La validation des modèles d'endommagement ductile, et plus généralement l'étude des phénomènes mécaniques qui y sont associés, passe fréquemment par la simulation numérique de la rupture de volumes élémentaires. Le plus souvent, ces volumes élémentaires sont constitués d'une matrice élastoplastique contenant un seul défaut (vide). Cependant, une telle cellule unitaire risque de représenter de façon trop simplifiée l'interaction possible entre les défauts. La ruine de cellules avec plusieurs défauts a cependant été peu étudiée. Ce chapitre se propose donc de simuler la rupture ductile de microstructures constituées d'une distribution aléatoire de défauts.

Nous présentons d'abord la méthodologie pour construire et mailler de telles microstructures, constituées d'une matrice élastique/parfaitement plastique incluant une population de défauts répartis aléatoirement. Pour la simulation par éléments finis, des conditions aux limites périodiques sont imposées. De plus, un chargement proportionnel est appliqué de telle manière que l'état de contrainte, caractérisé par la triaxialité des contraintes T et le paramètre de Lode L, reste constant tout au long de la simulation.

Identifier l'instant exact auquel la ruine de la cellule se produit est une tâche complexe. Cet instant peut être décrit par la survenue des phénomènes de coalescence (liée à la fusion de défauts proches, mais pouvant s'interpréter comme un état de déformation particulier) ou de localisation (liée à une perte d'ellipticité des équations mécaniques). Ces deux critères sont susceptibles de fournir des informations complémentaires. Un critère en coalescence, construit à partir du déterminant de la vitesse en déformation moyenne, sera utilisé ici (le prochain chapitre traitera de la localisation).

Les comportements en rupture ductile des cellules unitaires et des microstructures aléatoires peuvent alors être comparés. Lorsque la coalescence se produit, des bandes de concentration de la déformation plastique sont visibles dans les deux types de cellules. Ces bandes sont principalement planes, même si dans le cas des cellules aléatoires, les bandes peuvent serpenter pour joindre différents défauts. L'orientation de ces bandes est dictée par les conditions aux limites : orthogonale à l'axe de chargement principal (mode d'extension) ou à 45 • de celui-ci (mode de cisaillement). Par ailleurs, l'évolution de la déformation à rupture avec T est similaire pour les deux types de cellules, mais très différente quant à l'évolution en fonction de L. Pour un chargement en cisaillement généralisé (L ∼ 0), la déformation à rupture est fortement abaissée pour les microstructures aléatoires. Ceci est dû à un passage du mode de rupture en cisaillement plutôt qu'en extension. Ce changement est beaucoup plus restraint dans le cas de la cellule unitaire.

Par ailleurs, la dispersion due aux microstructures aléatoires peut être étudiée. Cette dispersion est élevée: il peut y avoir un facteur 4 sur la déformation à rupture pour deux réalisations de microstructures aléatoires soumises à un même chargement. Cette dispersion décroît avec le nombre de défauts mais reste élevée même pour de grandes microstructures avec une centaine de défauts.

Enfin, plusieurs aspects de la méthode de simulation des microstructures aléatoires sont discutés. Les résultats du critère en coalescence sont comparés à ceux d'un autre critère de la littérature, proposant une interprétation énergétique de la coalescence. Les deux critères fournissent des résultats similaires pour une rupture en extension, mais notre critère semble mieux réprésenter la rupture en mode de cisaillement. Par ailleurs, l'effet du comportement du matériau, et notamment de l'écrouissage est discuté: avec certains écrouissages, la différence entre cellules unitaire et aléatoires est réduite.

Introduction

Predicting the failure of a structural part subjected to monotonic loading requires a good understanding of the ductile fracture behavior of the material. However ductile fracture is a complex phenomenon involving a variety of mechanisms, strongly dependent on the material and involving large strain at least on a local level (Besson, 2004). Voids are first nucleated within the material, especially near second phase inclusions. Depending on the loading conditions, the voids may or may not grow. Finally the material fails when voids coalesce, either by internal necking or by the nucleation of secondary voids (mostly for shear-dominated loading). Softening due to void growth may also be sufficient to initiate failure without coalescence per se [START_REF] Tekoglu | On localization and void coalescence as a precursor to ductile fracture[END_REF]. Although a large body of literature on ductile fracture has already been published, accurate failure prediction is still a research problem, as evidenced by the Sandia fracture challenges [START_REF] Boyce | The Sandia Fracture Challenge: blind round robin predictions of ductile tearing[END_REF][START_REF] Boyce | The second Sandia Fracture Challenge: predictions of ductile failure under quasi-static and moderate-rate dynamic loading[END_REF][START_REF] Kramer | The third Sandia fracture challenge: predictions of ductile fracture in additively manufactured metal[END_REF]: besides the difficulty of calibrating modelling parameters from experimental data, predicting ductile failure requires to take into account many strongly nonlinear physical processes.

Experimental studies have shown that the failure behavior strongly depends on the stress state to which the material is subjected. The effects of the stress triaxiality (ratio of the von Mises equivalent stress to the mean stress) and the Lode parameters (reflecting the third stress invariant) have been extensively investigated (for instance by [START_REF] Helbert | The influence of stress triaxiality on the damage mechanisms in an equiaxed Ti-6AI-4V alloy[END_REF], [START_REF] Bao | On fracture locus in the equivalent strain and stress triaxiality space[END_REF], [START_REF] Barsoum | Rupture mechanisms in combined tension and shear-Micromechanics[END_REF], [START_REF] Gao | Effects of the stress state on plasticity and ductile failure of an aluminum 5083 alloy[END_REF], [START_REF] Dunand | On the predictive capabilities of the shear modified Gurson and the modified Mohr-Coulomb fracture models over a wide range of stress triaxialities and Lode angles[END_REF], [START_REF] Gilioli | Effect of Triaxiality and Lode Angle on the Plasticity Behaviour of a Ti-6Al-4V Titanium Alloy[END_REF], [START_REF] Zhai | Modeling the ductile damage process in commercially pure titanium[END_REF], [START_REF] Xiao | Effect of the Lode parameter in predicting shear cracking of 2024-T351 aluminum alloy Taylor rods[END_REF], and [START_REF] Zhang | Ductile fracture prediction using enhanced CDM model with Lode angle-dependency for titanium alloy Ti-6Al-4V at room temperature[END_REF]). Models representing ductile failure should therefore account for the effect of these two parameters.

Analytic and computational approaches at a micromechanical level have also been developed to investigate the mechanisms of ductile fracture, to model ductile fracture and provide microscaleinformed failure prediction. Following [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] results from limit analysis, increasingly precise analytic models have been developed by explicitly representing approximate strain fields near voids in a plastic material. [START_REF] Besson | Continuum models of ductile fracture : A review[END_REF] provides a review of such models but more recent ones have been developed to represent void growth and coalescence either by necking or in shear [START_REF] Benzerga | Effective Yield Criterion Accounting for Microvoid Coalescence[END_REF][START_REF] Morin | A unified criterion for the growth and coalescence of microvoids[END_REF][START_REF] Torki | A unified criterion for void growth and coalescence under combined tension and shear[END_REF][START_REF] Nguyen | A nonlocal approach of ductile failure incorporating void growth, internal necking, and shear dominated coalescence mechanisms[END_REF], and can be used for practical applications [START_REF] Keralavarma | Ductile failure as a constitutive instability in porous plastic solids[END_REF]. The limit analysis approach was also extended by [START_REF] Leblond | A theoretical approach of strain localization within thin planar bands in porous ductile materials[END_REF] to random distribution of voids. Computational studies of ductile fracture may help validate these models or provide valuable insights in the failure mechanisms on their own, for example by quantifying the effect of stress triaxiality and Lode parameter [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF][START_REF] Zhu | Combined effects of triaxiality, Lode parameter and shear stress on void growth and coalescence[END_REF] or by distinguishing strain localization from coalescence [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF][START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF][START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF]. However these studies are mostly carried out on unit cells: the global behavior of the material is summarized by that of a meshed cell containing a single void. Even though this approach was proven useful to analyze fundamental mechanisms at the void level at a low computational cost, it oversimplifies the interaction between voids, whose influence increases with porosity, by assuming that voids are regularly distributed as a cubic lattice.

Some studies have investigated the interaction of voids in simplified configurations, involving only a couple of voids. For instance [START_REF] Bandstra | On the influence of void clusters on void growth and coalescence during ductile fracture[END_REF] considered three-voids clusters with rotational symmetry in a hexagonal volume element; [START_REF] Tvergaard | Effect of void cluster on ductile failure evolution[END_REF] and [START_REF] Tvergaard | Nucleation from a cluster of inclusions, leading to void coalescense[END_REF] considered 2D clusters with three aligned pores, whereas Trejo [START_REF] Navas | Void growth and coalescence in a three-dimensional non-periodic void cluster[END_REF] systematically studied 3D three pore clusters. [START_REF] Khan | On the role of secondary voids and their distribution in the mechanism of void growth and coalescence in porous plastic solids[END_REF] investigated the interaction between two populations of voids, in the simplified context of a high symmetry periodic arrangement. However, in a real material, a large number of voids, with complex spatial distribution interact with each other. [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF] considered 2D microstructures with a random population of voids and showed that increased triaxiality accelerates coalescence. [START_REF] Shakoor | Ductile fracture of a metal matrix composite studied using 3D numerical modeling of void nucleation and coalescence[END_REF] also provided a very detailed description of the mechanisms of ductile fracture from nucleation up to coalescence, between randomly distributed voids. All these studies evidence the role of clusters but do not allow to compute coalescence properties depending on loading conditions, as a model of ductile fracture would require, because they investigate too few void configurations and loading cases.

Analytical approaches can take random void distributions into account. For instance, [START_REF] Leblond | A theoretical approach of strain localization within thin planar bands in porous ductile materials[END_REF] developed a limit analysis model coupling coalescence and shear band formation initially for a periodic distribution, but proposed a method to extend it to the random case by considering all possible orientations of the shear bands. Moreover, works by [START_REF] Danas | A finite-strain model for anisotropic viscoplastic porous media: I -Theory[END_REF] or [START_REF] Vincent | Porous materials with two populations of voids under internal pressure: I. Instantaneous constitutive relations[END_REF] for instance, considered random void populations within the context of a nonlinear variational homogenization scheme: the porous medium was compared to a linear composite, whose stiffness is based on Ponte [START_REF] Castañeda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF] bounds, an effective method to represent a population of random elliptical voids. This variational technique was subsequently used by [START_REF] Danas | Influence of the Lode parameter and the stress triaxiality on the failure of elasto-plastic porous materials[END_REF] to investigate the influence of stress triaxiality and Lode parameter. However, such analytical approaches should be compared to simulations to check the validity of their assumptions. For instance, [START_REF] Danas | Influence of the Lode parameter and the stress triaxiality on the failure of elasto-plastic porous materials[END_REF] predictions for the behavior at low triaxiality were found to be unrealistic by [START_REF] Hutchinson | Comment on "Influence of the Lode parameter and the stress triaxiality on the failure of elasto-plastic porous materials" by K. Danas and P. Ponte Castañeda[END_REF] FEM simulations on unit cells with the same loading conditions. Explicit simulations of random void distributions have been carried out in a limited number of works. [START_REF] Bilger | Effect of a nonuniform distribution of voids on the plastic response of voided materials: a computational and statistical analysis[END_REF] and [START_REF] Bilger | Bounds and estimates for the effective yield surface of porous media with a uniform or a nonuniform distribution of voids[END_REF] using Fast Fourier Transform then [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] with Finite Elements Analysis proposed a computational homogenization method to determine an effective yield surface. Several microstructures consisting of a random void distribution embedded in a plastic matrix are simulated up to overall plastic yield for several loading conditions. The results are averaged over the several microstructures to determine a homogenized yield surface (represented for [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] by a GTN criterion). This approach was extended to a Green-type porous matrix [START_REF] Fritzen | Computational homogenization of porous materials of Green type[END_REF], to multiple void populations of different size [START_REF] Khdir | Computational homogenization of plastic porous media with two populations of voids[END_REF] and to non-spherical voids [START_REF] Khdir | A computational homogenization of random porous media: Effect of void shape and void content on the overall yield surface[END_REF]. However these studies were focused on yield surface and did not address coalescence. Recently, Hure (2021) did perform FFT simulations on cells with multiple voids up to coalescence, and illustrated the influence of the number of voids on the stress at coalescence. Yet this study was limited to the simple case of axisymmetric loading.

To the authors' knowledge, a description of coalescence for various loading conditions and at the level of a representative volume element with multiple voids, has not been done yet. We therefore propose here to extend the methodology of [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] and [START_REF] Hure | Yield criterion and finite strain behavior of random porous isotropic materials[END_REF] to the study of coalescence under various stress states. We aim to assess the effect of the interaction between randomly distributed voids on the macroscopic failure response of a cell, depending of the stress state. The results should be compared to those of unit cells to identify how they differ from cells with multiple voids. Moreover the statistical dispersion in failure results linked to the random distribution should be quantified.

To this end, cells composed of a random population of identical spherical voids are generated and subjected to various loading conditions, characterized by constant stress triaxiality and Lode parameter levels, in finite element simulations with Z-set software [START_REF] Besson | Object-Oriented Programming Applied to the Finite Element Method Part I. General Concepts[END_REF]. We then propose a coalescence indicator based on the loss of full rankedness of the macroscopic deformation gradient rate. The identification of coalescence during the simulation allows to extract several quantities of interest at the onset of coalescence. Our main results show that the evolution of the onset of coalescence with respect to the Lode parameter is qualitatively different between random microstructures and unit cells. This difference is associated with a change of coalescence modes for random microstructures. Finally, dispersion of the results due to the randomness of the void distribution is studied.

Section 2.3 describes the methodology used to generate random microstructures and to prescribe the loading conditions within the FE simulation. In section 2.4 typical numerical simulation results are presented and an indicator is defined to identify failure. Section 2.5 applies the methodology of sections 2.3 and 2.4 to compare the response of random microstructures to that of a unit cell, both on the evolution of macroscopic (cell-level) quantities, and on plastic strain field patterns. The dispersion of the results is also investigated. Finally, we discuss in section 2.6 the simulation hypotheses chosen in this work, and verify to what extent the results can be generalized.

An intrinsic notation is used for tensors: vectors, as first order tensors, are represented as v = v i e i and second order tensors as A ∼ = A i j e i ⊗ e j , where (e i ) is an orthonormal frame. The subscript 0 in the notation A 0 refers to the value of A in the initial configuration at time t = 0. The position of a material point initially at x 0 evolves with time t as x = Φ (x 0 ,t); the deformation gradient is then defined as

F ∼ = ∂ Φ ∂ x 0 .
Quantities decorated with an overlying bar, such as Ā, refer to the macroscopic counterpart (at the level of a cell) of a quantity A defined locally. For instance F ∼ is the average deformation gradient (defined more precisely in section 2.3.3), and J = det F ∼ .

Methodology

Generation of random microstructures and finite element meshing

The methodology to create the elementary volumes follows that of [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF]. These cells consist of a cubic matrix containing a population of identical non overlapping spherical defects. As all the N de f ects spheres have the same radius r, the porosity of the cell (of size L cube and therefore of volume V 0 = L 3 cube ) is defined as:

f = 4π 3 N de f ects r L cube 3 (2.1)
The radius of the voids is fully determined once the porosity and the number of voids are chosen. The initial porosity was chosen as f 0 = 6%, to be compared to the range of porosity levels f 0 ∈ [0.1%, 30%] considered by [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF]. However unit cell analyses frequently study lower porosities, with f 0 ∼ 0.1% [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF][START_REF] Vishwakarma | Micromechanical modeling and simulation of the loading path dependence of ductile failure by void growth to coalescence[END_REF][START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF]. For low porosity values, interactions between defects can indeed be neglected [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF], at least for the growth phase. [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] showed for instance that unit cells and random microstructures with sufficiently low prosity levels have a similar growth behavior, which can be represented by a GTN criterion.

Nonetheless, high porosity levels of 6% are possible in sintered materials [START_REF] Becker | The effect of porosity distribution on ductile failure[END_REF], nodular cast iron [START_REF] Zhang | Ductile fracture of materials with high void volume fraction[END_REF], irradiated stainless steel [START_REF] Cawthorne | Voids in Irradiated Stainless Steel[END_REF]. Moreover, overall porosity of 0.5% to 2% can be found in weld joints [START_REF] Li | Influence of mechanical mismatching on the failure of welded joints by void nucleation and coalescence[END_REF][START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF][START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF], but porosity values defined at a smaller scale, near void clusters, can be higher. A high initial porosity level can also provide insight for coalescence at lower initial porosity levels. Coalescence starts after a sufficient phase of growth so that voids begin interacting with each other and can no longer be considered isolated, which means that the porosity is no longer negligible. Notwithstanding the change of void shape which will play a significant role, starting at high porosity is equivalent, to some extent, to considering the end of a simulation at lower porosity.

The position of the defects is chosen according to a Poisson sphere process [START_REF] Matern | Spatial Variation[END_REF]. As the target porosity is significantly lower than the jamming porosity levels that characterize such processes (around 38 % according to [START_REF] Gamito | Accurate Multidimensional Poisson-Disk Sampling[END_REF]), a dart-throwing method is sufficient for the sampling. The position of the center of a sphere is chosen according to a uniform distribution on the cube. If the distance between the resulting spheres and the already built defects is larger than 10% of the radius of a sphere, the new sphere is included in the list of defects. Otherwise it is rejected and a new possible center is chosen randomly. Introducing a repulsion distance between the defects allows a better mesh quality. During the FEM simulations, periodic boundary conditions are applied (see section 2.3.3). Therefore a periodic microstructure and in turn a periodic mesh should be used. In order to ensure the periodicity of the population of defects, each time a new defect intersects a side of the cube, it is copied on the other side (thus there are four copies if an edge is intersected, and eight if the defect contains a vertex of the cube). All of these copies are taken into account to determine intersections between defects. [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] verified several statistical properties of the representativeness of this process.

The cell with the preceding defect population is meshed with NETGEN software [START_REF] Schöberl | NETGEN -An advancing front 2D/3D-mesh generator based on abstract rules[END_REF]. This tool first meshes surfaces, then volumes, and generates a non structured tetrahedral mesh. Periodicity of the mesh is imposed so that opposite sides of the cube have identical surface meshes. A maximum element size of h cell = r is imposed globally on the cell, but on the surface meshes of the defects the maximum element size is reduced to h void = r/5. The mesh is thus refined on the part of the surface mesh corresponding to the surface of voids. Finally, tetrahedral second-order 10-node elements with reduced integration are used to limit volume locking (due to large strain plasticity) in the FEA simulations. An example of the meshing of a microstructure with 27 voids is shown in figure 2.1a.

As cells are cubic, they define a canonical orthonormal frame (O, e 1 , e 2 , e 3 ) where O is a vertex of the cube and the unit vectors e 1 , e 2 , e 3 are parallel to edges of the cube (in the initial configuration). All tensor components will be expressed in this frame.

Although a diversity of microstructures were used in this study, several are repeatedly referred to in this article; they are shown in figure 2.1. The microstructures R1 and R2 are two random microstructures with 27-voids of radius r = 0.08L cube, R1 . The one-pore cell unit is defined as a cubic matrix of size L cube, unit = L cube, R1 /3 containing a unique defect of radius r = 0.08L cube, R1 (same radius as before, and thus same volume fraction). It is meshed with the same procedure and same parameters as the larger cells with a random population. Finally the microstructure lattice consists in 3 × 3 × 3 defects arranged on a cubic lattice; it is meshed in the same way as the random microstructures, so the mesh is not the assemblage of 27 small identical meshes of the unit cells.

Material behavior law at finite strain

Finite element simulations are carried out using Zset software [START_REF] Besson | Object-Oriented Programming Applied to the Finite Element Method Part I. General Concepts[END_REF]Zset, 2022). As the matrix can undergo large deformation before coalescence, the simulations are performed in a finite strain framework. A local objective frame approach is adopted to formulate the constitutive law of the matrix [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]. The strain rate D ∼ and Cauchy stress σ ∼ tensors are convected in a corotational frame:

ė ∼ = Q ∼ T D ∼ Q ∼ (2.2) s ∼ = JQ ∼ T σ ∼ Q ∼ (2.3) where Q ∼ is a rotation matrix verifying -Q ∼ T Q ∼ = Q ∼ T Q ∼ = W ∼ (W ∼ = skew( Ḟ ∼ F ∼ -1
) is the material spin tensor). This choice of corotational frame is equivalent to using the Jaumann derivative of the stress in the hypo-elasticity law.

The constitutive law is then defined by a classical additive decomposition of convected strain rates in an isotropic elastic part and a plastic part. Isotropy and time-independent perfect plasticity (absence of hardening) with a von Mises yield criterion are assumed for the matrix material:

ė ∼ = ė ∼ e + ė ∼ p e ∼ e = 1 + ν E s ∼ - ν E (tr s ∼ )1 ∼ s vm = 3 2 s ∼ dev : s ∼ dev ψ(s ∼ ) = s vm -R 0 ≤ 0 ė ∼ p = ṗ ∂ ψ ∂ s ∼ (2.4) x y z
(a) Random microstructure R1 with 27 voids (FE mesh with 255628 nodes).

x y z (b) Random microstructure R2 with 27 voids (FE mesh with 176982 nodes).

x y z (c) Unit cell with one void. stress and ṗ = 2 3 ė ∼ p : ė ∼ p playing the role of the plastic multiplier. The Young modulus, the Poisson ratio and the yield strength are respectively chosen as E = 200 GPa, ν = 0.3 and R 0 = 500 MPa, hence R 0 /E = 0.0025. The cumulative plastic strain is defined as:

p = t 0 ṗ dt (2.5)
where t is actually a fictitious time in rate-independent plasticity, acting as an increasing loading parameter.

During the finite element analysis, this constitutive law is integrated at each quadrature point of the finite element mesh by an implicit Euler method, then the global static equilibrium is solved in total Lagrangian formulation by a Newton-Raphson scheme with a consistent tangent matrix.

Boundary and loading conditions

The boundary and loading conditions follow that of [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]. Periodic boundary conditions are applied on the sides of the cube [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]. The displacement field u should therefore have the form:

u = ( F ∼ -1 ∼ ) • x 0 + v (x 0 ) (2.6)
with F ∼ the average deformation gradient, and v a displacement fluctuation field, periodic and with zero average over the cell. The periodicity of v and the anti-periodicity of traction vectors mean that:

v (x + 0 ) = v (x - 0 ) (2.7) σ ∼ • n (x + 0 ) = -σ ∼ • n (x - 0 ) (2.8)
if x + 0 and x - 0 represent two homologous points on opposite sides of the periodic mesh and n (x 0 ) represents the outward-pointing normal to the mesh boundary at x 0 . In this formulation, the degrees of freedom are the three components of the displacement fluctuation field for each node of the mesh and the nine components of F ∼ (or rather of E ∼ = F ∼ -1 ∼ ). The macroscopic Boussinesq (or first Piola-Kirchhoff) and Cauchy stress tensors are defined by:

S ∼ = 1 V 0 V 0 S ∼ dV 0 = 1 -f 0 V 0 V matrix 0 S ∼ dV 0 (2.9) σ ∼ = 1 J S ∼ . F ∼ T
(2.10)

where J = det F ∼ and V 0 is the volume of the cell (matrix and defects) in the initial configuration. The integration on V 0 implicitly considers that stress is well-defined and identically zero in the voids.

The simulations are carried out at constant (macroscopic) stress triaxiality and Lode parameter. These quantities are here defined as:

T = tr σ ∼ 3 σvm (2.11) L = 2 σ2 -σ1 -σ3 σ1 -σ3 (2.12)
where σeq is the von Mises equivalent stress calculated from σ ∼ and σ1 ≥ σ2 ≥ σ3 (with σ1 > σ3 ) are the three eigenvalues of σ ∼ . L = -1, L = 0 and L = 1 respectively correspond to states of generalized tension, shear and compression. An alternative definition of a Lode parameter with L = 1 for tension and L = -1 for compression can also be found in literature (e.g. [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF] and [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF]).

To ensure that T and L remain constant during the simulation, a special macroscopic spring element was developed by [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]. It acts on the E i j degrees of freedom, and its reaction forces are calculated so that σ ∼ keeps the following diagonal form throughout the simulation:

σ ∼ =   σ1 0 0 0 σ2 0 0 0 σ3   = σ11   1 0 0 0 η 2 0 0 0 η 3   (2.13)
where η 2 = σ2 / σ1 and η 3 = σ3 / σ1 are prescribed constants which define the stress state. Therefore the eigenvectors of σ ∼ are collinear to the three axes of the cube. Unlike [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], [START_REF] Dunand | Effect of Lode parameter on plastic flow localization after proportional loading at low stress triaxialities[END_REF], [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF], and [START_REF] Zhu | Combined effects of triaxiality, Lode parameter and shear stress on void growth and coalescence[END_REF] but like [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF] and [START_REF] Guo | Crystal plasticity based investigation of micro-void evolution under multi-axial loading conditions[END_REF], we chose not to consider the effect of a shear stress component (for instance σ12 ) for computational cost reasons. However the cubic cell exhibits cubic symmetry and has an anisotropic behavior. The additional stress component could allow different loading orientations with identical T and L values. The consequences of this choice will be discussed in section 2.6.1.

To prevent degeneracy of solutions due to rigid body motion, a global translation and a global rotation of the cube should be fixed. The translation is taken care of by fixing a vertex of the cube. For the rotation, a possible method is to impose three additional constraints on the average deformation gradient F ∼ . For instance F ∼ can be supposed symmetric, as done by [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]:

F12 = F21 F23 = F32 F13 = F31 (2.14)
Another method is presented and discussed in appendix 2.B.2. With the aforementioned conditions, the simulation can be strain-controlled by specifying only the average strain along the first axis E 11 = F11 -1. Strictly speaking, the tensor E is not a strain tensor, but it will be used all the same to control the simulation and measure deformation.

We impose E 11 = εt, with ε an arbitrary strain rate (the value can be arbitrarily chosen, as the plasticity is time-independent). At the beginning of the simulation t = 0, the cell is undeformed, and F ∼ = 1 ∼ .

Description of a coalescence indicator

Typical results of a computation

With the method described in the previous subsection, simulations can be carried out for several loading conditions. In this study, we are interested in the evolution of several quantities at failure. However defining ductile failure and detecting it during the simulation is not straightforward. To illustrate this issue, some enlightening simulation results will be described first. Three simulations on the microstructure R1 were carried out in generalized tension (L = -1) at three triaxiality levels T = 0.8, T = 1 and T = 1.4. The figure 2.2a compares the macroscopic Cauchy and Boussinesq stress components along the main loading axis (the marker on the curves corresponds to the indicator described later in section 2.4.3). The three loading conditions lead to a similar evolution of stress. The stress maximum is reached shortly after the beginning of the computation (for E 11 < 0.01) then the stress decreases monotonically and almost linearly. However, at a critical strain that depends on the loading condition, the decrease in stress suddenly accelerates and the unit cell quickly loses all its load-bearing capacity (at approximately E 11 = 0.5, 0.35, 0.18 for T = 0.8, 1, 1.4 respectively). This event can be thought as the failure of the cell. Moreover, at the same strain as the onset of stress drop, the transverse strain E 22 stabilizes (fig. 2.2b). This macroscopic failure is also related to the behavior at a more microscopic level. Figure 2.2 shows the cumulative strain field p shortly after this failure, for T = 1: plastic strain is concentrated in a band, mostly parallel to a side of the cube (and perpendicular to the main loading axis) but its exact shape fits closely the distribution of voids.

Although the stress decrease acceleration is clearly visible on the stress-strain plots in generalized tension, it is difficult to define its exact location so as to determine the precise failure onset and compute relevant physical quantities at this instant. Moreover the stress decrease does not generalize to shear-dominated loading conditions. Therefore a more precise failure indicator remains to be determined.

Available failure indicators in the literature

Several criteria for ductile failure in a unit cell have been developed, and are reviewed for instance by [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF]. The earliest approaches were purely geometrical: [START_REF] Brown | Initiation and growth of voids at second phase particles[END_REF] criterion determines when strain bands are oriented at 45 • relative to the main loading axis, whereas [START_REF] Mcclintock | A Criterion for Ductile Fracture by the Growth of Holes[END_REF] and [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] (who modified [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] model) consider a critical porosity. Following [START_REF] Needleman | Analyses of plastic flow localization in metals[END_REF] work, a class of criteria determines the instant when strain is no more homogeneous and concentrates in the ligaments between voids. These criteria compare the norm of the strain rate in a localization band and its value outside the band (or the average value throughout the unit cell): if the ratio is higher than an arbitrarily chosen value, failure is said to have been reached. Such criteria are used for example by [START_REF] Barsoum | Rupture mechanisms in combined tension and shear-Micromechanics[END_REF] or [START_REF] Dunand | Effect of Lode parameter on plastic flow localization after proportional loading at low stress triaxialities[END_REF]. Similarly, [START_REF] Luo | On the prediction of ductile fracture by void coalescence and strain localization[END_REF] and [START_REF] Vishwakarma | Micromechanical modeling and simulation of the loading path dependence of ductile failure by void growth to coalescence[END_REF] consider unit cells composed of several layers and force strain localization to happen in the central one (because the external layers contain smaller voids or no voids at all): failure can then be monitored by comparing the behavior of the layers. Another class of criteria determines when a maximum stress or force is reached. Such criteria can be derived by limit analyses, for instance [START_REF] Thomason | Three-dimensional models for the plastic limit-loads at incipient failure of the intervoid matrix in ductile porous solids[END_REF], [START_REF] Benzerga | Effective Yield Criterion Accounting for Microvoid Coalescence[END_REF] or [START_REF] Morin | A unified criterion for the growth and coalescence of microvoids[END_REF]. [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] interpreted the maximum of an effective force in terms of [START_REF] Rice | The localization of deformation[END_REF] criterion on strain localization. Another approach, adopted by [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF] and used for example by [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF] defines coalescence as the transition to a specific strain state: in coalescence, ligaments are in a state of uniaxial straining (whereas the rest of the cell is rigid and hardly deforms). Coalescence could also be interpreted in terms of plastic and elastic energy, as done by [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF]. A last approach was proposed by [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF] and involves computing the macroscopic acoustic tensor in order to directly apply [START_REF] Rice | The localization of deformation[END_REF] criterion on strain localization. However, as pointed for instance by [START_REF] Tekoglu | On localization and void coalescence as a precursor to ductile fracture[END_REF], [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] or [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF], the above criteria actually described two different physical processes: strain localization and coalescence. During strain localization, strain concentrates in narrow bands, which can be interpreted as a loss of ellipticity, according to [START_REF] Rice | The localization of deformation[END_REF] analysis. As stated previously, [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] establish a link between strain localization (through Rice's criterion) and maximum force criteria. Nonetheless, the more direct application of Rice's criterion by [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF] detects localization significantly later than [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] interpretation. On the other hand, coalescence represents the fusion of several voids into a unique larger void during ductile failure. However the material model described in this article contains no ingredient to represent explicitly this process of coalescing voids. The state of coalescence can be deduced nevertheless from the FEM results: at some point in the loading, the cell stops thinning and the plastic flow inside becomes macroscopically uniaxial according to [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF].

As the figure 2.2a shows, the cell's failure, defined by the sudden acceleration of the stress decrease, is incorrectly predicted by the instant of maximum force applied on the cell (with our choice of periodic boundary conditions, this force is here represented by S 11 ). Due to the absence of hardening, the maximum of S 11 happens at the beginning in the simulation, much earlier than the sudden stress drop. On the other hand, this stress drop occurs simultaneously with the stabilization of deformation in the 2-direction transverse to the main loading 1-axis, and can thus be associated with coalescence: the stabilization of the average transverse deformation indicates that the macroscopic strain becomes purely uniaxial. Coalescence thus seems an accurate failure indicator in this situation. According to [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF], an ellipticity loss approach based on the computation of the macroscopic acoustic tensor could also give sensible values of failure strains. This criterion was found to predict slightly earlier failure than a coalescence indicator. However [START_REF] Morin | Numerical study of ductile failure under non-proportional loading[END_REF] tried to apply coalescence and strain localization approaches to match experimental results; both gave acceptable results, with slightly better results for coalescence. Therefore we will focus first on the coalescence approach.

Failure indicator based on the loss of full rankedness of Ḟ

∼

The criterion of the stabilization of transverse displacement, as used by [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF], suffers from two main drawbacks. In a random population of voids, strain localization bands might not be parallel to a face of the cube, so monitoring E 22 with respect to E 11 might not detect coalescence. Moreover, this criterion is limited to the detection of coalescence by internal necking where voids coalesce in the plane orthogonal to the main loading axis. However, for shear dominated loading conditions (when the Lode parameter is close to zero), coalescence is known to occur in shear bands [START_REF] Barsoum | Rupture mechanisms in combined tension and shear-Micromechanics[END_REF][START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF]. We generalize here the stabilization of transverse deformation by noting that for both internal necking and shear, deformation gradient has a specific form during coalescence. After coalescence, there exist orthogonal unit vectors e and e such that F ∼ = 1 ∼ + εt e ⊗ e for uniaxial straining and F ∼ = 1 ∼ + εt (e ⊗ e + e ⊗ e ) for pure shear. In both cases, det Ḟ ∼ = 0. Therefore, as coalescence takes place, det Ḟ ∼ should vanish.

This behavior of det Ḟ

∼ should be compared to the homogeneous plastic deformation case (which is an approximation, since strain may be concentrated in some ligaments). Let us then consider the function:

δ (t) = ε-3 (1 + εt) 3 det Ḟ ∼ (2.15)
which compares the evolution of det Ḟ ∼ to its theoretical evolution in the case of homogeneous compressible plastic flow. A derivation of the expression of δ and an example can be found in appendix 2.A. Therefore, if δ (t) → 0, det Ḟ ∼ decreases faster than expected by homogeneous plastic flow, and localization can be considered to have taken place.

The onset of failure can then be defined as the first instant t c such that:

δ (t c ) ≤ min(A max t<t c δ (t), B) (2.16)
where A = 0.05 is a threshold comparing the maximal and current values of δ and B = 0.005 is an absolute threshold. As shown by the α 2 factor in equation (2.33), δ keeps smaller values for simulations with L close to zero. In these cases, the relative threshold (depending on A) was found to be inappropriate due to numerical errors, and an absolute threshold B (consistent with the value of A) was implemented; it is only needed for loading conditions with |L| < 0.3. A sensitivity analysis with respect to the empirically chosen values A and B is carried out in appendix 2.A and shows that the results which will be presented in section 2.5 are not strongly influenced by the values chosen for A and B. The indicator is therefore robust with respect to the choice of these parameters.

As this criterion using the δ function relies only on macroscopic quantities (at cell-level), it is easy to compute and does not make any assumption on the position and orientation of the possible strain localizations. Moreover it can be used as a landmark in order to stop the simulations shortly after failure in order to spare computation time. However, the indicator detects a loss of full rank of the deformation gradient rate, and is therefore not adapted to loading conditions where the deformation gradient rate is intrinsically of rank 1 or 2. This is especially the case for L = 0 for which the material is initially in shear, so that the indicator is activated in the elastic regime and predicts an early failure. This is acceptable for a perfectly plastic von Mises matrix, but leads to an underestimation of the strain at failure for materials whose hardening behavior delays coalescence. Moreover it is not able to represent a third and rarer form of coalescence known as necklace coalescence. This form was studied by [START_REF] Gologanu | Theoretical models for void coalescence in porous ductile solids. II. Coalescence in columns[END_REF] for a cylindrical unit cell with an axisymmetric loading corresponding to our L = 1 situation. The coalescence between voids takes place along the cylinder axis, which corresponds in our situation to the third and least stressed axis. Necklace coalescence is not associated with a loss of full rank, so the δ indicator cannot be activated. However for the loading conditions involving overall stress triaxiality considered in this work, the proposed indicator has been found to be relevant in all cases.

The failure onset t c can be determined with this method for the different loading conditions, and allows to define several quantities at the onset of coalescence: deformation at coalescence E c = E 11 (t c ), stress at coalescence σ c = σ11 (t c ) and porosity at coalescence f c = f (t c ). In the following, the evolution of those quantities and their dispersion due to the randomness of microstructures will be studied with respect to T and L parameters.

Results

Response of a microstructure subjected to proportional loading with different stress triaxiality and Lode parameter values

The random microstructures constructed in section 2.3.1 have two main differences in comparison to standard unit cells: they contain several voids and these voids are located irregularly within the cell. In this section, the effect of these differences on the behavior of cells is investigated. Several microstructures are considered and subjected to various loading conditions (defined by T and L).

Their failure behavior (evolution of E c , f c and σ c with respect to T and L) are then compared.

The four microstructures shown on figure 2.1 are analysed: two random 27-void cells R1 and R2, a unit cell, and a 27-void cell, lattice, where the voids are distributed following a 3 × 3 × 3 cubic cell. The four microstructures have the same porosity 6%, the same pore size and are meshed with identical mesh size requirements in order to limit the influence of mesh convergence on the comparison (the mesh of the unit cell is thus composed of significantly fewer elements than the other three microstructures). Mesh size will be further discussed in section 2.B.1. The lattice cell allows to separate effects due to the presence of several voids in the cell from those due to the irregular void distribution.

Although the failure behavior on the whole T -L space should be explored, it is instructive to first consider constant triaxiality or constant Lode parameter slices of this space. Let us concentrate first on axisymmetric loading cases characterised by fixed L = -1, as in [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]. We focus on the zone of intermediate triaxiality T ∈ [0.7, 2.0], as usual in unit cell studies [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF][START_REF] Vishwakarma | Micromechanical modeling and simulation of the loading path dependence of ductile failure by void growth to coalescence[END_REF]. We did not study the very low triaxiality regime T < 0.4 where the phenomenon of void collapse takes place [START_REF] Bao | On fracture locus in the equivalent strain and stress triaxiality space[END_REF][START_REF] Liu | Void behaviors from low to high triaxialities: Transition from void collapse to void coalescence[END_REF]. Triaxiality levels T ∈ [0.4, 0.7] were not studied so as to limit the duration of simulations: coalescence generally happens with the same mechanisms as for T > 0.7, but at significantly higher strain values.

The evolution of E c , f c , and σ c with respect to T are shown on the left side of figure 2.3. The four microstructures display globally similar responses: E c monotonically decreases with increasing T while σ c increases linearly with T . f c behaves similarly to E c , except for the microstructure R2: f c is still a mainly decreasing function of T but a local maximum is found at T = 1.2. Note that the evolution of E c with respect to T is smoother and less noisy than that of σ c (for T = 0.8, the stress value for the lattice cell is for instance particularly low, when compared to the values at T = 0.6 or T = 1.0). A possible explanation is that, unlike E 11 which is linearly increasing with time, σ11 and f vary rapidly around the instant of coalescence: σ11 decreases sharply around the coalescence (as evidenced by figure 2.15). Note also that for R2, the porosity at coalescence for T = 1.2 is larger than for T = 1.1, in contradiction to the overall evolution. Coalescence is detected at approximately the same strain in these two conditions, but as the porosity grows faster with increasing triaxiality, the porosity at coalescence is larger for T = 1.2 than for T = 1.1. This slight deviation from the overall evolution with T seems due to the randomness of the void population.

The evolution of strain at coalescence E c was plotted in a logarithmic scale, so as to illustrate the exponential decrease for each microstructure. According to [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF] results, a spherical void typical growth rate varies as exp(3T /2). If we assume that coalescence happens at a given porosity (as for [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]), E c should vary as exp(-3T /2). The evolution of strain at coalescence for the random microstructures, the unit and the lattice cell can be well represented by this simple relation, as shown by the comparison with the straight line of slope -3/2.

The evolution of failure-related quantities with respect to triaxiality, at fixed L = -1 is thus very similar for the various studied microstructures, although some differences are visible. The situation is different if the triaxiality T = 1 is fixed and the coalescence behavior is studied with respect to the Lode parameter (the whole range L ∈ [-1, 1] is explored). The results of the simulations are shown on the right side of figure 2.3. Values for L = 1 are indicated with superimposed arrows but should be taken with caution because, for these loading cases, the simulations diverged or the failure indicator was not reached; the data for the last computed point is indicated to serve as lower or upper bounds for the real value at coalescence, if it exists. The case L = 1, which corresponds to an axymmetric loading where the two largest principal stress components are equal, is associated by [START_REF] Gologanu | Theoretical models for void coalescence in porous ductile solids. II. Coalescence in columns[END_REF] with the necklace coalescence. Our criterion described in section 2.4 is not able to represent this type of coalescence, which does not correspond to a loss of full rank of the deformation gradient rate. Examining the stress strain curve of the unit cell in the case L = 1 (not shown here for brevity) shows a stabilization of stress which could be linked indeed to a coalescence event, undetected by the δ indicator.

If we do not consider anymore the values for L = 1, the unit and lattice cells behave in a similar way (the difference between these two types of cells, which should represent the same void configuration, is due to the meshing): E c increases slowly with L. This type of evolution was reported by [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF] and [START_REF] Zhu | Combined effects of triaxiality, Lode parameter and shear stress on void growth and coalescence[END_REF] and by [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] localization indicator (when σ ∼ does not have shear components). However, [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF], [START_REF] Dunand | Effect of Lode parameter on plastic flow localization after proportional loading at low stress triaxialities[END_REF], [START_REF] Guo | Crystal plasticity based investigation of micro-void evolution under multi-axial loading conditions[END_REF], [START_REF] Zhu | Combined effects of triaxiality, Lode parameter and shear stress on void growth and coalescence[END_REF], and [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] (for the latter two, in more general loading cases), report that E c is a convex function of L, with a minimum near L = 0. Yet, in our case, a sharp decrease in ductility is observed for L close to zero. This behavior in generalized shear corresponds to the expected behavior of a perfectly plastic von Mises material which localizes immediately in shear. However the random microstructures R1 and R2 do not exhibit the same evolution as the unit cell. Three zones can be observed on the E c -L plot for the R1 microstructure (schematized in figure 2.4). The first zone corresponds to L ∈ [-1, -0.7], in which E c increases up to a maximum value on a cusp. For L ∈ [-0.7, 0.4], E c is convex in L and minimal for L = 0. The third zone corresponds to L ∈ [0.4, 1[, where E c decreases from its maximum at L = 0.4 and stabilizes.The zone boundaries correspond to local maxima (significantly higher than the rest of the data points) of E c . They correspond to slope discontinuities, although E c remains continuous. In the following, these three zones will be referred to as: Low Lode parameter Extension Mode Zone (LLEMZ), Shear Mode Zone (SMZ) and High Lode parameter Extension Mode Zone (HLEMZ); the rationale behind these names will be made clearer after section 2.5.2. A similar decomposition in three zones can be seen for the microstructure R2 although for different zone boundaries (L = -0.9 and L = 0.55). To the knowledge of the authors, such an evolution of E c with respect to L was not found in literature. Although [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF]'s coalescence criterion yields a non-smooth evolution of E c , there is only one local maximum, for L > 0 (with our definition of L). An explanation for the existence of these three zones will be proposed in section 2.6.1.

An asymmetry between positive and negative values of L can also be observed: ductility E c is higher in generalized compression than in tension. This asymmetry is present in previously mentioned studies, but the sign of the difference varies among them. Our results are consistent with [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF] and [START_REF] Dunand | Effect of Lode parameter on plastic flow localization after proportional loading at low stress triaxialities[END_REF] but [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF], and [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] have found cells more ductile in generalized tension than in compression (taking into account the different definitions of L). Similar behaviors and differences between the unit and lattice cells on the one hand and the random microstructures on the other hand can be observed on the results for porosity at coalescence. For the stress at coalescence, the asymmetry between L < 0 and L > 0 is clear. There is no significant difference between the microstructures at L > 0, for L < 0; σ c is almost constant for the lattice and unit cells, whereas it increases slightly with L for the random microstructures.

RESULTS

No zone boundaries can be easily identified. The theoretical values for σ c obtained for a von Mises material failing when σ vm = R 0 are also represented. The type of evolution agrees with the results for the unit cells, but due to the porosity and the complex coalescence behavior, stress levels are significantly lower for the cells, and the slope of σ c with respect to T for the simulations at constant L = -1 also differs.

In contrast to the unit and lattice cells, the random microstructures display several zones on their E c -L curve, which could be linked to different coalescence behaviors. The different zones for the microstructure R1 are also shown in the T -L space in figure 2.5. Multiple simulations were carried out for T ∈ [0.7, 1.1]. A simple interpolation using Gaussian Process Regression (as implemented in Scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]) is proposed and allows for an easier visualization in the E c -T -L space, although cusps at zone boundaries are smoothed. The results of the simulations are also projected in the E c -L plane. The triaxiality has two effects on E c : E c globally decreases with higher triaxiality levels, in agreement with the previous study at fixed L = -1, and the position of the zone boundaries is modified (at T = 1.1, the central zone is wider than at T = 0.7).

Relation to localization modes

Several ductility zones were identified on the strain at failure curves for the random microstructure cells, in contrast to unit cells. However the E c curves only give macroscopic information and shed no light on the mechanisms inside the cell responsible for the drastic changes in strain at coalescence. We now investigate the relation between the presence of these zones and the aspect of strain fields inside the cell.

The figure 2.6 shows, for each microstructure, the cumulative plastic strain field p shortly after coalescence (for E 11 1.1E c ). Each image corresponds to a different L value (T = 1 is fixed). The images are inserted on E c -L curves in order to better correlate macroscopic and field information.

All the p fields display zones of higher strain or even strain localization (as localization is known to happen before coalescence [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF]). These zones are organized along approximately planar bands. For both the unit and lattice cells, these bands are exactly planar and correspond to a crystallographic plane of the void lattice. In the lattice cells, the three rows of voids are equivalent, but this symmetry is broken after coalescence. For random microstructures, the bands are more complex: a base plane can be identified but bands are distorted by void distribution so as to include more voids.

For a given microstructure, the orientation of the bands is not constant with L. Two different orientations can be distinguished. In the first one the band is roughly parallel to a face of the cell (and perpendicular to the main loading axis). For the cases with L 1, the localization pattern is more complex and is composed of several bands. The second type of orientations is characterized by strain bands of overall direction approximately 45 • relative to the faces of the cell (although their precise shape is more complex). These two orientations are partly constrained by the periodic boundary conditions because strain localization bands should be compatible with the periodicity of the cell. Notice that bands oriented at 45 • are only found for Lode parameters close to zero (and only for L = 0 in the regular unit and lattice cells) whereas parallel orientation is found for higher values of |L|. Observing more carefully the relation between the orientation of the bands and the macroscopic E c -L curves for the random microstructures shows that strain band orientation is systematically associated with ductility zones: the 45 • orientation is only found in the SMZ whereas parallel orientations are found in the LLEMZ and HLEMZ. Therefore the transition between ductility zones can be linked to a change in strain localization mode: between extension mode, with strain bands at parallel orientation, and shear mode characterized by the 45 • orientation.

To better characterize the transition between ductility zones, as explained by p fields, the similarity between the p field at coalescence for a given value of L and three reference coalescence p fields for L ∈ {-1, 0, 0.9} is here quantified for the R1 microstructure. Each loading case is considered a paragon of its ductility zone (respectively the LLEMZ, the SMZ, and the HLEMZ). If two p fields are similar, they should represent a similar coalescence mechanism. A similarity indicator is defined as follows. The p field after coalescence (at strain E 11 = 1.1E c ), as produced by a FEM computation, is represented by the vector [P] of p values at all Gauss points (for the R1 microscructure, the [P] vectors are around 7 × 10 5 components long). The relative spatial position of Gauss points is irrelevant here. For two vectors [P] and [P ] representing p fields on the same mesh and with the same ordering of Gauss points, the similarity can then be defined as the angle (or rather its cosine) between [P] and [P ]:

cos(θ PP ) = [P] • [P ] ||[P]|| • ||[P ]||
(2.17 .6: Link between E c evolution with respect to L and cumulative plastic strain field shortly after coalescence (at strain E 11 = 1.1E c , which depends on the simulation). Fixed triaxiality T = 1

RESULTS

The evolution of the similarity indicator cos(θ ) to the reference strain fields L = -1, L = 0, L = 0.9 is plotted in figure 2.7. The three reference strain fields are not orthogonal, so significant overlap between the indicators is possible. As strain fields at L = -1 and L = 0.9 are similar (cos(θ ) = 0.85), their similarity indicator shows comparable behavior. However the evolution of the indicator for L = 0 is reversed. The three ductility zones defined earlier are apparent on the figure. For the LLEMZ L ∈ [-1, 0.7], the contributions of L = -1 and L = 0.9 are high and almost constant whereas the contribution of L = 0 is lower but increasing. On the contrary, in the SMZ [-0.7, 0.5], strain fields are predominantly linked with L = 0 and little with L = -1 or L = 0.9. In the last zone, HLEMZ, above L = 0.5, the similarity to the L = 0 strain field decreases, whereas L = -1 and L = 0.9 contributions are higher. Notice however that the L = -1 similarity indicator is high at L = 0.5 and decreases with L, unlike the L = 0.9 indicator. For L 0.5, the situation is close to that of L = -1, whereas at very high L, another mechanism could come into play, especially the competition between two perpendicular strain bands observed earlier at very high L. Around the ductility zone boundaries, strain fields quickly change from one mode to the other. This competition between modes could explain the cusps in strain at failure observed at zone boundaries. L=-1.0 L=0.0 L=0.9

Figure 2.7: For the R1 microstructure, similarity between the coalescence p field at varying L (T = 1 is fixed) and three reference p fields obtained at L = -1, L = 0, L = 0.9.

Dispersion due to microstructure sampling

In the previous two sections, two random microstructures were considered and the evolution of coalescence-related quantities with respect to loading conditions were studied, showing significant differences when compared to the unit cell. Rather than choosing fixed microstructures and varying T and L, another approach is to treat E c , f c and σ c as random variables (depending on the microstructure), and study their statistics. The present section thus focuses on quantifying coalescence onset quantities; the localization mode and the position of the strain localization bands depending on the microstructure are briefly discussed in section 2.C.

N = 20 microstructures with 27 voids and initial porosity f 0 = 6% ( f 0 is not a random variable) were randomly and independently generated. Each of them was subjected to the same loading conditions (T, L) ∈ {(1, -1), (1, -0.5), (1, 0.5), (1.5, -0.5)}. The results for E c , f c and σ c are shown in figure 2.8 as box plots, and are compared to the values for the unit cell. A strong relative dispersion is present for all loading cases: the ratio of the standard deviation to the average is 34%, 59%, 55% and 62% respectively. This indicates a strong influence of the microstructure on the coalescence behavior. The results from unit cells do not represent well the behavior of the multiple void cells, and lead for instance to an overestimation of the stress at coalescence. Dispersion also depends on the loading conditions: for T = 1, the case L = -1 shows lower interquartile range than the cases L = ±0.5. This can be linked to the proximity of zone boundaries for the latter cases, as E c was shown to vary sharply near those boundaries. Moreover, and especially for L = 0.5, some microstructures coalesce in tensile mode whereas others coalesce in shear mode (compare for instance the strain fields of R1 and R2 in figure 2.6); the possibility of different coalescence modes may increase dispersion. At higher triaxiality T = 1.5, L = -0.5, the dispersion is reduced for E c and f c when compared to T = 1, L = -0.5 but the relative dispersion is not. This is due to the overall effect of coalescence appearing earlier at high triaxiality. Besides, the interquartile range for σ c is comparable for both triaxiality levels.

The previous results dealt with a small number of loading conditions. In order to determine an effective model of coalescence in random multiple-void cells for all loading conditions, the T -L space should be explored more extensively, while still keeping a large enough set of microstructure realizations. As in section 2.5.1, multiple simulations were carried out for T ∈ [0.7, 1.1] and L ∈ [-1, 1] on five random microstructures among which R1 and R2 (keeping 20 realizations would have been computationally too expensive). The same loading conditions were tested for each microstructure. The results for E c are shown in figure 2.9. The minimal, maximal and average values are first plotted for T = 1. In agreement with the preceding results, significant relative dispersion is present, and its extent depends on L: dispersion is particularly strong near L = 0.5, whereas it is negligible for L = 0 (all the microstructures agree on almost immediate localization for generalized shear). Despite the dispersion, the overall aspect of the E c curve, as described in the previous section, and its decomposition in ductility zones, are still observable. An interpolation by Gaussian Process Regression of the results in the T -L space is also proposed, based on the average value of the five microstructures at each loading condition. The aspect is similar to that of figure 2.5.

Discussion

In this section, we discuss the significance of the results presented up to now, and assess how representative the results are and how far they can be generalized. First we compare the failure indicator proposed in section 2.4 to [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] coalescence criterion, in order to interpret the difference between unit and random cells. The dispersion due to the random microstructures is then statistically studied with increasingly large void populations. Finally the influence of a work-hardening material is also addressed.

Interpretation of the proposed failure indicator

In order to better understand the failure mechanism identified by the δ indicator and the observed difference between the unit cells and the random microstructures, the δ indicator is compared to another failure criterion reported in the literature. We focus on [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] energy-based coalescence indicator, although Zhu, Ben Bettaieb, et al.'s (2020) and Daehli, Morin, et al.'s (2020) approach with [START_REF] Rice | The localization of deformation[END_REF] criterion could also be useful. According to the former indicator, coalescence is associated with concentration of the plastic deformation in the ligament whereas elastic unloading takes place elsewhere. Therefore coalescence can be detected by monitoring the evolution of the plastic Ẇp and elastic Ẇe work rates and the onset corresponds to the minimum of the ratio Ẇe / Ẇp .

For our cells, the corresponding work rates can be computed by the following equations:

Ẇp = V σ vm ṗ dV (2.18) Ẇtot = V 0 S ∼ : Ḟ ∼ (2.19) Ẇe = Ẇtot -Ẇp (2.20)
The plastic power can be computed either on the cell (with voids) or more easily on the matrix, since stress is zero in the voids. The total power, sum of the plastic and elastic parts, can be computed by only using macroscopic quantities, according to the results of homogenization theory [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]. The figure 2.10 compares the failure onsets, as determined by the δ and the energy-based criteria, with respect to L, for the unit cell and the microstructure R1. For the unit cells, the energy criterion identifies a coalescence onset for all the simulations, and the trend is typical of E c vs. L curves in the literature (for instance [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF]). Moreover the two criteria yield similar values of E c , except for L = 0 where the δ -criterion predicts early failure as previously. The situation is more complex for the random microstructure. For the extension mode zones of the curves, the two criteria also yield very similar results, so they can be thought to represent the same failure mechanism. On the other hand, the energy-based criterion fails to activate in the SMZ, so no coalescence is detected, according to [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] definition. The evolution of the power ratio Ẇe / Ẇp for the unit and the R1 cells, at L = -1 and L 0.2 (in the SMZ for R1), is shown on the figure 2.11a: contrary to the unit cell and the L = -1 case for R1, the L = 0.2 does not show any minimum of the power ratio. In all simulations, the elastic power does not become negative because, unlike [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] unit cell, our microstructures do not possess large void-free regions, in which an elastic unloading can take place. Moreover the simulations also differ by the evolution of porosity (fig. 2.11b): unlike the other three cases, the random microstructure with L = 0.25 does not show any acceleration of void growth during failure, which is typically observed for coalescence. This comparison shows that the δ criterion acts as a coalescence indicator in the LLEMZ and the HLEMZ, and correctly predicts failure in the SMZ according to another mechanism: localization along 45 • bands in shear. As failure should happen quickly in the SMZ for a von Mises matrix, the results obtained for the δ indicator appear more accurate than those for a pure coalescence criterion. The existence of the three ductility zones and the lower ductility in the SMZ could actually be due to boundary conditions. As noted previously, two evolutions of E c with respect to L are reported in the literature: in [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], [START_REF] Dunand | Effect of Lode parameter on plastic flow localization after proportional loading at low stress triaxialities[END_REF], and [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF], strain at coalescence is minimal for L = 0 whereas it increases almost linearly for [START_REF] Zhu | Combined effects of triaxiality, Lode parameter and shear stress on void growth and coalescence[END_REF] and [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF]. The difference between these two groups of studies is that the former consider a shear stress component in equation (2.13). Several loading conditions therefore correspond to the same triaxiality and Lode parameter, and the reported strain at coalescence is the minimum value over all tests at a given (T, L) couple. Coalescence therefore happens earlier than in the absence of shear stress, and this might lead to different responses, as pointed by [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF]. Another point of view is that the cubic unit cells have an anisotropic localization behavior. Although the cubic cell paves space when periodic boundary conditions are enforced, the axes parallel to the sides of the cube remain privileged, and the response of the homogenized material displays anisotropy. As localization bands should be compatible with the periodic boundary conditions, they are always parallel or [START_REF] Coenen | Novel boundary conditions for strain localization analyses in microstructural volume elements[END_REF]. Adding a shear stress component amounts to changing the principal loading directions relatively to the cube, and coalescence occurs when the most favorable band activates.

In the present study, shear stress was not considered but for the random microstructures, it was shown that coalescence can happen either by a localization band perpendicular (corresponding to the LLEMZ and HLEMZ) or oriented at 45 • to the main loading axis (for the SMZ). Therefore the random microstructures appear softer than the unit cells in that they allow several localization band orientations. The resulting response of the cell is then due to a competition between a limited number of coalescence modes (instead of the theoretical infinity of orientations considered by [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF] for instance). The LLEMZ and the HLEMZ correspond then to the evolution shown by Zhu, Engelhardt, et al.'s (2018) study, whereas the response of random microstructures in the SMZ near L = 0 is closer to that of unit cells in [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF]'s study.

Influence of the number of voids

The microstructures considered in the above sections were composed of 27 voids. A small number of voids allows to investigate the effect of a cluster of pores whereas a sufficiently large number can provide results for an effective homogenized material. As pointed by [START_REF] Morin | A unified criterion for the growth and coalescence of microvoids[END_REF], the homogenization theory does not stricto sensu apply to coalescence, which takes place in a small area in the immediate vicinity of voids.

For computational homogenization with a volume element (VE) approach, random microstructures should contain enough voids to reduce the uncertainty due to sampling and limit the influence of boundary conditions (as there is no intrinsic length scale, the size of the VE is only determined by the number of voids it contains). However the computation power required to simulate large cells with many voids, which lead to FEM problems with millions of degrees of freedom, is prohibitive if carried on dozens of loading conditions and microstructures. This problem is in part mitigated by the use of periodic boundary conditions: [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] showed that homogenized properties converge faster with VE size in this case than with kinematic or static uniform boundary conditions. Their study dealt however with elasticity and the extrapolation to coalescence properties is not possible yet. [START_REF] Hure | Yield criterion and finite strain behavior of random porous isotropic materials[END_REF], who carried out simulations of cells with random voids up to coalescence, compared cells with different number of voids (up to 64) and reported that the maximum stress reached during the simulation stabilizes with the number of voids (indicating the existence of a representative volume element), but the stress at coalescence still shows dispersion between realizations. However only five simulations were performed for each number of voids, which is limiting for a statistical analysis of dispersion.

In a complementary approach, we compare the strain at coalescence results for cells with different numbers of voids: 27, 64, 125. All cells are generated with the process described in section 2.3.1 and their porosity is always 6%; the meshing parameters are however adapted so that the ratio between void radius and maximum element size remains constant for all cells. There are typically 2 × 10 5 , 6 × 10 5 and 1 × 10 6 nodes for meshes of cells embedding 27, 64 and 125 voids respectively. As the computational cost of the simulations increases with the number of voids, we only considered two loading conditions T = 1, L = -1 and T = 1, L = -0.5 and a smaller number of 125-void cells than the twenty 27-void cells already used in section 2.5.3. Examples of p fields after coalescence for a microstructure with 125 voids (fig. 2.12) display very complex localization paths between voids, but still show a principal direction parallel to or at 45 • from the faces.

Dispersion results are shown in figure 2.13. For the T = 1, L = -1 case, dispersion is comparable for the three types of cells: a Brown-Forsythe test [START_REF] Brown | Robust Tests for the Equality of Variances[END_REF] was carried out to verify the equality of variances for the 27, 64 and 125-void groups of cells (this test and the following one use the Scipy implementation [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF]). The statistical p-value is 0.19 so the hypothesis of equal variances cannot be rejected. The mean failure strain is significantly lower for 64 and 125-void cells than for 27-void cells, as proven by a one-way ANalysis Of VAriance [START_REF] Heiman | Understanding research methods and statistics: An integrated introduction for psychology[END_REF] between the three groups (p-value of 0.002). However for the T = 1, L = -0.5 loading case, the dispersion is significantly lower for the 64-void cell (Brown-Forsythe test between the three groups: p-value of 0.026). The average failure strain seems to decrease with the number of voids (an ANOVA test could not be performed due to the unequal variances) Therefore failure seems to begin earlier for cells with more voids. This could be explained by the higher probability of a favorable path for a localization bands when the number of voids grows. Variance remains high for all groups of cells, but it is possible that the number of voids reduces dispersion. The simulations evidence that the size of the volume element can exert an influence on the failure results. The above simulations therefore extend Hure's (2021) study with the results from larger and more numerous cells (allowing a statistical analysis) and are in agreement with his findings. More simulations at an even higher number of voids could be carried out to reinforce the statistical significance of the previous conclusions.

Influence of material behavior

The results previously described hold for a perfectly plastic material. However hardening can mitigate the effects of softening due to void growth, and delay coalescence. We here consider two other types of material behavior characterised by their flow stress functions R(p) replacing the constant R 0 used for perfect plasticity in equation (2.4): For the R1 microstructure, at fixed T = 1 and varying L, a comparison of the strain at coalescence E c between the three hardening behaviors is shown in figure 2.14b. On the one hand, for the power law hardening and the slow saturating exponential hardening b = 10, no central SMZ is observed (except a sudden drop near L = 0), and the evolution is quite similar to that observed for unit cells in section 2.5.1. On the other hand, if hardening saturates more rapidly, as for b = 200, the same response as in the perfectly plastic matrix case is obtained. Therefore, hardening seems able to prevent the change of coalescence mode for intermediate values of L, at least if it does not saturate too quickly so as to provide a stabilization effect throughout the deformation process.

R(p) = R 0 + K p n (power law hardening) (2.21) R(p) = R ∞ -(R ∞ -R 0 ) exp(-bp) (saturating exponential) (2.22) with R 0 = 350 MPa, R ∞ =

Conclusion

In the present study, random microstructures made of identical spherical voids within an elastoplastic matrix were generated, and simulated at constant stress triaxiality and Lode parameter with periodic boundary conditions. The FEM simulations were carried out in a large strain framework up to coalescence. The major findings are the following:

1. Failure was identified using an indicator based on the loss of full rank of the average deformation gradient rate, while taking into account the response in case of homogeneous deformation. The results of this indicator are consistent with other indicators reported by the literature but better captures shear dominated localization modes.

2. Random microstructures show two failure modes, that differ by the orientation of the localization band: perpendicular to the main loading axis for an extension mode, or oriented around 45 • for a shear mode. Unlike unit and lattice cells, the shear mode is not limited to the immediate neighborhood of L = 0. The competition between these two modes leads to a non-smooth evolution of the strain at coalescence with respect to the Lode parameter, showing three zones on the E c -L curve, with reduced ductility near L = 0. The difference between unit cells and random microstructures is reduced when the matrix is no more perfectly plastic, due to a stabilizing effect of hardening. However, the response with respect to T is similar for unit cells and random microstructures.

3. When applying the same loading state to microstructures with similar characteristics, a significant dispersion is found in the results (up to 60% of relative dispersion for strain at coalescence). This strong dispersion is also found in simulations with a higher number of voids.

If a model expressing coalescence quantities with respect to loading conditions is desired, using unit cells therefore appears to misrepresent the effective behavior of a material with a complex void distribution, with differences in the general evolution and oversight of the statistical aspects. Care should therefore be taken when applying results on unit cells to more complex applications.

The present work could be extended in several ways. Firstly larger population sizes will be considered based on parallel computing, in order to improve the statistical representativeness of the presented results. Secondly a broader description of the mechanisms of coalescence in random microstructures will be reached by adding a macroscopic shear stress component to the loading state, so as to explore a greater variety of loading paths. Moreover the link between the proposed coalescence indicator and strain localization criteria such as macroscopic or local loss of ellipticity will also be investigated. Finally this work can be the basis to develop and calibrate an effective damage and plasticity model for materials containing randomly distributed pores. [START_REF] Hure | Yield criterion and finite strain behavior of random porous isotropic materials[END_REF] proposed an example of such a homogenized model, but a new model could integrate the effects of the Lode parameter and the dispersion. However simulating enough loading cases and with sufficient statistical representativeness to completely explore the space of parameters is computationally expensive, especially as the effect of initial porosity should be taken into account. Therefore a strategy to construct a surrogate model with as reduced a number of required simulations as possible should be developed.

2.A About the δ indicator

This appendix provides several complements about the δ indicator. Its expression is first derived by computing the homogeneous deformation of Green matrix. An example of application is then presented. Finally a sensitivity analysis regarding the threshold coefficients is carried out.

2.A.1 Derivation of the expression for δ

Consider a perfectly plastic volume element (neglecting here the elasticity) which deforms homogeneously when subjected to the loading conditions (2.13). In order to simply represent the porous nature of the cell, the material behavior will obey [START_REF] Green | A plasticity theory for porous solids[END_REF] isotropic yield criterion (also used by [START_REF] Fritzen | Computational homogenization of porous materials of Green type[END_REF]):

f (σ ∼ ) = σ eq -R 0 σ eq = 3 2 σ ∼ dev : σ ∼ dev +C(tr σ ∼ ) 2 (2.23)
with C a constant (C = 0 corresponds to a von Mises material, and for C = 1/2, there is no lateral contraction of the cube in tension). The other equations in equations (2.4) are unchanged, but they are applied here to macroscopic quantities.

As the material behavior is isotropic, F ∼ stays diagonal in the diagonalizing basis of σ ∼ . Then F ∼ can be written as:

F ∼ =   1 + vt 0 0 0 b 2 (t) 0 0 0 b 3 (t)   (2.24)
where b 2 and b 3 are functions to be determined. As D

∼ = sym( Ḟ ∼ F ∼ -1
), and F ∼ is diagonal, D ∼ can be written as:

D ∼ = Ḟ ∼ F ∼ -1 = diag ε 1 + εt , ḃ2 b 2 , ḃ3 b 3 (2.25)
For a perfectly plastic Green material, the behavior law in (2.4) reads:

D ∼ = ṗ σ eq 3 2 σ ∼ dev +C(tr σ ∼ )1 ∼ = ṗ σ eq 3 2 σ ∼ + (C - 1 2 )(tr σ ∼ )1 ∼ (2.26)
D ∼ is diagonal so there are three constants α 1 , α 2 and α 3 such that:

D ∼ = diag(α 1 , α 2 , α 3 ) ṗ (2.27) α 1 + α 2 + α 3 = 3Ctr σ ∼ σ eq .
(2.28)

Combining (2.25) and (2.27) yields the system:

ε 1 + εt = ṗα 1 (2.29) ḃ2 = ṗα 2 b 2 (2.30) ḃ3 = ṗα 3 b 3 (2.31)
The plastic multiplier is then ṗ = ε/α 1 1+εt and the differential equations can be solved with the initial conditions b

2 (0) = 1, b 3 (0) = 1: b 2 = (1 + εt) α 2 /α 1 b 3 = (1 + εt) α 3 /α 1 (2.32) Finally, det Ḟ ∼ = ε ḃ2 ḃ3 = ε3 α 2 α 3 α 2 1 (1 + εt) α 2 +α 3 α 1 -2 = ε3 α 2 α 3 α 2 1 (1 + εt) -3-3C tr σ σeqα 1
(2.33)

The function comparing the behavior of det Ḟ ∼ and the homogeneous plastic deformation case is then:

δ C (t) = ε-3 (1 + εt) 3+3C tr σ σeqα 1 det Ḟ ∼ (2.34)
The δ criterion used throughout the article is recovered by setting C = 0, which corresponds to the simplified case of a von Mises material. In this case, the criterion depends no more on the applied σ ∼ . The evolution of δ for a simulation with T = 1 and L = -1 (coalescence in uniaxial strain state) is shown in figure 2.15, for two values of C: 0 and 1/2. For both values of C, the vanishing of δ C is simultaneous with the stabilization of transverse displacement. However the sharp drop of δ C allows a more precise numerical determination of the onset of coalescence than the more progressive stabilization of the transverse strain. For C = 1/2, δ 1/2 is approximately constant at the beginning of the simulation, so that the hypothesis of homogeneous flow in a Green volume element (taking into account the porosity) well represents the overall behavior of the cell with a von Mises matrix. However, with C = 0, δ 0 does not depend anymore on the stress state, while still keeping the sudden drop of δ C necessary for the determination of the coalescence onset.

2.A.2 Sensibility analysis regarding the threshold coefficients

Finally we verify that the δ indicator is a reliable indicator of failure by assessing its sensitivity to the choice of the empirically chosen threshold values. As the equation (2.16) shows, the determination of the onset of coalescence relies on two thresholds: a relative one A, which compares the current value of δ to its maximum, and an absolute one B mostly active in shear-like conditions. The values for those were chosen as A = 0.05 and B = 0.005 but a robust indicator should not be too sensitive to these values.

Figure 2.16 compares the effect of different A and B values on the E c -L curve (common triaxiality T = 1, microstructure R1). At constant B, the effect of A is only visible in the HLEMZ and the LLEMZ, and generally negligible. At constant A, B only affects the coalescence strain values in the SMZ. Although a change in B can modify the strain by 0.05, the global aspect of the curve is preserved. The determination of failure by the indicator therefore appears robust with respect to changes in the coefficients.

2.B Effect of different meshing parameters and boundary conditions

In this section, we review the simulation hypotheses and assess their influence on the results presented up to now, showing therefore how representative the results are and how far they can be generalized. First we verify that finite element discretization effects can be neglected, and investigate the effect of different boundary conditions.

2.B.1 Effect of the meshing parameters

All the simulations described up to now were carried out on meshes of cells with the same meshing parameter. To determine the influence of mesh size on coalescence results, the same microstructure R1 was meshed with different meshing parameters h cell /r 0 ∈ {1.25, 1, 0.875, 0.625} (with the notation of section 2.3.1). The maximum element size near the voids is also adapted to keep the ratio h cell /h void = 5 constant. The same loading condition T = 1, L = -1 is applied to the four meshes. Figure 2.17 shows that stress values during the simulations differ between the meshes, but the relative difference between the finest and coarsest meshes is about 5%, which remains acceptable. The onset of coalescence E c which is our main quantity of interest, is almost identical between the meshes, at E c = 0.33 ± 1%. Therefore the influence of mesh refinement for random microstructure cells appears limited (although there was only a ratio of 2 between the element sizes of the coarsest and the finest mesh), which justifies the value h cell = 0.08 adopted throughout this study.

0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35

E 11 1.0 1.1 1.2 1.3 1.4 11 /R 0 h/r = 1.25 h/r = 1.0 h/r = 0.875 h/r = 0.625

2.B.2 Effect of the boundary conditions

We here investigate the influence of boundary conditions. The results from section 2.5 are first compared to those obtained with different boundary conditions. Namely we investigate the influence of conditions on the average gradient, and of planar faces conditions. The consistency of results at L = -1 is also checked by a comparison with simulations on axisymmetric cells. The conditions imposed on the average gradient F ∼ to prevent rigid body motion are first investigated. In section 2.3.3 we imposed F ∼ symmetric, as for [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]. However another reasonable choice would be to fix some degrees of freedom at the vertices of the cubic cell, as depicted in figure 2.18, which is the standard method for boundary value problems. A vertex is already fixed in order to prevent translations, but by fixing two degree of freedom on a second one, and a last one on a third vertex, all rotations are fixed. This can be reformulated as:

F12 = F13 = F23 = 0 (2.35)
i.e. F ∼ is an upper triangular matrix. Due to the mixed conditions imposed by the macroscopic spring element, the results from the symmetric F ∼ case cannot be easily transposed to the triangular F ∼ case. These two choices lead to distinct proportional loading path classes and should therefore be compared.

1 2 3 • • • x x x x x x Figure 2
.18: Conditions on average deformation gradient obtained by fixing some degrees of freedom on vertices of the cubic cell.

On the microstructure R1, at fixed triaxiality T = 1, simulations were performed for several Lode parameters to compare the two sets of conditions on F ∼ (figure 2.19a). The evolution of E c is close between the two types of conditions, and the same ductility zones can be identified for the triangular gradient condition. However, in that case, cusps seem to be less pronounced than for a symmetric gradient; this may be due to the different treatment of shear components by the two conditions. Therefore the influence of the conditions on F ∼ remains limited. We then compare the effects of periodic and parallel faces boundary conditions. Parallel faces conditions mean that the cubic cell retains parallel flat faces throughout the computation (for instance all the points on the x 0 = 0 face have the same x-displacement). This condition is more constraining than periodic boundary conditions. As the comparison in figure 2.19b shows, the two types of conditions lead to qualitatively different responses. For the parallel faces, no separation between three ductility zones can be seen (except near L = 0) and the response of the random microstructure is closer to that typical of the unit cell. Moreover no decrease of ductility near L = 0 is observed for the unit cell. Results for the unit cell differ between the parallel faces and periodic boundary conditions, because in the periodic case, faces are allowed not to remain strictly parallel and planar. On the contrary, boundary conditions made of parallel sides strongly hinder the shear mode failure and only the extension mode remains possible. The competition between these two modes tends to postpone failure (see the cusps on figure 2.4). Therefore, the reduced competition between modes may explain an earlier coalescence for parallel unit cells. The preceding results show that boundary conditions exert a strong influence on the response of the cell.

Finally, the consistency of results obtained at T = 1 is checked. As this type of loading is axisymmetric, a computation with a 2D axisymmetric unit cell was also performed for comparison. Such unit cells are frequent in ductile fracture studies [START_REF] Morin | Coalescence of voids by internal necking: Theoretical estimates and numerical results[END_REF] for instance). The diameter and the height of the cylinder were chosen equal to L cube . The porosity is still 6%, so the radius of the void was modified to 0.22L cube . The boundary conditions for this cell differ slightly from those described in section 2.3.3: they are no more periodic and are replaced by straight edges conditions. Besides the virtual constant triaxiality element is not linked to the average deformation gradient but to the displacement of the top left node.

Figure 2.20 compares results for the unit cell and the 2D axisymmetric cell at varying stress triaxiality for L = -1. In this type of loading, the unit cell was shown in section 2.5 to exhibit the same behavior as random microstructures, compatible with a Rice-Tracey evolution. At initial porosity f 0 = 6%, the axisymmetric cell presents however a significantly higher exponent (in absolute value) for the evolution of E c with respect to T . This effect seems due to the relatively high porosity in the axisymmetric cell: as depicted in figure 2.20, the evolution of E c with respect to T for the low porosity f 0 = 1% axisymmetric cell is much closer to the one predicted by Rice and Tracey.

2.C Position of voids and localization bands

The above study focused on computing the failure strain and therefore quantifying the resistance to ductile fracture of random microstructures for various loading conditions. This section takes a more qualitative approach and provides a brief discussion of the strain concentration bands by which failure takes place, and especially the position of these bands. As noted in subsection 2.4.1, the strain localization bands are mostly orthogonal to the main loading axis, but are not strictly planar for they are attracted by the voids. Their position is therefore guided by the void distribution, but is there a direct relationship between the two? More precisely, can the position of the localization band, and therefore the weakest point of the cell for ductile fracture, be predicted without simulation from the position of the voids. Planes orthogonal to the main loading axis, and with a high concentration of voids, appear especially prone to form localization bands.

This idea is discussed by comparing strain concentration bands in the 20 realizations of random microstructures 27 voids from Fig. 2.9, to the planes with highest porosity given by the distribution of their voids. For each realization, the surface porosity of each plane x = x 0 is computed from the void distribution: This porosity is equal to the surface of the intersection of the plane with the population of defects. Furthermore, the strain localization band for the loading condition (T = 1, L = -1) is considered, and its apparent width is manually extracted. This width is defined as the difference between the abscissas of the leftmost and rightmost points of the failure band (accounting for periodicity). This process is illustrated in figure 2.21. The position of the failure band agrees with the idea of planes of maximum porosity, if the plane with maximum surface porosity belongs to the localization band. This agreement can be quantified by a score:

min x b ∈band 1/2 -1/2 1{ f s (x b ) ≥ f s (x)} dx (2.36)
where f s corresponds to the surface porosity, and the cell is located between the planes x = -1/2 and x = 1/2. This score corresponds to computing the cumulative distribution function of f s and taking the minimum over the failure band. If the band contains the plane of maximum porosity, the score will be zero. For several realizations, the comparison between the strain localization band and the surface porosity level is shown in Fig. 2.22. Fig. 2.23 also shows box plots for the apparent failure band width (relative to the size of the cube) and the score of the 20 realizations under consideration. For most of the realizations, failure takes place near the plane of maximum porosity (for instance in figs. 2.22a and 2.22b). However the failure band is often of significant apparent width, which shows that is not really planar but zigzags and connect several voids. Moreover in several microstructures, the plane of maximum porosity is not contained in the localization band, because a cluster of voids elsewhere provides a more favorable failure path (as in figs. 2.21b and 2.22c). Although an increased surface porosity appears to promote coalescence, failure does not systematically take place where surface porosity is maximal. Predicting the failure path of a given microstructure without simulating the complete response appears therefore challenging.

Chapter 3

Strain localization analysis of materials with randomly distributed voids

Strain localization is often a precursor to the ductile failure of the material. The chapter investigates this phenomenon in the case of random microstructures, namely cubic cells made of an elastic-perfectly plastic matrix embedding distribution of identical spherical voids. Random microstructures allow a better representation of the interaction between voids and greater diversity of failure modes than single-void (or unit) cells. The cells are simulated by means of the finite element (FE) method for proportional stress loading paths with controlled stress triaxiality and Lode parameters. Strain localization is detected with Rice's criterion computed at the cell level. This criterion is shown to accurately capture the onset of localization and the type of failure mode, either extension or shear banding. Moreover, the influence of the loading orientation, i.e. the orientation of the principal frame of the applied stress with respect to the microstructure, is systematically studied. Significant anisotropy of failure behavior is observed, especially in the case of single void unit cells, which can be attributed to the intrinsic anisotropy of the simulation cells. Finally minimal failure strain values at localization on all loading orientations are found. A zone of reduced ductility is observed under generalized shear loading conditions. This chapter is mainly based on the homonymous article, published in the Journal of the Mechanics and Physics of Solids [START_REF] Cadet | Strain localization analysis in materials containing randomly distributed voids: Competition between extension and shear failure modes[END_REF]. A new appendix on the localization in porous balls was added.

57

Résumé détaillé

Le chapitre précédent a montré l'utilité de considérer des microstructures aléatoires, car des comportements différents entre cellules unitaires et aléatoires sont possibles. Cependant, les chargements étudiés jusqu'à présent étaient diagonaux, ne permettant pas de choisir l'orientation du chargement par rapport à la cellule cubique. Comment évolue cette différence de comportement entre les deux types de cellule lorsque des chargments d'orientation quelconque sont considérés ? Par ailleurs, un critère de coalescence, principalement phénoménologique a été utilisé jusqu'ici. Un critère de localisation, davantage justifié théoriquement et plus adapté à des chargements complexes, est donc également considéré. L'objectif de ce chapitre est donc d'étudier la localisation en rupture ductile dans des microstructures aléatoires, soumises à des chargements d'orientation quelconques.

La méthodologie reprend globalement celle du chapitre précédent, mais les conditions de chargement sont généralisées de façon à permettre des orientations quelconques : le tenseur des contraintes moyennes n'est plus nécessairement diagonal, c'est à dire que les directions principales du tenseur ne coincident plus nécessairement avec les axes de la cellule cubique. Par ailleurs, la ruine de la cellule est désormais détectée par le critère de localisation de Rice. Ce critère nécessite de déterminer l'opérateur tangent macroscopique reliant les déformations et contraintes moyennes à l'échelle de la cellule.

Le critère de Rice est d'abord testé pour des cellules soumises à un chargement diagonal. Comme dans le chapitre précédent, la déformation à rupture en fonction de la triaxialité T et du paramètre de Lode L est étudiée : elle décroît avec T et est minimale pour L ∼ 0 lorsque le mode de rupture est en cisaillement (à la fois pour les cellules unitaires et aléatoires). Contrairement au critère en coalescence utilisé au chapitre précédent, le critère de Rice permet de déterminer précisément la nature du mode de rupture : extension, cisaillement ou une combinaison des deux. Une zone de ductilité réduite est associée au mode de cisaillement.

L'effet de l'orientation du chargement est ensuite étudiée. Pour des chargements de mêmes T et L, des orientations différentes peuvent induire une grande dispersion de la déformation à rupture. Ce phénomène se produit dans les cellules unitaires et de façon un peu plus réduite dans les cellules aléatoires. Cette dispersion des résultats avec l'orientation du chargement est le signe d'une anisotropie en rupture ductile, due à la cellule cubique et aux conditions aux limites périodiques. Pour obtenir un modèle de rupture ductile isotrope, il convient de s'intéresser à l'orientation de chargement minimisant la déformation à rupture. Cette déformation minimale est toujours associée à une rupture en cisaillement.

Une dernière partie approfondit la discussion du critère de localisation en le comparant plus précisément aux critères utilisés dans le chapitre précédent. De plus l'anisotropie des cellules est précisée, en montrant qu'elle ne concerne pas seulement la déformation à localisation, et qu'elle reste présente, à un niveau plus réduit, pour des cellules de porosité initiale plus faible.

Introduction

In the preceding chapter, random microstructures were simulated up to coalescence and the scatter of failure strain due to the random distribution of voids was investigated. A significant scatter was found even for large populations. Moreover, by considering the influence of the Lode parameter, differences between unit and random cells were identified in failure behavior. Random cells can more freely express shear-dominated failure modes. Therefore random microstructures can better represent failure mechanisms and they should be favored for determining failure models.

As already pointed, a significant issue of cell studies is to operationally detect the onset of failure. Many indicators have been proposed, focusing on different physical mechanisms and characteristic features. These different choices are not equivalent and a given microstructure's resistance (strain at failure) can significantly depend on the considered indicator. [START_REF] Tekoglu | On localization and void coalescence as a precursor to ductile fracture[END_REF] showed that cell failure could actually be linked to two distinct processes: strain localization and void coalescence. Void coalescence can be identified by a transition to a specific strain state, where ligaments are uniaxially strained and the rest of the cell becomes rigid. Operationnally, the cell ceases thinning in the directions transverse to the main loading axis, and reaches a state of simple extension. This approach was used by [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF] and [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]. The previous chapter reformulated this indicator in terms of the macroscopic deformation gradient rate (and especially its determinant) and proposed an extension of this indicator to handle shear modes of failure.

On the other hand, strain localization corresponds to the concentration of deformation within bands, generally containing voids, as an increased local porosity facilitates localization. This process can be linked to the loss of ellipticity of the governing partial differential equations. This mathematical characterization of localization can be traced to [START_REF] Hadamard | Application à la théorie de l'élasticité[END_REF]'s works on the stability of elastic materials. His interpretation of a discontinuity surface as a zero velocity wave was then generalized to increasingly complex constitutive relations [START_REF] Thomas | Plastic Flow and Fracture in Solids[END_REF][START_REF] Hill | Acceleration waves in solids[END_REF][START_REF] Mandel | Conditions de Stabilité et Postulat de Drucker[END_REF] including non normal relations or yield surfaces with vertices [START_REF] Rudnicki | Conditions for the localization of deformation in pressuresensitive dilatant materials[END_REF][START_REF] Stören | Localized necking in thin sheets[END_REF]. The generalization of the strain localization criterion to arbitrary constitutive laws was then performed by [START_REF] Rice | The localization of deformation[END_REF]. This criterion, henceforward referred to as Rice's, is mostly used for structural computations (see for instance Al [START_REF] Kotob | A general and efficient multistart algorithm for the detection of loss of ellipticity in elastoplastic structures[END_REF]) but has also been recently used for unit cell studies by [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF]. [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] also claimed their maximum applied force criterion is equivalent to Rice's analysis. [START_REF] Tekoglu | On localization and void coalescence as a precursor to ductile fracture[END_REF] and [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] showed that coalescence, if it happens, is always preceded by strain localization. [START_REF] Morin | Numerical study of ductile failure under non-proportional loading[END_REF] tried to match experimental ductile failure results with localization and coalescence and found a slightly better agreement with coalescence, so that structural failure may be better linked to the coalescence process. However localization criteria may provide a relevant lower bound for ductile failure resistance. The coalescence of random microstructures was studied in the previous chapter, but the localization behavior of such microstructures remains to be investigated.

If a full response surface for ductile fracture is desired, sufficiently general loading conditions should be used. For an isotropic material model subjected to a proportional loading, stress triaxiality and Lode parameter are mainly used. Yet conflicting results with respect to the influence of the Lode parameter on the failure strain are reported in literature. [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] and [START_REF] Luo | On the prediction of ductile fracture by void coalescence and strain localization[END_REF] found minimal ductility for generalized shear. On the other hand failure strain grows from generalized tension to generalized compression, without a minimum in shear, according to [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] and [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF]. In chapter 2, the behavior was mostly of the former type for random microstructures, and mostly of the latter type for unit cells. The discrepancy, which can be interpreted as a strong indication of the intrinsic anisotropy of the simulation cell, is due to different ways of applying the loading conditions, as pointed out by [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF] and in chapter 2. However, a reliable failure model should not present this anisotropy and this dependence on the methodology. In order to avoid this effect, [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], [START_REF] Dunand | Effect of Lode parameter on plastic flow localization after proportional loading at low stress triaxialities[END_REF] and [START_REF] Tekoglu | On localization and void coalescence as a precursor to ductile fracture[END_REF] carried out the simulations by rotating the principal directions of the applied stress. The failure strain of the cell is then found by minimizing over all loading conditions with the same stress invariants. A minimal failure strain for generalized shear is then found. However only one axis of rotation was considered in these studies. This choice is motivated by [START_REF] Rudnicki | Conditions for the localization of deformation in pressuresensitive dilatant materials[END_REF] results on localization, which show that localization bands in homogeneous materials should be orthogonal to the eigendirection associated with the middle principal stress. However it was never verified numerically for unit cells. Moreover the situation for random microstructures could be more complex than for unit cells, which have higher symmetry. Fully general three-dimensional loading orientations, i.e rotations of the principal axes of the applied stress with respect to the cell's axes, should therefore be considered.

The present study therefore aims to investigate ductile fracture and more precisely strain localization behavior in random microstructures, while fully accounting for the effect of the loading orientation. To the authors' best knowledge, this is the first time that localization is studied in cells with a distribution of voids. Moreover, no previous studies had considered the effect of loading orientation with such generality. To this end, we generate cells consisting of a perfectly plastic matrix embedding a random distribution of identical spherical pores. Using the finite element software Zset (2022), simulations are performed in a large strain formulation. Various proportional loading conditions are applied up to the failure of the cell. Following the methodology by [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF], this failure is detected by Rice's localization criterion, computed at the global scale of the cell with a macroscopic tangent operator. The failure behavior for random microstructures and unit cells is then compared depending on the loading conditions. It is shown that Rice's criterion successfully describes not only localization onset, as performed by [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF], but also failure mode. A region of low ductility is observed in generalized shear. A strong anisotropy of failure behavior is observed for unit cells but is reduced for random microstructures. Minimal failure strains depending on the loading conditions can then be determined. This chapter is organized as follows. We first describe the methodology for generating random microstructures, applying loading conditions, and detecting localization with Rice's criterion. Secondly, the performance of the criterion is assessed on simplified loading conditions, in which the principal axes of the applied stress coincide with those of the cubic cell. Thirdly, more general loading orientations are investigated. Finally, with the help of further computations, we discuss the preceding results, especially the failure indicator used, and the impact of cell anisotropy. Appendices present the numerical validation of the failure indicator, the effect of the simulation parameters (mesh size, and temporal discretization), and check the limited role of local loss of ellipticity in the matrix.

METHODOLOGY OF MICROMECHANICAL FINITE ELEMENT SIMULATIONS 61

Methodology of micromechanical finite element simulations

This section presents the general methodology for micromechanical simulations, carried out with Zset software [START_REF] Besson | Object-Oriented Programming Applied to the Finite Element Method Part I. General Concepts[END_REF]Zset, 2022). The generation and meshing of the random microstructures, the material behavior used for the matrix, the boundary and loading conditions are described. This methodology follows initially that of chapter 2, it is here summarized for reference and practicality. The subsection 3.3.3 presents however an extension of the previous methodology allowing a greater variety of loading conditions when compared to the previous chapter.

Generation of microstructures and meshing

The microstructures are made of a periodic population of N voids identical non-overlapping spherical voids of radius r within a cubic matrix of size a cube . The total initial porosity

f 0 = 4π 3 N voids (r/a cube ) 3 (3.1)
and the number of voids N voids determine together the radius of the spheres. A sphere intersecting the cube is copied on the other side (two copies for an intersected face, four for an edge and eight for a vertex). The periodic microstructure with its spherical voids is thus able to pave space. The positions of the spheres are chosen randomly according to a simple iterative dart-throwing process: A position of the center is chosen uniformly on the cube, and the new sphere is added to the collection if it (or its periodic copies) does not intersect any already accepted sphere.

The microstructure is then meshed by the NETGEN software Schöberl, 1997, ensuring finer elements close to the voids. Most of the simulations on random microstructures will be performed on the cells R1 and R2 visible in Fig. 3.1a and 3.1b. They were constructed with a porosity of f 0 = 6% and N voids = 27 voids and have already been used in chapter 2. This rather high void volume fraction, which can be found in some materials such as nodular cast iron [START_REF] Zhang | Ductile fracture of materials with high void volume fraction[END_REF], induces significant interaction between voids even during the growth phase [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF]. Considering such a high porosity level is also useful when studying the coalescence of cells with lower initial porosity after sufficient void growth. A unit cell with a single void, is also used for comparison (Fig. 3.1c). All meshes use quadratic elements with reduced integration.

Material behavior law at finite strain

Simulations use a finite strain framework to account for the large deformation that can take place in the matrix. The strain rate D ∼ and Cauchy stress σ ∼ tensors are convected in a corotational frame [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]:

ė ∼ = R ∼ T D ∼ R ∼ s ∼ = JR ∼ T σ ∼ R ∼ (3.2) where R ∼ is a rotation matrix verifying -R ∼ T Ṙ ∼ = Ṙ ∼ T R ∼ = skew( Ḟ ∼ F ∼ -1
) (material spin tensor). The matrix is chosen elastic-plastic with isotropic elasticity and von Mises plasticity without This choice of an elastic-perfectly plastic constitutive law plays a key role in the ductile failure of the cell. Such a material is prone to early coalescence for loading in generalized shear. The presence of a slow or non-saturating hardening can reduce this effect.

ė ∼ = ė ∼ e + ė ∼ p e ∼ e = 1 + ν E s ∼ - ν E (tr s ∼ )1 ∼ s vm = 3 2 s ∼ dev : s ∼ dev f (s ∼ ) = s vm -R 0 ≤ 0 ė ∼ p = ṗ ∂ f ∂ s ∼ (3.

Boundary and loading conditions

Periodic boundary conditions are applied on the sides of the cube [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]. The displacement field u has the form:

u = ( F ∼ -1 ∼ ) • x 0 + v (x 0 ) = E ∼ • x 0 + v (x 0 ) (3.4)
with F ∼ the average deformation gradient, E ∼ = F ∼ -1 ∼ and v a periodic displacement fluctuation field with zero average gradient over the cell. The periodicity of v and the anti-periodicity of traction vectors mean that:

v (x + 0 ) = v (x - 0 ) (3.5) σ ∼ • n (x + 0 ) = -σ ∼ • n (x - 0 ) (3.6)
if x + 0 and x - 0 are two homologous points on opposite sides of the periodic mesh and n (x 0 ) is the outward-pointing normal to the mesh boundary at x 0 . These periodic boundary conditions are numerically imposed by multi-point constraints on homologous nodes. The components of the macroscopic deformation gradient therefore appear as explicit degrees of freedom. In order to control rigid body motion, a node of the mesh is fixed, and the tensor F ∼ (or equivalently E ∼ ) is kept upper triangular:

E 21 = E 31 = E 32 = 0 (3.7)
This condition was investigated in the previous chapter and found to yield equivalent results to the condition imposing F ∼ to remain symmetric [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]. In order to facilitate the ulterior computation of Rice's criterion, this condition is applied by penalization. For (i, j) ∈ {(2, 1), (3, 1), (3, 2)}, a high stiffness spring connects the DOF E i j to a fixed one Êi j = 0. The force applied on the DOF E i j is then:

R E i j = k( Êi j -E i j ) (3.8)
where k is a stiffness (chosen sufficiently high compared to the stiffness of the simulated cube).

The macroscopic Boussinesq (or first Piola-Kirchhoff) S ∼ and Cauchy stress tensors σ ∼ are defined by:

S ∼ = 1 V 0 V 0 S ∼ dV 0 = 1 -f 0 V 0 V matrix 0 S ∼ dV 0 (3.9) σ ∼ = 1 J S ∼ . F ∼ T (3.10)
where J = det F ∼ and V 0 is the volume of the cell (matrix and defects) in the initial configuration. The integration on V 0 considers that stress is well-defined and identically zero in the voids. The components of S ∼ also correspond, up to a factor V 0 , to the reaction forces conjugate to the degrees of freedom E i j .

The loading conditions are applied by an extension of Ling et al.'s (2016) method. A special spring element with ten degrees of freedom is used: Nine are connected to the E i j degrees and the last one, Ê11 , corresponds to the displacement of a ghost node. This element guarantees the following form for the macroscopic Boussinesq stress tensor:

S11 = k( Ê11 -E 11 )/V 0 (3.11) S ∼ = J σ11 η ∼ F ∼ -T (3.12)
where k is a stiffness which can be chosen equal to that of Eq. (3.8), and η ∼ is a constant symmetric tensor with η 11 = 1. Eq. (3.12) can be rewritten in terms of E i j , η i j and Ê11 and K by using Eq. (3.11) and explicitly inverting F

∼ . The derivation is straightforward but tedious so the resulting expressions were obtained with the SymPy computer algebra software [START_REF] Meurer | SymPy: symbolic computing in Python[END_REF] and are not reproduced here. With this spring element, the average Cauchy stress tensor remains proportional to the constant tensor η ∼ throughout the simulation:

σ ∼ = σ11 η ∼ (3.13)
The loading is then driven by Ê11 , which essentially coincides with E 11 due to the high stiffness k.

Zhu, Ben Bettaieb, et al., 2020 and chapter 2 studied the case where η ∼ is diagonal (i.e. η 12 = η 13 = η 23 = 0). [START_REF] Barsoum | Rupture mechanisms in combined tension and shear-Micromechanics[END_REF], [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] and [START_REF] Liu | Void behaviors from low to high triaxialities: Transition from void collapse to void coalescence[END_REF] used conditions equivalent to η 22 = η 33 and η 12 = η 23 = 0. Following [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF] example, several studies have considered rotating the applied stress σ ∼ around e 2 . The problem is here considered with more generality, with a potentially full η ∼ tensor. As η ∼ is a symmetric tensor, it can be diagonalized in an orthonormal frame and its components can be more usefully expressed as:

[η ∼ ] =   1 η 12 η 13 η 12 η 22 η 23 η 31 η 32 η 33   = Q   σI / σ11 0 0 0 σII / σ11 0 0 0 σIII / σ11   Q T (3.14)
where σI ≥ σII ≥ σIII are the principal stresses and Q is a rotation matrix defining the loading orientation (that is, the orientation of the eigendirections of η

∼

) with respect to the frame canonically linked to the cube.

From the principal stresses, which are linked by the condition η 11 = 1, two stress invariants can be extracted. The stress triaxiality T and the Lode parameter L are here defined as:

T = tr σ ∼ 3 σvm (3.15) L = 2 σII -σI -σIII σI -σIII (3.16)
where σeq is the von Mises equivalent stress. In our convention, the values L = -1, L = 0 and L = 1 respectively correspond to states of generalized tension, shear and compression. Therefore a loading condition can be defined by T , L and a 3D orientation matrix Q. Conversely, given a triplet (T, L, Q), the correct tensor η ∼ can easily be found by computing the principal stresses [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF], applying the rotation Q, and normalizing so that η 11 = 1.

Formulation and implementation of the macroscopic Rice criterion

The simulations are primarily aimed at determining the random microstructures' resistance to ductile failure. As the material model defined in subsection 3.3.2 does not account for damage, failure should be investigated at the cell level. An indicator is required to identify failure. In this work, the localization approach is used by applying Zhu, Ben Bettaieb, et al.'s (2020) methodology to the random microstructures. Rice's indicator could also be better suited to complex loading conditions, as illustrated in section 3.5.

Let L be the macroscopic tangent operator linking the macroscopic Boussinesq stress and deformation gradient rate tensors:

Ṡ ∼ = L : Ḟ ∼ (3.17)
Rice's criterion predicts strain localization when, for some direction n the acoustic tensor along n becomes singular: det n k Lik jl n l = 0 (3.18)

In the following, the expression det n L n will refer to the above equation with the correct choice of indices. Whereas the local tangent operator L is directly given by the material behavior, the macroscopic operator L should be computed at the level of the cell in order to link Ṡ ∼ and Ḟ ∼ . It therefore combines the local material behavior and structural effects. Zhu, Bettaieb, et al., 2020 compared several numerical techniques to compute this macroscopic tangent operator. A condensation method derived from [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF] but tailored to our finite element formulation is here used. The general idea of the condensation technique is to compute a Schur complement on the finite element matrix separating macroscopic and local degrees of freedom.

More precisely, from subsection 3.3.3, three types of DOF can be distinguished in the simulation : (i) E i j , macroscopic DOF corresponding to the component i j of F ∼ -1 ∼ ; (ii) v n,i , nodal DOF corresponding to the displacement fluctuation in the direction i of the node n; (iii) four ghost nodes Ê11 , Ê21 Ê31 , Ê32 useful to control the applied stress or the rigid body motion. The vector collecting the E i j (resp Êi j , v n,i ) DOF is noted E (resp. Ê, V ).

The stiffness matrix K computed during the Newton iterations of the FEM simulation, can thus be decomposed in blocks so as to verify:

  K Ê Ê K ÊE 0 K E Ê KEE K EV 0 K V E K VV     ∆ Ê ∆E ∆V   =   ∆R Ê ∆R E ∆R V   (3.19)
The right hand side represents the external forces on each DOF. R V corresponds to the nodal forces, whereas [R E ] i j = V 0 Si j with S ∼ the macroscopic Boussinesq stress tensor applied to the cell, and V 0 the initial volume of the cell. The blocks K ÊV and K V Ê are zero because the ghost nodes are not linked to the displacement fluctuation DOF. Moreover, KEE is obtained by combining a contribution K EE from the tetrahedral elements of the mesh and contributions from the springs:

KEE,i j,kl = K EE,i j,kl + ∑ s spring K (s) E i j ,E kl (3.20)
where K (s) E i j ,E kl is the elementary stiffness matrix of the spring. Up to a factor equal to the volume of the cell, K EE represents the average local tangent operator on all tetrahedral elements. For each spring controlling the rigid body motion, only one coefficient in the block EE is non-zero:

K (s) E 21 ,E 21 = k for the E 21 spring, for example. Similarly K (s) E 31 ,E 31 = K (s) E 32 ,E 32 = k.
The expressions for the extension of [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF] special spring element are more complicated but always computable, provided Ê11 and F ∼ are known. The tilde on KEE represents the modification of the quantity K EE linked to the physical mesh by contributions from the springs.

The effect of Dirichlet boundary conditions and multi-point constraints on K should now be considered. Dirichlet boundary conditions, which are applied on the ghost nodes, lead to the elimination of DOF. For multi-point constraints applying the periodic boundary conditions of Eq. (3.6), rows related to homologous DOF should be summed and combined in a single new row (and the same for columns). Eliminating the Ê DOF due to Dirichlet boundary conditions and combining rows and columns for homologous V DOF is equivalent to modifying the matrix:

K = KEE K EV K EV K VV ∆E ∆V = K -1 ∆R E ∆R V (3.21)
The block KEE has not been modified and the DOF Ê have disappeared. Had the conditions E 21 = E 31 = E 32 = 0 not been imposed with springs, these DOF would have disappeared from the stiffness matrix. Computing the related coefficients in the macroscopic tangent operator would have been impossible.

Computing the macroscopic tangent operator means determining the linear relationship between ∆ F ∼ and ∆ S ∼ , and therefore between ∆E and ∆R E , when no other force is applied. In this context, the nodes of the meshed cell are free to move so that ∆R V = 0 at equilibrium. In Eq. (3.19), this implies K V E ∆E + K VV ∆V = 0, hence:

∆R E = K∆E with K = KEE -K EV K -1 VV K V E = ([K -1 ] EE ) -1 (3.22)
In the last expression, the EE block of K -1 is extracted. The macroscopic tangent operator should then be computed from K by removing the contribution from the springs. Otherwise, even for a homogeneous cube, it would not coincide with the local tangent operator. L is then obtained as:

Li jkl = KE i j E kl /V 0 with K = K EE -K EV K -1 VV K V E (3.23)
K/V 0 is therefore different from K EE /V 0 , which is the average value of the tangent operator over all the elements of the cell. It includes the structural effect -K EV K -1 VV K V E . In practice, at the end of each increment of the Z-set simulation, a Python script with references to Z-set's internal variables is called. This script can access the global stiffness matrix with boundary conditions K . The quantity [K -1 ] EE can then be computed by solving the following equation for the nine E i j DOF:

K α∆E α∆V = α∆R i j E 0 (3.24)
where all the coefficients of ∆R i j E are zero, except the coefficient related to E i j : R i j E,i j = 1. The quantity α is a very large value, allowing to neglect the modifications of the right hand side made by the Dirichlet boundary conditions. The matrix [K -1 ] EE is then constructed from the nine solutions by taking for each column, only the coefficients from the E block. Solving these nine equations with Z-set's solver is computationally efficient. No inversion of the matrix K is necessary, and only one additional factorization of K needs to be performed at the end of each time increment.

The matrix [K -1 ] EE is then inverted by means of the Scipy solver. This inversion is inexpensive as [K -1 ] EE is only of size 9 × 9. The resulting matrix K can then be stored in a separate file. This process is repeated for each increment of the simulation. Finally, as a post-processing operation, K is computed from K by subtracting the contribution of the springs.

The macroscopic tangent operator L has now been computed, and Eq. (3.18) remains to solve. Rather than solving it directly, the minimum of det n L n is found for n on the unit sphere. By continuity, if this minimum is nonpositive, Eq. (3.18) has a solution, and the cell exhibits strain localization. The BFGS minimizer from SciPy [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF] is used with the multi-start procedure described by Al [START_REF] Kotob | A general and efficient multistart algorithm for the detection of loss of ellipticity in elastoplastic structures[END_REF].

For a typical simulation of a random cell, the Rice analysis can lead to an increase of CPU time of approximately 60% when compared to a simulation without Rice analysis. This is principally due to the additional factorization of K , whereas the overheads due to the forward and backward substitutions in Eq. (3.24) and the minimization of det n L n are negligible for sufficiently large meshes.

Appendix 3.A validates the methodology for computing L by comparison with theoretical results and other methods. Note than L only computes the consistent macroscopic tangent operator used in numerical computations, and not directly the instantaneous tangent operator corresponding to the material formulation [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]. However the validation from appendix 3.A proves that the difference can be neglected for sufficiently small steps.

Investigation of localization with Rice's criterion for diagonal loading

Unit cell studies frequently consider the application of diagonal macroscopic stress tensors. To compare with literature, the analysis is first restricted to diagonal η ∼ loading. The loading is then characterized only by T and L. Rice's criterion's is shown in the following to efficiently detect localization onset and failure mode for random microstructures.

Strain at localization depending on T, L loading conditions

The methodology for Rice's criterion application is first verified on the simple case of the unit cell loaded with T = 1, L = -1. This loading condition corresponds to generalized uniaxial tension. For such a loading condition, the failure onset directly corresponds to an inflection point in the stress-strain curve. It is also well understood by applying [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF] criterion: the cell stops thinning and the macroscopic transverse strain E 22 (or equivalently E 33 ) stabilizes. Fig. 3.2 compares the evolution of three quantities throughout the computation, parametrized by the macroscopic strain E 11 : the transverse strain E 22 , the macroscopic stress component σ11 along the loading axis, and the Rice quantity min det n L n throughout the computation. For E 11 0.28, min det n L n = 0, which corresponds to Rice's criterion activating. This localization is simultaneous with the stabilization of the transverse strain and the inflection point of the strain-stress curve. The three failure indicators therefore yield the same failure onset, which confirms the soundness of the methodology using Rice's criterion. This allows defining the strain at localization E loc as:

E loc = arg min E 11 {E 11 | min det n L n ≤ 0} (3.25)
The dependence of the strain at localization on T and L is depicted in Fig. 3.3 for the unit cell and the random cells R1 and R2. Each subfigure considers a slice of the (T, L) space: A parameter is fixed whereas the other one varies. At fixed L = -1 (axisymmetric, generalized tension) or L = 0.5, failure strain decreases almost exponentially with stress triaxiality, for all cells. Such a decrease mostly agrees with the evolution E loc ∝ exp(-3T /2) from [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF], although the agreement, especially regarding R2, is slightly reduced compared to the δ coalescence indicator of chapter 2.

For fixed T = 1 and varying L (Fig. 3.3c), the situation is more complex and differs between cells. For the unit cell, three zones can be distinguished. For L < -0.5 or L > 0.3, failure strain evolves slowly with L and the plot forms a plateau. A slow increase of E loc with L is frequently found in unit cell studies with diagonal loading [START_REF] Zhu | Combined effects of triaxiality, Lode parameter and shear stress on void growth and coalescence[END_REF][START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF]. On the contrary, between L = -0.5 and L = +0.3, E loc is a convex function of L, minimal for L = 0 where localization is almost immediate. Early localization in shear is expected for a homogeneous perfectly plastic material. The role of the matrix's behavior in the localization is studied in more detail in 3.C. For the random cells, the central zone of reduced ductility with a minimum for L = 0 is wider than for unit cells. The plateaus at high |L| are reduced and even disappear for the cell R2. For the cell R1, the distinction between the plateaus and the reduced ductility zone are still present. A local maximum associated with a strong slope discontinuity is present at L ∼ 0.5 and marks the transition between the central low ductility zone and the upper plateau. For random cells, the reduced ductility zone is therefore dominant. Such an evolution is reminiscent of the results by [START_REF] Barsoum | Rupture mechanisms in combined tension and shear-Micromechanics[END_REF], [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF] and [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] or [START_REF] Luo | On the prediction of ductile fracture by void coalescence and strain localization[END_REF], who considered non diagonal loading conditions.

The difference with [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF] who also used a localization criterion and a diagonal loading can be explained by the present use of a cubic cell (instead of a parallelepipedic cell by Zhu, Ben Bettaieb, et al.) which is more sensitive to shear localization.

The distinction between plateaus and a central reduced ductility zone was already predicted by the δ coalescence criterion. However, the localization criterion used in the present study, predicts a wider reduced ductility zone, especially for the unit cell. In the previous article, that zone was limited to the immediate vicinity of L = 0. Moreover, the plateaus at high |L| are not necessarily present with the localization criterion. A more precise comparison between these failure indicators will be performed in subsection 3.6.1.

Localization bands and failure mechanisms

In the previous section, Rice's criterion was used as an indicator for failure onset. A complex dependence of E loc with respect to L was found with several failure zones. chapter 2 correlated these zones with differences in failure modes, identified by their plastic strain fields. However the methodology of Rice's criterion can help quantify and distinguish failure mechanisms.

Let n be a unit vector such that det n L n = 0, and call g the unit eigenvector of n L n corresponding to eigenvalue 0. Rice's criterion then predicts a localization band of normal n , and the two sides of the band have a relative velocity along g . Two extremal cases can be distinguished. The case where g and n are collinear represents a perfect extension mode: the band opens without tangential movement. If g and n are orthogonal, the two half spaces delimited by the band have a pure tangential displacement: this represents pure shear. Therefore the analysis Furthermore the value of g • n distinguishes extension and shear failure modes. This analysis can be performed not only at localization onset, but more generally throughout the computation.

To this end, the vector n minimizing det n L n is found, and the eigenvector for the lowest eigenvalue of the associated acoustic tensor is computed. The solution (n , g ) only describes a localization band at the instant of failure. Nonetheless, this provides insight on likely bands before localization, and the evolution of such bands after localization. Fig. 3.4, 3.5, 3.7 and 3.8 depict the characteristics of the localization band for the unit cell and the random cell R1 for two loading conditions (T = 1, L = -1) and (T = 1, L = 0). Each figure shows first the plastic strain rate field at the instant of localization. Plastic strain rate is normalized by the strain rate Ė11 , which is required to drive the computations but plays no role for a rate-independent plastic material. The evolution of g • n and of the principal stress σ 11 is then represented for the whole simulation. Finally the evolution of the det n L n landscape is shown at several time steps during the simulation. The minimizing value is shown with a blue triangle (resp. red dot) before (resp. after) localization. White areas correspond to negative values of det n L n . The unit vector n is parametrized as:

n = (sin(φ ) cos(θ ), sin(φ ) sin(θ ), cos(φ )) (3.26)
For the unit cell loaded with (T = 1, L = -1), localization takes place near E loc ∼ 0.25. The localization band is a plane orthogonal to the first axis, as shown by the plastic deformation rate field. This is consistent with the evolution of the det n L n landscape. Initially, while the cell is still completely elastic, det n L n does not depend on n . As soon as plasticity begins, a minimum is found near (but not exactly) θ = 0 [π], φ = π/2, i.e. n = e 1 . The quantity g • n increases until localization and stabilizes shortly after at a value close to 1. The failure mode is thus extension, which is consistent with a localization band orthogonal to the tensile axis. On the contrary, for (T = 1, L = 0), the localization takes place immediately after yield. Two equivalent localization bands at 45 • from the cube's first axis are present. These bands can also be seen on the det n L n landscapes. At the localization onset, g • n ∼ 0.1: this low value is characteristic of the shear failure mode. However, during the rest of the simulation, g • n continues to grow, and the local minima in the det n L n come closer one to another. They even merge for E 11 ∼ 0.3, which also correspond to a stabilization of g • n at a value close to 1. This can be interpreted as a transition from an initial strain localization in shear mode to an extension mode (Fig. 3.6).

For the random cell R1 with T = 1, L = -1, results are globally similar to the situation of the unit cell for the same loading condition. However the localization band has a more complex pattern, as it tries to connect several voids. The overall orientation of this band remains orthogonal to the first loading axis. For T = 1, L = 0, localization takes place early, in shear mode. Contrary to the case of the unit cell, no stabilization of g • n and no transition to extension mode are witnessed. Moreover, due to the lower symmetry of the random microstructures, the det n L n landscape has lost its symmetry with respect to φ = π/2. The oscillations on Fig. 3.8b correspond to an alternation between two local minima.

The link between localization mode and failure strain is shown in Fig. 3.9. For all cells, the zones of reduced ductility identified in Fig. 3.3c are systematically associated with shear failure mode. On the other hand, the peaks and plateaus correspond to extension modes. These results are consistent with those obtained previously with a coalescence criterion. However, in the previous study, modes were only distinguished by the inspection of the plastic deformation field. On the contrary, the present study's analysis of the acoustic tensor and its eigenvectors allows to distinguish quantitatively the two localization modes. Moreover this analysis is straightforward, as the acoustic tensor has already been computed in order to apply Rice's criterion.

Influence of the loading orientation with respect to the cube

In the previous section, Rice's criterion was found to effectively predict failure by localization and was able to distinguish between extension and shear modes. Only simple loading conditions with a diagonal η were considered. The evolution of the strain at localization with respect to L in Fig. 3.3 showed plateaus and cusps due to the hesitation between failure modes. It is questionable whether this evolution still holds for different loading orientations with respect to the cube. General loading conditions with an arbitrary η ∼ are thus now considered. They allow investigating the influence of the loading orientation. The effects of the intrinsic anisotropy of the simulation cells can then be explored.

Parametrisation of the general loading orientation

As shown in Eq. (3.14), the tensor η ∼ , defining the average Cauchy stress tensor up to a proportionality constant, is equivalent to the combination of T , L, and a rotation matrix Q. This rotation represents the orientation of the loading with respect to the cube. In order to define the computations and report on the results, a suitable parametrisation for Q must be adopted. As the rotation group SO 3 is three-dimensional, Q can be represented by three parameters, for instance by the three components of a Rodrigues vector. If Q corresponds to the rotation of angle θ around the direction given by the unit vector m , with ||m || = 1, the associated vector is:

d = tan(θ /2)m (3.27)
In order to explore all possible loading orientations, the whole rotation group SO 3 is not necessary due to symmetry conditions. The periodic boundary conditions naturally entail cubic symmetry. The unit cell exactly displays this symmetry so no generality is lost when reducing loading orientations by cubic symmetry. In order to limit the computational expense, we suppose that it is also the case for random microstructures. Equivalently, the image of a random microstructure by a symmetry of the cube can be thought as a distinct microstructure, whose behavior is not studied in this work. On the other hand, the symmetric tensor η ∼ is characterized by three orthogonal axes bearing distinct principal stresses. It can thus be associated with orthotropic symmetry. Accounting for symmetries, the set of all possible loading orientations with respect to the simulation cell is the cubic-orthotropic disorientation space (using [START_REF] Heinz | Representation of orientation and disorientation data for cubic, hexagonal, tetragonal and orthorhombic crystals[END_REF] terminology). Inversion of a principal axis does not change η ∼ , so the disorientation space is only one half of the one found by [START_REF] Heinz | Representation of orientation and disorientation data for cubic, hexagonal, tetragonal and orthorhombic crystals[END_REF]. The appropriate disorientation space can be expressed in Rodrigues formalism:

0 ≤ d i ≤ √ 2 -1 for i = 1, 2, 3 d 1 + d 2 + d 3 ≤ 1 (3.28)
A sample of loading orientations in this space is represented in Fig. 3.10. They are repeatedly used in the following section. Red points correspond to the vertices of the disorientation space, whereas green points were randomly sampled. The latter were first sampled on the unit quaternion sphere with a uniform probability distribution, then were converted to Rodrigues vectors. The representative in the disorientation space is then computed with the method of [START_REF] Grimmer | Disorientations and coincidence rotations for cubic lattices[END_REF].

From T , L and a Rodrigues vector d , the macroscopic stress tensor η ∼ can be computed with the methodology described in subsection 3.3.3. However the 1-axis plays a special role in the spring element applying boundary conditions (Eq. (3.12)). The computations were found to be more stable when η 11 is the maximal component on the diagonal of η ∼ . The principal axes of η ∼ are thus cyclically permuted until this condition is reached. This operation is compatible with cubic symmetry and does not change the equivalence class of loading orientations.

The component E 11 does not take into account rotations of the loading axes, and cannot be used to compare simulations in different loading orientations. By extension of the diagonal case, the following definition of the failure strain is proposed. At the onset of localization t loc (where t is the fictitious time driving the simulation), the macroscopic deformation gradient F ∼ (t loc ) is extracted. If e I is the normalized eigenvector of η ∼ , associated with its largest eigenvalue, E loc can be defined as:

E loc = ( F ∼ (t loc ) -1 ∼ ) : (e I ⊗ e I ) (3.29)
This is equivalent to rotating the deformation gradient in the principal frame of the applied stress, and extracting the new component E 11 :

E loc = Q ∼ T • ( F ∼ (t loc ) -1 ∼ ) • Q ∼ 11
(3.30)

Figure 3.10: Disorientation space representing all loading orientations with respect to the cube, along with simulation points frequently used in this study where the notation of Eq. (3.14) is used. In the case of diagonal loading, Q is the identity, and the previous definition of E loc is recovered.

Results for general loading orientations

With the preceding formalism, simulations with any loading orientation can be performed. Before exploring the full orientation space, the simpler case of a rotation with respect to one axis of the cube is studied (Fig. 3.11). For the condition (T = 1, L = -1), simulations are performed on the unit cell and two random cells. The applied stress is rotated with various angles around the second or third axes of the cubes. In the Rodrigues formalism, d is collinear with e 2 or e 3 . As L = -1 is an axisymmetric case, rotation around e 1 has no effect on the applied stress. For the unit cells, results for the second and third axes are identical, due to symmetry. Failure strain decreases as the angle of rotation increases, until a minimum is reached at approximately 35 • . Notwithstanding the angular discretization, this is fully consistent with the results by [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], who found a minimum for an angle of 37 • . On the other hand, the two random microstructures show a more complex dependence of the failure strain with respect to the angle of rotation. The symmetry between the axes 2 and 3 is lost with the random cell. Moreover, results from R1 and R2 are quite different one from another. Minimum failure strain is not necessarily found for an angle of 37 • (10 • for R1, for instance). A clear peak is found for R1 and R2 for a rotation of 40 • around e 2 . The analysis of plastic deformation fields (not presented here) shows a transition of modes between the two sides of the peak. This phenomenon is similar to the cusps in Fig. 3.3c.

In order to explore the effect of orientation for given values of T and L, 150 different rotations were chosen in the disorientation space. Failure strain values were computed for these rotation points on the unit cell and the random cells at (T = 1, L = -1). Results are displayed in the disorientation space (Fig. 3.12). The case (T = 1, L = 0.5) is also shown for the unit cell. In the rest of this section, a coarser mesh with 7682 nodes (the same as in chapter 2) was chosen for the unit cell in order to limit the computational expense. For all cases, failure strain strongly depends on the loading orientation. Note that the minimal failure strain is systematically lower than in Fig. 3.11. In order to find this minimal failure strain, restricting to rotations around the axes is therefore not justified. Values of failure strain appear to depend in a rather smooth way on the loading orientation: Peaks corresponding to hesitation between modes as in Fig. 3.11a are scarce. Very clear trends can be seen for the unit cell. For instance, at (T = 1, L = 0.5), E loc appears to mainly depend on the value of d 2 . However results from different cells or different values of L do not display the same trends with respect to the loading orientation.

After having understood the structure of the dependence with respect to Q for a given value of L, we consider the loading condition as mainly a function of (T, L) and study the scatter which can be attributed to the effect of the loading orientation. Computations were performed on the three cells for various values of L, at fixed T = 1. For each value of L, 30 simulations with different loading orientations (visible on Fig. 3.10) were performed. Results are shown in Fig. 3.13 as box plots, representing the minimum, first, second (median), third quartiles and maximum for the 30 different orientations. For each cell, loading orientation is responsible for a large scatter. An order of magnitude can sometimes be found between the minimal and maximum failure strains for a given L. Results from the unit cell are even more scattered, which is the sign of greater anisotropy for the localization behavior. For L = 0, this scatter is however minimal, because the behavior of the material imposes early localization in shear. For random cells, minimal failure strain is a symmetric function of L, with minimal value for L = 0. The curve of the minimal failure strain is however flatter for the unit cells, due to the higher anisotropy. Counter-intuitively, localization is thus found to take place earlier in the unit cell than in the random cells, because the former is more sensitive to the loading orientation. The localization mode can also be investigated depending on the loading orientation. Fig. 3.14 shows the plastic deformation rate fields just after localization for the random cells R1 (Fig. 3.14a to 3.14d) and R2 (Fig. 3.14e to 3.14h). Two loading cases are considered: (T = 1, L = -1) (first and third rows) and (T = 1, L = 0.4) (second and fourth rows). For each loading case, the results for the most (resp. least) resistant orientations corresponding to the maximal (resp. minimal) failure strain value are shown on the left (resp right). The failure mode for the least resistant orientation is always associated with a shear failure mode. This is even true for the cell R2 at L = -1, for which the localization band is parallel to the sides of the cube. For the most resistant, shear or extension modes are possible. This shows that L is not sufficient to distinguish between the two failure modes. For a given L, the failure mode can change depending on the loading orientation.

The link between early localization and shear failure mode can be quantified with g • n , as in subsection 3.4.2. Fig. 3.15 correlates g • n at localization and E loc for all the simulations on random cell R2, seen in Fig. 3.13b. E loc is almost a linear function of g • n . A small failure strain is thus systematically linked to shear failure mode, whereas late localization is only found with extension mode.

The results of this section show the strong influence of the loading orientation. Yet the response of an isotropic material should only depend on the stress invariants T and L. The dependence with the loading orientation is therefore the sign of the anisotropic behavior of the cubic simulation cells. This anisotropy cannot be avoided but could be reduced with larger simulation cells. The random microstructures already display a smaller anisotropy than the unit cell. As proposed by [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF] and in order to obtain conservative results, the application of simulations to isotropic models should focus on the critical loading orientation with the minimal failure strain. This critical orientation may vary with T and L, but always corresponds to a shear failure mode.

Discussion

The previous results are now discussed with help from complementary simulations. Rice's criterion is first compared to other criteria from literature. Then the link between cell anisotropy and the scatter of results is further analyzed. The anisotropy in yield strength is studied, and the influence of a reduced porosity on localization is investigated.

Comparison of failure criteria for diagonal loading

In this study, Rice's localization criterion was used to detect the failure of the cells. Its predictions in the case of diagonal loading are here compared to those of two coalescence criteria (Fig. 3.16). The energy criterion from [START_REF] Wong | On the energetics of tensile and shear void coalescences[END_REF] compares the macroscopic plastic power Ẇp and elastic power Ẇe . In the coalescence state, elastic unloading occurs far from the voids, so the ratio Ẇe / Ẇp is minimal. The δ criterion assumes that for shear or extension coalescence modes, the macroscopic deformation gradient rate Ḟ ∼ becomes singular (rank one for extension, and rank two for shear). The evolution of det Ḟ ∼ can then be compared to that of a homogeneous loading. Coalescence is thus detected when this ratio reaches zero (within a given precision). For both the unit and random cells, localization (from Rice's criterion) takes place earlier than or almost at the same time as coalescence (from the two other criteria). This is consistent with the findings by [START_REF] Guo | Void-sheet analysis on macroscopic strain localization and void coalescence[END_REF] and [START_REF] Zhu | Investigation of the competition between void coalescence and macroscopic strain localization using the periodic homogenization multiscale scheme[END_REF]. For the unit cell, the three criteria give very similar predictions for L < -0.5 and L > 0.3. However, in the region L ∈ [-0.5, 0.3], localization takes place very early, whereas failure strain from the energy criterion remains high and grows with L. The behavior of the δ criterion is intermediate, with results close to that of the energy criterion, except in the immediate vicinity of L = 0, where failure strain is minimal. This division of the E loc -L curve in three zones matches the results from Fig. 3.9. The three criteria give similar results when the unit cell fails in extension, and quite different results for the shear failure mode. The case L = 0 can be better understood from Fig. 3.6. There is early strain localization in shear, but the failure mode ends switching to extension mode, which is recognized by the energy criterion. For the random cell, at low and high L, coalescence criteria give the same results. Coalescence and localization are almost simultaneous at high L, whereas localization occurs significantly earlier at low L (for generalized tension). In the region L ∈ [-0.5, 0.4], the energy criterion fails to detect any coalescence (effect already seen in chapter 2). The failure strain from the Rice and δ criteria both form a U-shaped curve. The curve is wider for localization, as localization takes place earlier than coalescence. In this region, failure takes place in shear mode and no switch to extension mode was seen in Fig. 3.8. To summarize, localization always precedes coalescence. The three criteria yield similar results for extension modes. For shear modes, the energy criterion fails to detect coalescence whereas the other criteria can recognize this mode. The zone of low ductility is wider for the localization criterion.

Anisotropy of the cells

The results of subsection 3.5.2 showed a major dependence of the failure strain on the loading orientation with respect to the cell axes. This can be interpreted as an evidence of an anisotropic localization behavior of the simulation cells. In order to precise this anisotropy interpretation, we show that the dependence on loading orientation does not only concern localization strain, but also the initial yield stress. The simulations from Fig. 3.13 are analyzed again, and the maximum value of the principal stress σI is extracted. Values of σI can be directly compared between simulations at different loading orientations. Fig. 3.17 shows the scattering of stress values with respect to the orientation for the unit cells and the random cells R1 and R2. This scattering is significant but is not as marked as for failure strain values. It confirms the anisotropic behavior of the simulation cells. Compared to unit cells, anisotropy is significantly reduced for the random microstructures, in agreement with Fig. 3.13. However, the scatter of maximum stress is maximal near L = 0, whereas the scatter of localization strain was minimal there. The anisotropic behavior, both for localization and for initial yield stress, is due to the geometric anisotropy of the cells. With the periodic boundary conditions, the material is a cubic array of voids and therefore displays cubic anisotropy. The anisotropy is smaller for random microstructures because they correspond to larger volume elements, or equivalently because the period of the cubic array created by boundary conditions is larger (relatively to the size of the voids).

Cubic unit cells thus exhibit major anisotropy. In order to determine an isotropic response from cell simulations, computations should be performed for many loading orientations. For random microstructures, this anisotropy is smaller but remains significant. With larger microstructures containing more voids, the anisotropy could probably be further reduced, at the expense of a higher computational cost.

Another method to reduce anisotropy would be to use different boundary conditions, for these conditions play a key role in localization. For instance, [START_REF] Daehli | Ductile failure predictions using micromechanically-based computational models[END_REF] compared several unit cells with different boundary conditions and found very distinct onsets of localization. However, in micromechanical analyses, there are no realistic boundary conditions especially in the presence of instabilities. More pragmatically, boundary conditions should play as small (or at least as local) a role as possible on the behavior of the representative volume element and should not hinder too much global failure modes. Periodic boundary conditions are satisfactory for strain localization analyses because they allow the development of (compatible) localization bands, and are found (in elasticity) to allow faster convergence of results with cell size [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]. Parallel faces boundary conditions for random and unit cells were tried in subsection 2.B.2 but were found to strongly hinder the development of shear failure modes. [START_REF] Coenen | Novel boundary conditions for strain localization analyses in microstructural volume elements[END_REF] also proposed a new type of boundary conditions to counteract the influence of the periodic boundary conditions, which are here responsible for the anisotropy. The most likely localization band is detected from the strain field. The boundary conditions are then adjusted to allow the development of this band (and only this one). Yet this approach does not seem adapted to the complex localization patterns seen in Fig. 3.14 or to possible transitions of failure modes (Fig. 3.6).

-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 L

Effect of porosity

A single fixed porosity of 6% was considered up to now. In order to assess the influence of porosity, a unit cell of porosity 1% was generated with the same procedure as previously. The only difference resides in the radius of the embedded spherical void. This cell was subjected to the same loading conditions as in Fig. 3.13 and the results are presented in Fig. 3.18. When compared to the 6%-porosity unit cell, a significant influence of the loading orientation is still present. The minimal failure strain still draws a U-shaped curve, with a minimum near L = 0. However the values of failure strain are generally higher. This is logical as lower porosity induces higher resistance to ductile fracture. More interestingly, the relative scatter with orientation is reduced, and the curve of minimal failure strain is less flat than for the 6%-porosity cell. As the void is less susceptible to interact with the periodic boundary conditions, the anisotropy is reduced for low volume fractions.

Conclusion

Random microstructures consisting in a perfectly plastic matrix containing a random distribution of voids were simulated up to failure. A wide range of loading conditions, characterized by the stress triaxiality and the Lode parameter, was studied. The orientation of the principal directions of the applied stress tensor with respect to the cubic cell's axes was systematically investigated.

Failure was detected with a strain localization indicator, namely Rice's criterion computed at the macroscopic level. The major findings of this study follow:

• Rice's criterion successfully detects the onset of strain localization and the orientation of the main localization band for both unit cell computations and random distribution of voids. Moreover it allows distinguishing two types of failure modes: shear and extension modes.

• For diagonal loading, in which the directions of the applied stress coincide with those of the cubic cell, a central zone of reduced ductility is found both for the unit cell and the random cells. Shear failure mode is systematically associated with reduced ductility.

• The loading orientation strongly affects failure strain. For given stress invariants, the values of failure strain for the least and the most resistant orientation can differ by more than one order of magnitude. The failure mode may depend on the loading orientation.

• The sensitivity to loading orientation is due to the intrinsic anisotropy of the cell. Although it affects both random and unit cells, the scatter is reduced for random microstructures, as they are less anisotropic.

• Minimal failure strain on all loading orientations was obtained both for unit cells and random populations of voids. It is a U-shaped function of L, minimal for L = 0. In order to reach this minimum, considering only a rotation along a fixed axis of the cube, as commonly done in literature, is not sufficient and general 3D rotations are required.

In order to develop reliable models of ductile fracture, minimal failure strain should be sought among all loading orientations. Moreover, random microstructures should be favored over unit cells because of their lower anisotropy in localization and inclusion of void interaction. However the simulation of random microstructures is significantly more computationally expensive than unit cells. An extension of this work should be to obtain a surrogate model for the ductile fracture of random cells. This model should require as small a number of simulations as possible by leveraging the link between random microstructures and unit cells.

3.A Validation of the computation of macroscopic tangent operator

The method for computing the macroscopic tangent operator is here validated by comparison with other methods. Homogeneous cubes of 3 × 3 × 3 elements without voids are first considered. A purely elastic material and an elastoplastic material are considered. The formulation of Eq. (3.3) is used, but R 0 is replaced either by R elas = +∞ or R hardening (p) = R 0 + H p (with H = 2000 MPa). Periodic boundary conditions and a large strain formulation are used. A diagonal loading condition corresponding to (T = 1, L = -1) is applied. As the cubes are homogeneous, the macroscopic tangent operator is equal to the local tangent operator, which is directly computed by the software (Zset, 2022) within the element. Fig. 3.19 compares the evolution of the Rice criterion using macroscopic tangent operators obtained from three methods: (i) the condensation method of subsection 3.3.4; (ii) the local operator within an element; (iii) taking K EE (in Eq. (3.23)) as a macroscopic tangent operator. This last method is equivalent to averaging the local tangent operator over all elements, and neglecting their interactions. The three methods give identical results for the two types of materials, as should be. However, the averaging method only works here because of the homogeneous cube. For the elastic material, the Rice criterion is not constant because of the large strain formulation. A more complex case is then considered with an elastic porous unit cell. Strain, and therefore the local tangent operator, is no more homogeneous within the cell. The reference results are now computed with a perturbation method (Zhu, Bettaieb, et al., 2020). Regularly during the main simulation, the state of the finite element problem is saved. Secondary computations are then launched from the saved states by applying a small perturbation ∆ F ∼ = e i ⊗ e j , with i, j ∈ [ [START_REF] Keralavarma | Ductile failure as a constitutive instability in porous plastic solids[END_REF]3]]. The evolution of the Boussinesq stress tensor ∆ S ∼ is then observed. The macroscopic tangent operator can then be estimated by Li jkl = ∆ Si j /∆ Fkl . This method is less practical than the condensation method because it requires launching nine secondary computations for each output time of the main computation. More problematically, it is not adapted to multibranch tangent operators, which is automatically the case in elastoplasticity. Perturbations ∆ F ∼ = ±e i ⊗ e j may correspond to either plastic loading or elastic unloading. Distinguishing which one corresponds to the correct evolution in the primary computation is not easily possible a priori. This problem does not affect pure elasticity. For (T = 1, L = -1) on the elastic unit cell, the condensation and perturbation methods give identical results ( Fig. 3.20). However only taking the average of the local operators is not enough, for it implies neglecting structural effects. (with Q = 1 MPa and otherwise the same characteristics as in subsection 3.3.2). For a von Mises material, according to [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF], the localization should take place when dR/d p = -E /4. This is effectively the case, as shown by Fig. 3.21.

3.B Sensitivity of Rice's criterion to simulation parameters

This section studies the influence of the temporal and spatial discretization on the localization results (Fig. 3.22). Simulations were carried out with diagonal loading conditions, at T = 1 and various values of L. The fine mesh from Fig. 3.1c and the coarse mesh from Fig. 2.1 were compared. The effect of fine time steps was also studied. As the simulation is driven by the evolution of E 11 , time steps correspond to the increase ∆E 11 between two iterations of the implicit Euler method. In the main study, increasing time steps from ∆E 11 = 0.0001 to ∆E 11 = 0.005 were used. This allowed a sufficient precision even in the case of early localization. and finely meshed unit cells, the overall appearance of the curves are similar with a plateau at low L and a zone of reduced ductility around L = 0. The behavior at high L > 0.3 is however different with a clearer plateau for the fine mesh. This justifies preferring the finely meshed cell in this study. Nevertheless, the unit mesh was used for Fig. 3.12 and 3.13 in order to limit the computational cost. For the fine cell, there is almost no effect of the time step (except at the lower end of the central reduced zone). The effect is more pronounced for the coarse cell, especially in the L > 0.3 region. The finer the mesh and the time steps, the narrower the zone of reduced ductility.

3.C Loss of ellipticity in the matrix

This section elucidates the role of the matrix in the macroscopic loss of ellipticity. For a von Mises plasticity and especially in a large strain formulation, a loss of ellipticity in the matrix for L 0 cannot be a priori excluded.

To investigate this matter, the same computations (diagonal loading with T = 1 and L = -1) as in Fig. 3.3c were performed on a homogeneous perfectly plastic cube of size 5 × 5 × 5 elements without any void, with the parameters of the matrix material. During the simulations, deformation remains homogeneous until the onset of localization. The localization analysis was performed for these homogeneous cubes, in the same way as for the voided cells. The results are shown on Fig. 3.23 superimposing results previously obtained for unit and random cells. For |L| > 0.2, the simulation were performed until a strain of E 11 = 1 but no localization was observed. The reported values of E loc = 1 represent a lower bound of possible later localization (which might never take place); they are represented by upwards-facing triangles. For |L| ≤ 0.2, localization was found in a homogeneous cube but significantly later than for voided cells (except for L = 0 where localization is so early that it is simultaneous for all cell types). In all cases but L = 0, the role of the matrix's loss of ellipticity appears negligible compared to the roles of the voids, or at least cannot be captured by a homogeneous defect-free cell. Moreover the conservation of ellipticity in the localization bands between the voids was checked, because localization can lead to mesh size dependency. Using the methodology developed by (Al [START_REF] Kotob | A general and efficient multistart algorithm for the detection of loss of ellipticity in elastoplastic structures[END_REF], a local Rice analysis was performed at each Gauss point of the matrix material for several simulations of the unit cell under diagonal loading conditions (more precisely, for T = 1, L ∈ {-1, -0.5, 0, 0.5}). Fig. 3.24 shows the det n L n field at the global localization onset for the (T = 1, L = 0) case (which is prone to ellipticity loss); other loading cases provide similar results. Some negative values are found, but in isolated elements. This appears a numerical artifact rather than a real pattern of ellipticity loss inside a band. The matrix's loss of ellipticity can therefore be neglected. Moreover the influence of mesh size was investigated in subsection 2.B.1, which showed that similar values of strain at coalescence and localization patterns are found for finer and finer meshes of the same geometry. Mesh sensitivity therefore appears small and plays no role on the results of the present study.

3.D Localization analysis on porous balls

Cubic unit cells were shown to exhibit major anisotropy for localization. In order to determine an isotropic response from cell simulations, computations should be performed for many loading orientations. For random microstructures, this anisotropy is smaller but remains significant. With larger microstructures containing more voids, the anisotropy could probably be further reduced, at the expense of a higher computational cost. Requiring many simulations for the same loading conditions (T, L) is a hindrance to the development of models, especially for computationally expensive computations on random microstructures. Cells with lower localization anisotropy could be therefore sought for. Another method would be to modify the boundary conditions instead of the cell. Such a modification was proposed by [START_REF] Coenen | Novel boundary conditions for strain localization analyses in microstructural volume elements[END_REF] to counteract the influence of the periodic boundary conditions, which are here responsible for the anisotropy. The most likely localization band is detected from the strain field. The boundary conditions are adjusted to allow the development of this band (and only this one). This approach does not seem adapted to the complex localization patterns seen in Fig. 3.14 or to possible transitions of failure modes (Fig. 3.6). This section investigates the localization in a hollow ball with periodic boundary conditions. [START_REF] Glüge | Comparison of spherical and cubical statistical volume elements with respect to convergence, anisotropy, and localization behavior[END_REF] investigated this type of cells and found a reduced lower elasticity anisotropy when compared to cubic cells. Moreover they were able to study localization in homogeneous spherical cells, made of a damaging material.

The ball under consideration contains a single central void representing a porosity level of 6%. The structured mesh of a hemisphere is shown in Fig. 3.26. Although the ball is invariant by rotation, the mesh displays a slight geometrical anisotropy. Following [START_REF] Glüge | Comparison of spherical and cubical statistical volume elements with respect to convergence, anisotropy, and localization behavior[END_REF], periodic boundary conditions are applied exactly like on cubic cells (Eq. (3.6)). The pairs of homologous nodes are now formed by diametrically opposed nodes on the exterior surface. Nevertheless there is a mathematical difference between cubic and spherical cells. Cubic cells with periodic boundary conditions form a representation of a three-dimensional torus, whereas the ball associating diametrically opposed points represent a projective space.

With this spherical cell, the same computations as described in subsection 3.3.3 can be performed. Results from Rice's indicator for various values of L and 30 different orientations are shown in Fig. 3.25. The overall aspect of localization strain curves are reminiscent of the results with cubic cells with a minimum failure strain near L = 0. However localization seems to take place later than for cubic cells. Moreover there is a significant and spurious dependence on the loading orientation. Fig. 3.27 focuses on the simulations from 30 different loading orientations at (T = 1, L = -1). All the stress-strain curves share a common inflection point near E I = 0.2. This point should be clearly considered as the onset of failure of the cell in the simple axisymmetric case L = -1. After this failure point, the stress strain curves may diverge one from another. However the Rice criterion curves for the different simulations do not agree on the value of the localization strain. These differences appear due only to the slight geometrical anisotropy of the spherical mesh. For some loading orientations, Rice's criterion may even never detect localization. Therefore Rice's criterion seems unable to reliably detect localization in the unit spherical unit cells. This can be explained by the failure mode seen in Fig. 3.26. For this loading orientation, no localization was detected by the Rice criterion. A complex pattern of strain localization can be seen, but it does not form a band. Contrarily to a torus, a planar band is impossible in a projective space. As the Rice criterion is derived from the hypothesis of such a band, it logically fails in the case of a spherical unit cell. Similar negative results (not shown here) were obtained with the energy and δ criteria from subsection 3.6.1. Therefore hollow balls with periodic boundary conditions do not appear reliable for unit cell studies of ductile fracture, as they hinder the formation of consistent localization bands. Chapter 4

Multifidelity ductile damage model for random microstructures

A nonparametric surrogate model for ductile failure is developed from simulation results on cubic representative volume elements with a random distribution of voids. This model fully takes into account the anisotropy due to the simulation conditions. The metamodelling strategy uses Gaussian Process Regression coupled with a multifidelity approach involving simulations on cells with a single void. Through cokriging and metamodel parameter transfer, information can be transferred from the unit cell simulations, which are relatively inexpensive, to the model on random cells containing 27 spherical holes. This allows an increased accuracy for a given computational budget. Strategies for adaptive experimental design are also investigated.

Résumé détaillé

Les chapitres précédents ont montré l'utilité des simulations sur des microstructures contenant une distribution aléatoire de défauts, en permettant de prendre en compte l'interaction entre défauts. La résistance de telles microstructures à la rupture ductile a pu être étudiée en analysant la localisation à l'échelle globale de la cellule. La déformation à rupture E loc est alors une fonction des conditions de chargement et notamment de leur orientation. L'importance de cette orientation témoigne d'une anisotropie due aux cellules de simulations cubiques et aux conditions aux limites périodiques.

Cependant, les simulations sur microstructures aléatoires ne permettent pas seulement d'étudier les phénomènes de rupture ductile mais peuvent également fournir un modèle de résistance à la rupture ductile en fonction des conditions de chargement. Ce modèle pourrait alors être utilisé dans des simulations de pièces réelles.

Pour construire ce modèle de rupture ductile, il convient donc de déterminer la déformation à rupture pour toute condition de chargement (T, L), en utilisant les résultats d'un nombre fini de simulations. Afin de se rapprocher de l'application visée, les simulations seront réalisées avec une matrice présentant un comportment matériau plus réaliste. L'exemple académique du matériau élastique-plastique parfait utilisé dans les chapitres précédents est donc remplacé par une loi de 97 comportement d'un matériau élastoplastique avec écrouissage se rapprochant de celui de l'alliage de titane TA6V dans un cordon de soudure. La question de l'anisotropie du comportement en rupture ductile reste présente avec ce nouveau comportement (section 4.3). La méthodologie d'obtention d'un modèle de rupture ductile sera cependant illustrée en se concentrant sur une unique réalisation de cellule aléatoire et donc en ignorant la dispersion résultant de la distribution aléatoire des défauts.

Afin de reproduire l'ensemble des résultats au sein d'un modèle de rupture ductile, une approche par surface de réponse non paramétrique (métamodèle) est utilisée, ce qui permet une plus grande expressivité par rapport aux modèles paramétriques disponibles. La méthodologie d'obtention de cette surface de réponse se fonde sur la régression de processus gaussien ou krigeage. Cette approche consiste à identifier de façon bayésienne une réalisation d'un processus gaussien à partir des résultats de simulations. L'approche permet donc une interprétation fournissant des niveaux de confiance, et les paramètres qu'elle fait intervenir ont un sens physique clair. Cependant, une exploration fine de l'espace des conditions de chargement possibles est impossible à réaliser pour les microstructures aléatoires, car les cellules aléatoires impliquent des simulations longues et coûteuses numériquement. Une possibilité est donc d'utiliser les résultats sur les cellules unitaires, qui peuvent fournir des informations utiles à un coût de simulation réduit. Une méthodologie de transfert d'information entre cellules unitaires et cellules aléatoires est présentée en section 4.4. Elle s'appuie d'une part sur un transfert de paramètres entre les processus gaussiens identifiés pour les cellules unitaire et aléatoire, et d'autre part sur une réprésentation de la corrélation entre les résultats, par la méthode du cokrigeage.

Grâce à cette stratégie de transfert d'information, il est possible d'obtenir une surface de réponse représentant correctement les résultats de simulations, et prenant en compte l'anisotropie due aux conditions de simulations (section 4.5). Pour un même nombre de simulations sur cellules aléatoires, l'utilisation des informations des cellules unitaires permet une augmentation de la précision, en comparaison à une surface de réponse construite uniquement à partir des résultats de cellules aléatoires. Cette meilleure performance apparaît davantage due au transfert des paramètres de processus gaussiens qu'à l'utilisation du cokrigeage.

Enfin, la régression de processus gaussiens peut être utilisée pour décider de façon adaptative du plan d'expériences ou de simulations. Les résultats de simulations précédentes peuvent être utilisés pour déterminer les nouvelles conditions de simulations susceptibles d'apporter le maximum d'information, que ce soit pour les simulations sur cellules unitaires ou aléatoires. Cette approche de plan d'expériences bayésien est discutée en section 4.6. Si cette méthode permet d'identifier plus rapidement les tendances générales de E loc qu'une méthode naïve, elle ne semble pas permettre d'identifier plus rapidement les niveaux absolus et précis de déformation à rupture.

Introduction

By simulating the failure of unit and random cells, the preceding chapters provided useful information on the ductile fracture phenomenon. They showed the relevance of simulations on microstructures with randomly distributed voids so as to take into account interactions between defects. The resistance to ductile failure was studied by analyzing strain localization at the cell level for various loading conditions. The simulations on random microstructures can therefore also be used to identify a homogenized ductile failure model taking into account loading conditions. Such a ductile failure model could then be used in simulations of industrial parts.

Ductile failure models were built by [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] or Hure ( 2021) who calibrated yield surfaces on a limited set of loading conditions. Such yield criteria can then be used in a homogenized model coupling plasticity and damage. This model will provide the initial damaging behavior and an estimation of the strain at failure. However this failure strain, which is the true quantity of interest when determining resistance to failure, can also be extracted from the simulation. It could be used to calibrate a damage criterion representing the failure of the random cells directly.

Many parametric models for a ductile damage criterion have been proposed [START_REF] Lou | New ductile fracture criterion for prediction of fracture forming limit diagrams of sheet metals[END_REF][START_REF] Mohr | Micromechanically-motivated phenomenological Hosford-Coulomb model for predicting ductile fracture initiation at low stress triaxialities[END_REF][START_REF] Defaisse | Ductile fracture of an ultra-high strength steel under low to moderate stress triaxiality[END_REF][START_REF] Zouari | The effect of strain biaxiality on the fracture of zirconium alloy fuel cladding[END_REF], but they are often phenomenological. Using a parametric model is naturally easier but reduces the complexity of patterns that the model can express. On the other hand, data driven approaches have been used to represent in a non parametric way the ductile fracture of voided unit cells. Some methodologies include kriging [START_REF] He | An improved shear modified GTN model for ductile fracture of aluminium alloys under different stress states and its parameters identification[END_REF], model order reduction with clustering [START_REF] Shakoor | Data science for finite strain mechanical science of ductile materials[END_REF], neural networks [START_REF] Guo | Crystal plasticity based investigation of micro-void evolution under multi-axial loading conditions[END_REF], support vector machines [START_REF] Beluch | Response Surfaces in the Numerical Homogenization of Non-Linear Porous Materials[END_REF]. Kriging [START_REF] Wackernagel | Multivariate Geostatistics[END_REF], also known as Gaussian Process Regression [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], appears especially useful for response surfaces thanks to its interpretable parameters and its Bayesian formulation. This formulation allows a quantification and a propagation of uncertainty, and provides a natural model selection method by maximum marginal likelihood.

However obtaining a response surface from simulations on random microstructures is complex. For random cells, there is a strong dependence of the strain at failure on the loading orientation (for unit cells, the anisotropy is even stronger). Loading orientation should thus be explicitly considered. Therefore parametrizing all loading types and loading orientations requires a five-dimensional input space (see subsection 4.3.1). A fine sampling of loading conditions is impossible in such a high dimensional space. The anisotropy in ductile failure of the cell is actually linked to the cubic shape of the cell and its boundary conditions. Yet the ductile failure model should be isotropic for an isotropic material behavior. As discussed in chapter 3, the physical strain at failure to consider should be the minimum over all loading orientations. This minimization operation implies that the results of a computation for a given loading orientation cannot be directly used for the final isotropic model. Nevertheless, cells should be finely discretized for finite element analyses. Simulations on random cells therefore involve a high number of unknowns and are computationally expensive. Only a limited number of simulation runs may be performed without exceeding R&D computational budgets.

A solution to the data scarcity for the random cell is to use results from a similar problem, but for which data is much easier to collect. The auxiliary data provides the general trends for the metamodel, which can then be refined by a more limited number of training examples from the primary problem. This methodology of transferring information between related tasks is generally known as multifidelity metamodeling [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF][START_REF] Yang | Physics-informed CoKriging: A Gaussian-process-regression-based multifidelity method for data-model convergence[END_REF][START_REF] Mell | Multifidelity adaptive kriging metamodel based on discretization error bounds[END_REF] and used within a kriging framework. In our application, the role of the auxiliary problem can be played by simulations on the unit cell, for which the reduced number of nodes per mesh leads to faster simulations. Simulations on random cells are still needed to assess the influence of void distribution. Moreover, the metamodel obtained through kriging can provide a basis for a more efficient experimental design [START_REF] Chaloner | Bayesian Experimental Design: A Review[END_REF][START_REF] Frazier | Recent Advances in Optimization and Modeling of Contemporary Problems[END_REF]. By using the predictions from the metamodel, new input points can be chosen so as to maximize the information brought by a new simulation.

In this chapter, we therefore aim to apply the methodology of the previous chapters to detect the ductile failure of a random microstructure. The objective is to obtain a metamodel of the failure strain depending on the loading conditions. A kriging methodology is used and integrates information from simulations on the unit cell. Adaptive experimental designs are also explored. To the author's best knowledge, this methodology will provide the first nonparametric model for the fracture of random microstructures, fully accounting for the anisotropy of the cell. Moreover, the possibility and usefulness of transferring information from unit cells to random cells is evidenced.

The methodology for localization analysis in random microstructures is essentially the same as in the preceding chapter, but the differences are described in section 4.3. The main difference consists in a different material constitutive law, closer to the industrial application. The procedure to generate data for unit and random cells is then explained, and the metamodeling task is formally framed. Section 4.4 describes the methodology of Gaussian Process Regression and several information transfer strategies. These strategies are applied to the task at hand and compared in section 4.5. Section 4.6 finally investigates the possibility of an adaptive experimental design method.

The notation of the preceding chapters is extended: a vector is still noted v ; tensors and matrices are noted A ∼ . An overlying bar will represent the macroscopic average of a quantity defined at the microscopic level: for instance F ∼ for the average deformation gradient. The Gaussian Process estimating a physical quantity E will be noted Ê. A star, as in E * , will indicate a minimum over all loading orientations.

Problem formulation

Methodology of micromechanical finite elements simulations

The methodology for simulating the ductile failure in random microstructures closely follows the methodology of chapter 3. Two main differences should however be noted.

Firstly, the simulations are performed on two cells: a random one and a unit one. The random cell is the cell R1 from preceding chapters, with 27 voids Fig. 4.1a. The methodology for the determination of a failure strain metamodel is thus presented for a single realization of a random cell. The dispersion due to the random distribution of voids (which was studied in chapter 2) is therefore not taken into account. The unit cell uses the coarse mesh from chapter 2 (Fig. 4.1b), with the same discretization fineness relative to the void size as for the random cell. Results can therefore be less precise than with a finer mesh as in chapter 3. However using a coarser mesh strongly reduces the number of degrees of freedom, and therefore the computational cost of a simulation on unit cell. This better fits the role of using unit cell simulations as low fidelity computations, providing approximate information at a negligible computational cost. Secondly, it should be possible to use the metamodel identified on random microstructures in the context of the industrial application to weld joints. Contrary to the preceding chapters, a material model derived from Sarre ( 2018) is therefore used. It represents the behavior of the titanium alloy TA6V in a weld joint (see next chapter for more information). The same overall formalism as in preceding chapters is used, but an elastoplastic behavior with saturating exponential isotropic hardening is used (with the same notations as Eq. ( 3.3), and expressing strain and stress in the Jaumann frame):

ė ∼ = ė ∼ e + ė ∼ p e ∼ e = 1 + ν E s ∼ - ν E (tr s ∼ )1 ∼ ė ∼ p = ṗ ∂ f ∂ s ∼ s vm = 3 2 s ∼ dev : s ∼ dev f (s ∼ ) = s vm -R(p) ≤ 0 R(p) = R 0 + R 1 (1 -e -b 1 p ) + R 2 (1 -e -b 2 p ) (4.1)
with s ∼ dev the deviatoric part of the rotated Cauchy stress tensor s ∼ , s vm the equivalent von Mises stress and ṗ = 2 3 ė ∼ p : ė ∼ p playing the role of the plastic multiplier. The elastoplastic coefficients are summarized in Table 4.1.

The model from Sarre ( 2018) is actually more complex, with viscous and damaging behaviors, which are here not considered. The viscosity is not used in the present chapter because Rice's criterion is not applicable to rate-dependent plasticity. Moreover the damaging behavior would imply local loss of ellipticity (and not only global loss of ellipticity). The results from simulations would then be strongly sensitive to discretization level. Regularization approaches are possible [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF] but were not used here due to their complexity. The strain at localization results from different values of L, still considering T = 1 and 30 different loading orientations, are shown in Fig. 4.3. For each value of L the dispersion of results due to loading anisotropy is still present. As in chapter 3, dispersion and thus anisotropy are reduced for the random cell compared to the unit cell. However, as noted before [START_REF] Barsoum | Micromechanical analysis on the influence of the Lode parameter on void growth and coalescence[END_REF], a model for the ductile failure of a material with an isotropic material should depend only on T and L, and not on the loading orientations. The anisotropy is due to the computation conditions on the unit cell, and only the minimum of E loc over all directions is a quantity of interest. Note that this minimum has a similar trend for the unit cell and the random cell. It is decreasing with L for L < 0, and increasing for L > 0, with thus a minimum for L = 0.

A similar trend was found for a non hardening material.

-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 L

Training data and objectives

We now describe the strategy for generating the data used to train the ductile failure metamodel.

For practicality, all this training data is generated in advance, before any Gaussian Process regression. A precise formulation of the objective of the metamodel is also given. The Scipy implementation of Sobol sequences with Owen scrambling [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF][START_REF] Sobol | The distribution of points in a cube and the approximate evaluation of integrals[END_REF][START_REF] Owen | Scrambling Sobol' and Niederreiter-Xing Points[END_REF] is used to generate 10000 input points corresponding to different loading conditions, in a five-dimensional hyperparallelepipedic input space (T, L, d ), where d is the Rodrigues vector corresponding to a rotation matrix Q ∼ (Eq. (3.27)). L is chosen in [-1, 1], which corresponds to all possible values. T is chosen in [0.7, 3.0], which corresponds to rather high values of stress triaxiality. However a logarithmic scaling for T is used to favor lower triaxiality values for which there is significant deformation before localization. Each component of d is chosen in [0, √ 2 -1] (Eq. (3.28)). The points with d 1 + d 2 + d 3 > 1 are eliminated (Eq. (3.28)), which leaves 9649 input points. For each of these input conditions, a FEA simulation is launched for the unit cell according to the methodology from subsection 4.3.1, and the strain at localization E loc is computed. A typical computation on a unit cell requires a CPU time of about 10 3 s.

The same procedure is applied to the random cell. 475 loading conditions (from 500 Sobol points) are used. Due to larger meshes and a higher number of degrees of freedom, a typical computation usually requires a CPU time of about 10 5 s. Simulations on random cells are thus significantly more computationally expensive than those on unit cells (even if we consider there are twenty times more of the latter).

All the training examples used, with associated E loc , are shown in Fig. 4.4. A strong dispersion of strain values, even for close values of T and L, and thus an anisotropy, is visible, especially for lower values of T . Moreover, the number of examples for the random cell is clearly limited compared to the data for unit cells. However, the representation in two dimensions underestimates the data scarcity in five dimensions. If data was organized on a regular grid, there would actually be only about three examples per dimension for the random cell, and six for the unit cell. Note that the data for T = 1 and 30 different loading orientations, presented in Fig. 4.3, are not used for metamodeling and will only be used to evaluate the resulting model. The metamodeling task is thus defined as follows. As a minimal failure strain over all directions should be taken, the quantity of interest to model is:

E R * loc (T, L) = min Q ∼ ∈SO 3 E R loc (T, L, Q ∼ ) (4.2)
The available data consists in strain at failures for the random cell (R)

E R loc (T i , L i , Q ∼ i
), and for the

unit cell (U) E U loc (T i , L i , Q ∼ i ), where (T i , L i , Q ∼ i
) or equivalently (T i , L i , d i ) corresponds to some loading condition. This task is challenging because data (and especially for the random cell) is scarce in a high-dimensional input space. Moreover, a datum E R loc (T i , L i ,

Q ∼ i
) does not directly lead to a datum E R * loc (T i , L i ) since a minimization is required. Data from unit cells can be used but there is good correlation only for E * loc values and not directly for E loc , as Fig. 4.2 showed different patterns for the two types of cells.

A coarse estimation of E R * loc or E U * loc is shown in Fig. 4.5. The T, L space is partitioned in a grid, and each cell displays the minimal failure strain over all input points with T, L in this cell. The failure strain appears to be decreasing with T , and minimal for L 0 at given T for both types of cells. However the results remain noisy and difficult to interpret. A smooth metamodel should be preferred.

Multifidelity metamodeling methodology

In order the solve the ductile failure metamodeling task with the aforementioned training dataset, a multifidelity kriging approach is proposed. Generalities about Gaussian Process Regression

Generalities on Gaussian Process Regression and cokriging

Kriging or Gaussian Process Regression is a metamodeling methodology originally developed to estimate mineral deposits [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF]. It was mathematically formalized by [START_REF] Matheron | Traité de géostatistique appliquée[END_REF] and became one of the fundamental tools of geostatistics, before being also used in machine learning contexts. We here recall the general principles [START_REF] Wackernagel | Multivariate Geostatistics[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. Kriging is naturally expressed as a Bayesian approach. The core idea is to use Gaussian processes as the prior. Those are random fields f : X × Ω → R, where X and Ω are the input and sample spaces. For any input points x 1 , x 2 , . . . x n ∈ X , the random vector ( f (x 1 ), .... f (x n )) should be Gaussian: any finite linear combination of components has a normal distribution. The distribution of a Gaussian process is summarized by its mean and covariance functions. In the following, the mean will be chosen identically zero, and the covariance function verifies, for all x 1 , x 2 ∈ X :

Cov( f (x 1 ), f (x 2 )) = K(x 1 , x 2 ) + ς 2 δ (x 1 -x 2 ) K(x 1 , x 2 ) = AK Matern 5/2 ∑ d |x 1d -x 2d | 2 λ d (4.3)
where the sum is taken on all components of x 1 and x 2 , δ is the Dirac distribution and K Matern 5/2 is a Matérn kernel with length scales λ d (one for each component) [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]:

K Matern 5/2 (r) = 1 + √ 5r + 5 3 r 2 exp - √ 5r (4.4)
This kernel is stationary: only the distance between input points is relevant, not their positions. The parameters λ and A specify how values of the process at points close one to another are correlated: λ describes the typical correlation lengths and A the amplitude of the Gaussian process. Note that a Matérn kernel produces continuous realizations of the Gaussian process. Finally ς represents noise or epistemic uncertainty, i.e. any phenomenon that is ignored during modelling. If f has the above distribution, we will note f ∼ GP(λ , A, ς 2 ). An equivalent way of considering f is to decompose it :

f (x) = f (x) + ς ε(x) (4.5)
where f is a Gaussian process of covariance function K with continuous realizations and ε is a white noise of unit variance.

In the Bayesian approach, observations are obtained for several inputs ( f (x 1 ) = y 1 , f (x 2 ) = y 2 , . . . , f (x N ) = y N ) and the distribution is updated to reflect the new information. More precisely, the posterior is computed by conditioning the prior by the observations. The predictions for new input values x # i given ( f (x i ) = y i ) have the Gaussian distribution:

f (x # ) | { f (x ) = u } ∼ N (µ ; Σ ∼ ) µ = K ∼ x # x K -1 xx y Σ ∼ = K ∼ x # x # -K ∼ x # x (K ∼ xx + ς 2 I ∼N ) -1 K ∼ xx # + ς 2 I ∼N (4.6)
where I N is the N-dimensional identity matrix, the covariance matrix K ∼ x # x is defined by:

K x # x, i j = K(x # i , x j ) (4.7)
and similarly for the other matrices K ∼ . Practically, in order to use the results from the Gaussian process Regression, a Cholesky decomposition is performed on the covariance matrix K ∼ xx . The expectation of the posterior can then be computed from the observed data, providing a prediction for new inputs. The diagonal of the posterior covariance also provides confidence intervals. Finally, realizations from this posterior can be sampled. The parameters λ , A, ς can be fitted from data by maximizing the log-likelihood:

ln P(y) = 1 2 ln det K ∼ xx + ς 2 I ∼N - 1 2 y T (K ∼ x,x + ς 2 I ∼N ) -1 y - N 2 ln(2π) (4.8)
Multiple Output Gaussian Process (MOGP) Regression is a generalization of the Gaussian Process Regression for predicting several quantities of interest at the same time. A correlation between the outputs can be obtained by introducing a linear model of coregionalization [START_REF] Fricker | Multivariate Gaussian Process Emulators With Nonseparable Covariance Structures[END_REF] or cokriging [START_REF] Wackernagel | Multivariate Geostatistics[END_REF]. The random variables g 1 and g 2 representing the desired outputs are modeled as:

g 1 g 2 = L ∼ f1 f2 + ς 1 ε 1 ς 2 ε 2 (4.9)
where f1 and f2 are independent identically distributed Gaussian processes of unit amplitude:

f1 , f2 ∼ GP(1, λ , 0) (4.10)
and ε 1 and ε 2 are independent white noises of unit variance. The 2 × 2 matrix L ∼ correlates the output from the independent Gaussian processes f1 and f2 . Note also that f1 and f2 are continuous and noise-free, since this noise is directly expressed at the level of g 1 and g 2 . For given input x ∈ X , the covariance matrix of g 1 (x) and g 2 (x) is given by:

Cov ∼ (g(x)) = L ∼ L ∼ T + ς 2 1 0 0 ς 2 2 (4.11)
Therefore, with the Cholesky decomposition, it is enough to choose L ∼ lower diagonal to express any correlation between g 1 and g 2 . Using equations Eq. 4.11 and Eq. 4.3, a joint covariance matrix for g 1 and g 2 can be computed for any (x i ) ∈ X :

Cov(g α (x i ), g β (x j )).
(4.12)

With this covariance matrix, expectation, variances and samples for the posteriors of g 1 , g 2 given observations can be computed as for simple Gaussian processes (Eq. 4.6). Moreover, the cokriging parameters L ∼ , λ , ς 1 and ς 2 can be optimized by maximizing log-likelihood as in Eq. 4.8. The distribution of g will be noted MOGP(L ∼ , λ , ς 2 1 , ς 2 2 ). Linear models of coregionalization actually offer more modeling choices. For instance, there are not necessarily as many latent processes ( f1 and f2 ) as outputs; and latent processes do not necessarily have the same correlation length. However the above description will be sufficient for the following sections. For both kriging and cokriging, the implementations in GPflow [START_REF] Matthews | GPflow: A Gaussian process library using TensorFlow[END_REF] will be used.

Application of simple kriging to the ductile failure problem

Building a metamodel for E * loc from the data E loc (T i , L i , Q

∼ i
) is complex because the inputs are initially in a 5-dimensional space, but a minimization over the orientations is required in order to produce a model with a 2-dimensional input space. This section presents the general strategy for obtaining a metamodel, first without considering information transfer between the unit and random cells.

The dataset is first preprocessed. The feature T is replaced by its logarithm (in accordance with the sampling strategy). All features (T , L, d ) are normalized between 0 and 1, according to the theoretical ranges used for Sobol sampling. The output to model is chosen as the logarithm of E loc in order to more easily represent possible change of magnitude order. The output log(E loc ) is also centered by subtracting its mean. This centering is performed independently for unit cells and random cells, and only from data used for training (not necessarily from all available data).

The general methodology is then composed of three main steps. Firstly a simple kriging is performed in 5D. From the simulation data E R (T i , L i , d i ) for the random cell, a Gaussian Process

ÊR ∼ GP(A R , λ R , ς 2 R ) (4.13)
is trained (in order to simplify notation, the subscript loc is dropped from failure strain quantities).

During training, the amplitude A R , the correlation lengths λ R , the noise amplitude ς R should be fitted. However, in standard kriging approaches, the Gaussian Process can be made interpolatory between observations by setting ς R to a small but nonzero value. A negligible noise allowing to keep an acceptable condition number for the covariance matrix can be reached with ς R = 10 -4 . The fitting of parameters is performed by maximum likelihood and a direct minimizer from Scipy [START_REF] Virtanen | SciPy 1.0: fundamental algorithms for scientific computing in Python[END_REF] is used. Secondly, from the 5D posterior for the random cell, a minimization is performed over the three orientation features. More precisely, N T L = 64 samples (T i , L i ) of stress invariants and N d = 128 orientation samples d j are created by Sobol sampling, so as to build a 5D evaluation grid ((T i , L i , d j )) i j . Several (N real = 100) realizations of the posterior are sampled. Each realization r is evaluated on the grid: ÊRr (T i , L i , d j ). The hyperparameters were verified high enough not to affect significantly the determination of the minimum. An estimation of this minimum over

Q ∼ j is then computed: E * Rr (T i , L i ) = min j ÊRr (T i , L i , d j ). (4.14)
Then the average prediction can be extracted:

E * R (T i , L i ) = ÊRr (T i , L i ) r . (4.15)
(with • r denoting the average over all realizations r). A local indicator of the dispersion of the predictions can also be found by using the empirical variance Var ÊR (T i , L i ) over all realizations.

To simplify and obtain a single dispersion indicator, the average variance will be used:

V R = Var ÊR (T i , L i ) T i ,L i . (4.16)
The last step is to build a 2D kriging. The input data for this new kriging are the minimum failure strains E * T (T i , L i ) over all orientations (or similarly for the unit cell). The prior is chosen as:

Ê * R ∼ GP(A, λ ,V R ) (4.17)
where V R is the previously defined indicator for the dispersion of minimum. The amplitude of the Gaussian process and two length scales should then be fitted by maximum likelihood. The posterior of the Gaussian process finally yields a metamodel for E * R as a function of the stress invariants T , L.

Although the metamodel for the random microstructure is our main objective, the same process can be performed separately for the unit cell (all notations are conserved, replacing the subscript R by U). However the first kriging operation ÊU ∼ GP(A U , λ U , ς U ) poses a specific problem. As there are thousands of training examples for the unit cell, the Cholesky decomposition of the covariance matrix can become expensive. The unit cell kriging will therefore use a Sparse Variational Gaussian Process (SVGP), an extension of standard Gaussian processes better suited to handle large datasets (this extension is presented briefly in appendix 4.A).

Information transfer strategies

The previous subsection described a Gaussian process regression methodology for each dataset. Kriging operations for random and unit cells were completely separate and independent: The unit cell was not used at all for the model on random microstructures. We now consider several ways of coupling these regressions and therefore to transfer information between the unit and random cells: first by transferring parameters for the 5D regressions, then by performing cokriging at the T, L level.

The transfer of parameters at the 5D level is first considered: Parameters like the length scales of the amplitude of the Gaussian priors play a key role for kriging and can be learned from data by maximizing the likelihood. On the other hand, they can be thought to be characteristic of the metamodel to find. For the random cells, the optimal value found on a very limited number of examples may not correspond to the correct value in the high data limit. However, the unit cell may be more accurate in determining these parameters thanks to the larger dataset.

A parameter transfer strategy can then be defined as follows. The parameters for the unit cell are fitted by marginal likelihood on unit cell training examples. The noise ς U = 10 -4 is constrained so as to obtain an interpolatory Gaussian process. Then, for the random cell, some parameters are set to the values found for the unit cell and are fixed, i.e they are not fitted by maximum likelihood. The remaining parameters for the random cell are learned as usual by likelihood maximization. Their values are completely independent from those found on the unit cell. Three variants of this parameter transfer strategy will be considered.

• The strategy "No Transfer" (NT for short) corresponds to the standard uncoupled methodology described in the previous subsection. It will be used as a control.

• The strategy "Length Scale" (LS) reuses only the correlation length scales from the unit cell, i.e the distance at which there is correlation, and thus influence, between inputs. We therefore set λ R = λ U , fit A R , and keep the Gaussian process for the random cell interpolatory (ς R = 10 -4 ).

• Finally, the strategy "Gaussian Process" (GP) consists in reusing the two main characteristics of the Gaussian process E U : its length scales λ U and its amplitude A U . By setting λ R = λ U and A R = A U , we try to explain the results for the random cells as much as possible by the results from the unit cell. The rest should be explained by noise, so ς R is fitted.

Secondly and independently from parameter transfer, cokriging can be used in order to represent directly the correlation between results from unit cell and random cell. Cokriging is made at the 2D level after minimization over the orientations, because results in 2D are well correlated, but not those in 5D (see subsection 4.3.2). Cokriging is performed with the data

E * U (T i , L i ), and E * R (T i , L i ).
The prior is chosen as:

Ê * U Ê * R = MOGP a 0 a b , λ ,V U ,V R (4.18)
In accordance with Fig. 4.3, the correlation between the outputs Ê * U and Ê * R is expressed by a common multiple (coefficient a) of the first latent process. The difference between the two cells is explained by the second latent process (coefficient b). Different noise amplitudes are chosen, corresponding to the dispersion indicators for the minimum. If V R V U , the Gaussian process for the random cell is allowed not to fit closely training data, in order to better reproduce the overall trend given by the unit cell. The noise amplitudes are fixed, so only the parameters a, b and the two components of λ need training. The use of this cokriging methodology will be noted CK, whereas the use of simple kriging at the 2D level will be noted SK.

Application of the multifidelity metamodeling strategy

This section applies the previously defined metamodeling strategy to the data on the unit and random cell. The strategy is first illustrated by presenting results for the GP-SK variants. The evolution of the metamodels with respect to number of training examples is shown. The performance of different information transfer strategies is finally compared.

Typical results with the strategy GP

The metamodeling strategy is illustrated with the results from the GP-CK or GP-SK variants, using increasing numbers of training examples N R , always chosen among the 475 available for the random cell. However all 9649 examples from the unit cell are systematically used. When increasing N R , the dataset is growing in an inclusion sense : for instance between N R = 10, and N R = 20, the first ten points are kept. Data is used in the ordering given by the initial Sobol sampling (see 4.3.3). The different metamodels for different values of N R are constructed independently.

The evolution of the metamodels Ê * U/R (T, L) with respect to N R is shown in Fig. 4.6. The metamodels are only represented for T = 1, but the metamodels do make predictions in the 2D input space (T, L). The plotted confidence interval corresponds to the variance explained by the Gaussian process, excluding the noise. The metamodels are convergent with N R , since the predictions with increasing N R get closer to the predictions from the N R = 475 metamodel. The difference between cokriging and kriging is mostly found at low data. Cokriging allows to reproduce the trend given by the unit cell (decreasing then increasing with L, with a minimum near L = 0), which is missed by simple kriging for low values of N R . This overall trend is quickly found by cokriging, more quickly than the average level. For N R = 10, results for random cell are for instance lower than for unit cells, contrary to what happens for high values of N R . This underestimation at low data appears mostly due to chance: Other samples of N R = 10 points were investigated, and higher E R values were found.

The expectation of the full metamodel (with N R = 475) is represented in Fig. 4.7. Failure strain values are generally decreasing with T , and form a valley with L, with a minimum near L = 0 at given T . This is consistent with typical parametric models, and the coarse estimations found in Fig. 4.5. However a clearer and smoother evolution is found with the kriging metamodel than for the coarse minimization.

The metamodel should also be compared to some ground truth. A real ground truth is difficult to obtain, as FEA simulations can only give values for E R (T, L, d ) for specific loading conditions, and not the minimized values E * R (T, L). However The latter metamodel received no information from the unit cell. The GP-CK and GP-SK metamodels closely follow simulation results. They are actually slightly below "ground truth" which is not necessarily a problem. The minimum over 30 orientations is by definition a slight upper bound of the minimum over all orientations. On the other hand, the fully uncoupled strategy NT -SK is significantly lower than test data, and does not show the same trends, especially for L > 0. Therefore information transfer from unit cells can help reach satisfying agreement with reference computations.

Comparison of information transfer strategies

The above results with GP -CK showed that the unit cell can provide useful information. The different information transfer strategies are here more carefully compared. There are six possible variants: three choices for parameter transfer (NT, LS, GP) at the 5D level and the use of either cokriging (CK) or simple kriging (SK) at the 2D level.

To compare the strategies, two metrics comparing the metamodel E * R (T, L) to a reference E * R,re f (T, L) are defined. The reference is the metamodel using GP-CK on all 475 examples, which was shown in Fig. 4.8 to accurately reproduce direct results from simulations for T = 1. In order to compute the metrics, a regular 50 × 50 grid (T i , L i ) is introduced. The two Gaussian processes E * R and E * R,re f are evaluated on this grid of size N eval = 2500. The two metrics are defined as follows:

RMSE = 1 N eval ∑ i (µ i -µ re f i ) 2 1/2 (4.19)
KL = 1 N eval ∑ i KL (N (µ i , s 2 i ); N (µ re f i , (s re f i ) 2 ) (4.20)
where µ i and s i are the expectations and standard deviation of the random variable E * R (T i , L i ). The first metric only compares the expectation of the two Gaussian processes. The second metric is the average of the pointwise Kullback Leibler (KL) divergence [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], and allows to compare not only the expectations but the whole Gaussian distributions. All points are however considered independent and no correlation is taken into account. For Gaussian distributions, the KL divergence is given by:

KL (N (µ, s 2 ); N (µ , (s i ) 2 ) = 1 2 s 2 + (µ -µ ) 2 (s ) 2 + ln s 2 (s ) 2 -1 (4.21)
Moreover, in subsection 4.5.1, the overall trends were seen to be differently recovered than the average level. It is therefore worthwhile to consider normalizing the predictions of the metamodel. This normalization is performed by subtracting the mean of the predictions on the grid and dividing by their range. The pointwise variances can also be scaled in accordance. The pointwise expectation and standard deviation are thus given by:

µ scaled (T, L) = (µ(T, L)-< µ >)/(µ max -µ min ) s scaled (T, L) = s(T, L)/(µ max -µ min ) < µ > = 1 N eval ∑ i j µ(T i , L j ) µ max = max i j µ(T i , L j ) µ min = min i j µ(T i , L j ) (4.22)
With this normalization, the effects of the average level and of the range over the (T, L) space are neutralized: Two Gaussian process related by an affine transformation (in the output space) are considered equivalent. Only the overall trends and the shape of E * R (T, L) up to a homogeneous dilation are relevant. Errors will be thus presented both for these scaled results and the direct results (without the normalization).

In the following, as in subsection 4.5.1, the training is performed on a subset of the 475 available examples, considered in a given order. N R = 20 means, that the first 20 examples, according to this order are used for training. The set of considered inputs is increasing (not only in size, but according to inclusion). A strategy is efficient if its predictions quickly converge to those of the reference, that is if the distance measured by the metrics is low, even for a low data number.

Results for the initial Sobol sampling order are given in figure 4.9, comparing the six strategies for the two given metrics. Firstly, the results from the RMSE and KL metrics are very similar and yield the same global evolution. The expectation of the Gaussian process is sufficient to compare the strategies and the local variance does not provide significant additional information. With increasing N R , the errors are mostly decreasing, although dispersion can lead to local increases such as for N R = 40, if new examples change significantly the model.

The different strategies are not equally efficient. The strategy GP, both for simple kriging and cokriging, leads to smaller errors than the strategy NT without parameter transfer. However, the strategy LS is worse than the two other types of the parameter transfer. This is true both for direct results and the scaled results. This difference between the two strategies with parameter transfer indicates that the minimum E * R over all orientations is more accurately determined when the overall correlation patterns in 5D are taken from the unit cell, even if it means not interpolating between observations and explaining the difference by noise.

Moreover, using cokriging leads to increased efficiency, and especially for scaled results. A large difference between simple kriging and cokriging is noticeable at low data. The information from the unit cell allows to recover the overall trend, which can be hard to find at low data. However for direct results, the difference between simple kriging and cokriging is very small. Errors are then mostly due to average failure strain level, which is harder to determine, as illustrated in 4.6.

In order to strengthen these conclusions, the same investigation is performed for ten different shuffles of the random cell dataset. The order in which examples are used is therefore different, which affects especially models at low data, but all metamodels should still converge to the same one. Fig. 4.10 represents, for the GP and NT strategies, the minimal, maximal and the average errors over all ten shuffles. Results for the LS strategy, which were consistently worse than NT , are not represented. As above, the strategy GP is significantly better than NT especially with higher data. Cokriging leads to a slight reduction of errors compared to simple kriging. However there is a large dispersion of errors due to the different ordering of the examples.

The above results evidence the usefulness of the multifidelity approach using results on the unit cells to better approximate those on the random cell.

Adaptive sampling strategy

Results from 475 simulations may represent a low amount of data in 5D, obtaining this data is still computationally expensive. It is thus worth investigating whether the number of computations could be reduced further. Last section showed that the order in which examples are considered lead to a large dispersion of errors. Therefore a better ordering of data or more generally a more efficient choice of inputs for the simulations, using the results from the metamodel, could lead to lower errors with reduced numbers of simulations. This is the principle of Bayesian experimental design [START_REF] Chaloner | Bayesian Experimental Design: A Review[END_REF][START_REF] Frazier | Recent Advances in Optimization and Modeling of Contemporary Problems[END_REF]. Several adaptive experimental designs are presented and their performances are evaluated.

Definition of the adaptive sampling strategies

Two main strategies for an adaptive experimental design are here presented. They principally deal with the 5D kriging level. The posteriors of E * R (T, L) and E R (T, L, d ) have been computed on the account of some observations, and a new data point (T opt , L opt , d opt ) should be determined so that a new simulation with this loading condition is likely to provide as much information as possible.

The approach of "maximal variance" consists in finding the new point where uncertainty is maximal, that is:

(T opt , L opt , d opt ) = argmax Var E R (T, L, d ). (4.23)
This approach does not strictly require the values from E R (T, L, d ), but only the locations of the observations in the input space. This approach tends to find new points far from the observed data, and therefore to explore a greater portion of the input space. The approach of "Maximal Expected Improvement" tries to find locations for which the failure strain is likely to be lower then the currently predicted minimum, by maximizing the expected improvement compared to the current minimum:

(T opt , L opt , d opt ) = arg min E (E R * (T, L) -E R (T, L, Q ∼ )) + (4.24)
where x + = max(0, x). As E R (T, L, Q ∼ ) is a Gaussian variable of, say, mean µ, and standard deviation s:

E (E R * (T, L) -y) + = (y -µ)Φ y -µ s + sφ y -µ s (4.25)
where φ and Φ are respectively the probability density and cumulative distribution functions of the standard unit normal distribution. The maximal expected improvement approach refines the estimation of the currently found minima, and focuses on (or exploits) zones where a minimum is likely to be found. Other similar approaches could have used the probability of improvement or the lower bound of the confidence interval [START_REF] Frazier | Recent Advances in Optimization and Modeling of Contemporary Problems[END_REF]. Both criteria are applied on a number of evaluation points, which are compared. If a batch of N new points is desired, the N points, among the evaluated ones, with the highest variance (or expected improvement) should be chosen. In order to balance exploitation and exploration, a combination of the two strategies can be used. For a batch of N points, this mixed strategy leads to taking N/2 points with the highest variance and then N/2 points with the highest expected improvement. The computations can then be launched on the N new input points. The next metamodeling operation then uses the union of the data previously considered. This process can then be repeated until the desired number of data has been reached.

Comparison of strategies

The different methods for adaptive experimental design are now compared. As the information transfer method GP -CK was found to be the most efficient, only this method will be used in the following. For practicality, the possible new evaluation points are chosen among the 475 data already available for the random cell. All the results from the required simulations have thus already been produced, and no new simulation is needed at each phase. However, to find a new optimal evaluation point, the experimental design strategies consider only the location of the candidates in the input space and has no access to the associated failure strain. Therefore restraining the search to the available data, and only using failure strain results after they have been chosen by the strategy, can emulate a true experimental design with sequential batches of simulations.

Strategies are initialized by selected 10 points at random from the 475 points available. They are compared to a more naive strategy randomizing the ordering of examples (as in subsection 4.5.2), while keeping the first ten ones. The experimental design and the naive random strategy thus all start from the same initial data. The most efficient method will lead to metamodels converging faster, with lower errors at lower number of examples, to the reference, chosen as previously as the metamodel obtained by GP -SK on all 475 random cell examples.

In Fig. 4.11, the different methods are compared on three random choices of the ten initial points. The three adaptive experimental design methods perform similarly: they globally outperform the naive method for the scaled results, but are underperforming if direct results are considered. Cokriging and information transfer therefore appear to find effectively the overall trends, but they are less efficient at finding the average failure strain level, due to the complex operation of extracting a minimum. For scaled results, the method of maximum expected improvement appears the most efficient at finding the precise trends of the response surface. For instance, for N R ∼ 100, it has the same error as the random method with two or three times more examples. However, the difficulty of finding the absolute failure strain level may limit its usability.

Conclusion

The ductile failure of cells embedding a random distribution of voids was simulated for very general loading conditions, and showed a strong sensitivity to loading orientation indicating a significant anisotropy. A metamodeling strategy was proposed to construct an isotropic model from these computationnally expensive simulations and random cells. The main findings follow:

• The transfer of information from simulations on unit cells allows to build an accurate model by Gaussian Process Regression for random cells, despite a limited number of simulations for the latter, and the complexity of extracting a minimum over all loading orientations. The comparison is performed for three choices of first ten points (each choice of initial points corresponds to a row). The results are compared with the results from ten random orderings (with fixed 10 first points): the average is shown with a red line, while the minimum and maximum values delimit a red area.

• Several information transfer strategies between the two types of cells were considered.

Reusing the same Gaussian Process parameters from the unit cells and explaining the difference with random cells by epistemic uncertainty was found to be the most effective strategy.

• Cokriging always improves accuracy but is mostly effective for a very low number of examples and has a more limited effected at higher number.

• Strategies derived from Gaussian Process Regression were proposed to improve experimental design and obtain faster convergence of the metamodel with respect to number of examples. They allowed a better determination of the overall trend of the model, but were not found to improve the determination of the average level, limiting their efficiency.

The obtained response surface aims to be used directly as a non parametric uncoupled ductile failure model. Moreover, the methodology could be extended in several ways. For instance several different random cells could be used, so as to represent and quantify the dispersion of the strain at localization for different random microstructures. Although the metamodeling strategy was presented for simplicity on a single realization, an efficient metamodel should represent the average behavior of all possible realizations. Moreover the influence of the initial porosity could be assessed. Section 4.C presents some elements relative to this question, and especially shows that the effect of initial porosity can not be decoupled from the effects of other loading conditions. The effect of porosity could be more exhaustively studied by adding this porosity as a new dimension to explore by kriging. This requires a straightforward extension of the above methodology but also significantly more simulations on unit and random cells.

4.A Sparse Variational Gaussian Process

For the unit cell, thousands of examples are available. The standard Gaussian process methodology defined in subsection 4.4.1 leads to untractable computational times for maximum likelihood optimization, since an inversion (or at least a Cholesky factorisation) of the covariance matrix K ∼ xx is required (Eq. 4.8). All examples should however be used, so as to represent all the information available for the unit cell. Several adaptations to Gaussian Process Regression to big data are available [START_REF] Liu | When Gaussian Process Meets Big Data: A Review of Scalable GPs[END_REF]. The Stochastic Variational Gaussian Process (SVGP) method [START_REF] Hensman | Gaussian processes for big data[END_REF] was preferred, thanks to its formulation as a straightforward extension of the standard Gaussian Process, that can be used in its GPflow implementation [START_REF] Matthews | GPflow: A Gaussian process library using TensorFlow[END_REF]. The principles of this methodology are here briefly presented.

In a sparse Gaussian process [START_REF] Titsias | Variational Learning of Inducing Variables in Sparse Gaussian Processes[END_REF], the data from the examples is summarized by the value of the Gaussian process at a set of inducing points (z i ) i≤M ∈ X in the input space. Given the inducing points, the predictions of the Gaussian process for new points or existing examples (x i ) are given by:

f (x ) | f (z ) ∼ N (K ∼ xz K ∼ -1 zz f (z ); K ∼ + ς 2 I ∼M ) (4.26) where K ∼ = K ∼ xx -K ∼ xz K ∼ -1
zz K ∼ zx . The covariance matrix K ∼ xz is naturally defined by K xz,i j = K(x i , z j ) and similarly for the other matrices. The prediction is thus performed in exactly the same way as in Eq. 4.6. If the inducing points correspond exactly to the locations of the examples, the sparse Gaussian process is equivalent to the standard one. However, only the covariance matrix between the inducing points needs to be inverted (or decomposed) so with a lower number of inducing points, the computational cost of the prediction is decreased. The choice of the number of inducing points is determined by a trade-off between accuracy and computational cost (see section 4.B).

The SVGP is a sparse Gaussian process whose optimization is performed by a stochastic variational approach. The objective is to use a family of variational distributions to find a lower bound for the log-likelihood of the Gaussian process. [START_REF] Hensman | Gaussian processes for big data[END_REF] showed that for any Gaussian distribution N (m , S ∼ ) (in dimension M), the likelihood verifies:

ln P(y) ≥ N ∑ i=1 ln φ y i -k T i K -1 zz m ς - N ∑ i=1 Kii 2ς 2 + 1 2 tr(SΛ i ) -KL (N (m , S); N (0 , K zz )) (4.27)
where K is specialized to the examples (x i ), φ is still the probability distribution of the standard Gaussian distribution, k i is the i th column of K ∼ zx and

Λ i = 1 ς 2 K ∼ -1 zz k i k T i K ∼ zz (4.28)
The optimization of the log likelihood can thus be replaced by the optimization of the righthand side of Eq. 4.27 which acts as a lower bound. The parameters that can be optimized are those from the standard GP formulation (A, λ , ς ) and additional ones like the position of the inducing points (z i ), and the parameters m and S of the variational parameters. For the optimal solution, the values taken by the GP at the inducing points satisfy f (z ) ∼ N (m , S ∼ ) so they do not need to be directly optimized.

Only the two sums in the right hand side of Eq. (4.27) deal with the numerous training examples, and each term can be considered independently. This allows stochastic optimization and processing the examples by batches, reducing the computational cost of optimization. This stochastic optimization is performed with a hybrid approach [START_REF] Matthews | GPflow: A Gaussian process library using TensorFlow[END_REF], alternating steps of natural gradient optimization [START_REF] Hensman | Gaussian processes for big data[END_REF] only for the variational parameters and Adam [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF] stochastic optimization for all the remaining parameters (A, λ , ς , (z i )).

4.B Choice of hyperparameters for the 5D kriging on the unit cell data

This section discusses the choice of hyperparameters for the Gaussian Process Regression methodology which are used in the study. The main choices to make dealt with the optimization parameters for the SVGP on the unit cell, is number of inducing points, and the kernel used for Gaussian processes.

For the optimization of the parameters of the SVGP on unit cell data by Adam minimizer, batches of 200 examples are used and optimization is terminated when the lower bound on the log likelihood did not decrease by more than 10 % on the last 100 iterations. A learning rate of 1% is used for the Adam optimizer. This choice produces reasonable Gaussian processes; the resulting metamodels are not strongly dependent on the value of the above hyperparameters.

The influence of the kernel and the number of hyperparameters were more systematically analyzed. A three-fold cross-validation was performed on the 9649 unit cell examples. The examples are randomly partitioned in three folds, the model is trained on two folds and its predictions are evaluated (with a RMSE metric) on the remaining examples. This process is performed for each choice of held-out fold. The holding out of data is performed only in this section for hyperparameter selection, not in the rest of the study. The reference hyperparameters were M = 200 and a Matérn kernel of order 5/2.

Three choices of kernel were tried for the Gaussian Kernel: Matérn kernels of order 3/2 and 5/2 (Eq. 4.3) and a squared exponential [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]:

K Matern 3/2 (r) = 1 + √ 3r exp - √ 3r (4.29) K sq.exp. (r) = exp(r 2 /2) (4.30)
The squared exponential kernel favors smoothest realizations of Gaussian processes whereas Matérn-3/2 kernel leads to roughest realizations. Matérn kernels lead to significantly lower errors than the squared exponential (Fig. 4.12a). However, the difference between the two types of Matérn kernels is not significant. The Matérn-5/2 was favored in this study owing to the increased smoothness of its realizations.

Several choices of number of inducing points M are compared in Fig. 4.12b and 4.12c. In accordance with the discussion in section 4.A, the error decreases with the number of inducing points but the training time and computational cost increases with M. The value of M = 200 was chosen, because the error does not decrease significantly past this value, and the increased computational cost is therefore of little use.

4.C Influence of initial porosity

The methodology developed in this study allows to assess systematically the influence of the stress triaxiality, the Lode parameter and the loading orientation. However this does not exhaust all the parameters defining a simulation on a cell. For both unit and random cells, the initial porosity can influence localization. This was already illustrated in chapter 3 for an elastic-perfectly plastic material: a cell with initial porosity 1% is still liable to strain localization anisotropy, and strain localization values are higher than for an initial porosity of 6%. This appendix extends the analysis to a wider range of initial porosity values, in the context of the hardening material described in subsection 4.3.1. The objective is to discuss the integration of initial porosity within the surrogate model for strain localization.

The analysis focuses on unit cells. Seven unit cells with varying initial porosity f 0 ∈ {0.1%, 0.2%, 0.5%, 1%, 2%, 5%, 6%, 10%} are used in simulations. They are meshed with the same meshing parameters as in Fig. .12: Sensitivity analysis for the kernel and the number of inducing points of the SVGP. A 3-fold cross validation was performed. The RMSE score compares the predictions of the SVGP trained on two folds when compared with the reference values from the the third fold. For test score and fit time, the average and the standard deviation over the three folds are reported.

4.D. A SIMPLER PARAMETRIC MODEL

123 same methodology as previously, for loading conditions (T = 1, L ∈ {-1, -0.4, 0.4}), and the 30 loading orientations from Fig. 4.3 are tried. The results for f 0 = 6% therefore correspond to those already discussed. Fig. 4.13 then reports the minimum localization strain over all 30 loading orientations, depending on the initial porosity. For L = -1 and L = -0.4, localization strain mostly decreases with initial porosity which is consistent with the results from chapter 3. However, in the regime of low initial void volume fraction, failure strain stabilizes and appears to converge to a value independent of the initial porosity. For L = 0.4, failure strain is essentially independent on f 0 . Therefore the influence of porosity strongly depends on the loading conditions. The effect of initial porosity should be included in the failure strain model. However, the influence of initial porosity cannot be decoupled from the effect of the loading conditions (T, L). The model for f 0 = 6% can therefore not be simply extended with a tensorial form E * (T, L, f 0 ) = g( f 0 )E * (T, L, f 0 = 6%), where g is some function to be identified. The methodology of section 4.4 could admittedly be generalized by introducing f 0 as a new input dimension. This would require a new study, with a significant computational cost. For simplicity, the failure strain model identified with f 0 = 6% will be kept by remarking that failure strain is a decreasing function of f 0 and therefore the values obtained with f 0 = 6% can play the role of a lower bound for the true failure strain values.

4.D A simpler parametric model

In this chapter, a metamodel for E * (T, L) was obtained by Gaussian Process Regression. An advantage of the metamodeling approach is to provide a non-parametric model, which does not require hypothesizing an analytical form for the dependence on T and L. However the downside is that the metamodel is harder to use and implement in other contexts or software.

This short section is therefore dedicated to simplify the metamodel by fitting a prexisting parametric model. The resulting parametric model can not express the response surface with as much complexity as the metamodel, but may be easier to reuse. The metamodel should however be preferred whenever possible. Note that the parametric model could not be directly obtained from training examples, and the metamodeling approach remains necessary to perform the minimization over all loading orientations, and to take advantage from multifidelity.

The parametric model by [START_REF] Lou | Extension of a shear-controlled ductile fracture model considering the stress triaxiality and the Lode parameter[END_REF] was chosen, because it has a physical interpretation and naturally expresses failure strain. The model is given by:

E * (T, L) = C 3 2 max(0, 1 + 3T ) C 2 √ L 2 + 3 2 C 1 (4.31)
where C 1 , C 2 , C 3 are three parameters to determine. The metamodel obtained in subsection 4.5.1 for 475 points (depicted in Fig. 4.7), or rather the expectation of the Gaussian process, was evaluated on a grid (T , L). The coefficients C 1 , C 2 , C 3 were obtained by least square regression of log 10 (E * ). Due to Eq. (4.31), this corresponds to a simple linear regression. The optimal parameters were found to be:

C 1 = 3.9 C 2 = 1.9 C 3 = 0.47 (4.32)
The resulting E * (T, L) surface is shown in Fig. 4.14. The comparison with the original metamodel (Fig. 4.7) shows that the parametric model keeps the correct range of values but simplifies the dependence to T and to L. Notably the parametric metamodel cannot express asymmetry in L.

The original metamodel remains therefore preferrable.

Application to weld defects

This chapter investigates ductile failure in real TA6V weld joints. Contrary to the random microstructures considered in previous chapters, there is no separation of scales between defects and weld joint. The ductile behavior of welded specimens is experimentally investigated with in-situ tensile tests. The population of defects is identified and followed with synchrotron X-ray computed tomography. The analysis of macroscopic mechanical behaviors, of fracture surfaces and of the evolution of defects evidence ductile failure of the specimens. The experiments are also simulated by finite elements, taking into account the presence of defects. A good agreement is reached between experiment and simulation, both for the overall mechanical behavior and the local evolution of defects. The localization strain response surface obtained in earlier chapters for random microstructures is finally converted to a damage criterion, which is used to post-process and analyze the simulation results.

Résumé détaillé

Les chapitres précédents concernaient la simulation du comportement en rupture ductile de microstructures aléatoires statistiquement homogènes. Le présent chapitre s'intéresse désormais à la rupture de pièces soudées réelles. Ces pièces comportent un cordon de soudure qui, en plus d'être métallurgiquement (et donc mécaniquement) différent du métal de base, possède des défauts de soudure (c'est-à-dire des pores). Il n'est pas nécessairement possible d'identifier dans le cordon de soudure un volume élémentaire représentatif de la population de défauts, si bien qu'il est légitime d'étudier dans quelle mesure les résultats des chapitres précédents s'appliquent.

L'influence des défauts sur la rupture ductile est testée expérimentalement par des essais de traction in-situ sur des éprouvettes avec cordons de soudure. Les défauts de soudure sont identifiés et surtout suivis au cours de l'essai par tomographie aux rayons X (sur la ligne PSICHE du synchrotron SOLEIL). Les méthodes expérimentales et d'analyse des images sont présentées en section 5.3.

L'analyse des courbes expérimentales force-déplacement et des surfaces de rupture en section 5.4 met en évidence la rupture ductile. De plus, les images de tomographie permettent de suivre la déformation des défauts au cours de l'essai: les défauts, initialement quasi sphériques, s'allongent dans la direction de traction, tandis que leur volume augmente modérément. 125 Finalement, la section 5.5 se propose de reproduire les essais de traction avec des simulations par éléments finis. Les informations de tomographie sont prises en compte en maillant explicitement les défauts de soudure visibles sur les images au sein de l'éprouvette. La loi de comportement du matériau reprend le modèle identifié par Sarre (2018) pour le TA6V, en distinguant le comportement de la zone fondue et celui du métal de base. La comparaison simulation-expérience est effectuée aux niveaux macroscopique et local.

Au niveau macroscopique, un bon accord entre courbes force-déplacement est observé entre la simulation et l'expérience, même si la faible raideur de la machine de test complexifie l'analyse. Cependant les défauts semblent avoir un rôle limité sur le comportement mécanique global puisque des comportements mécaniques très similaires sont obtenus en négligeant les défauts. À un niveau plus local, l'évolution de la morphologie des pores est bien reproduite par la simulation. Une bonne représentation numérique de la surface de rupture nécessite la prise en compte explicite des défauts.

Finalement, la rupture ductile des éprouvettes soudées est reliée à celle des microstructures aléatoires des chapitres précédents. À partir de la surface de réponse de la déformation à localisation, un modèle découplé de rupture ductile peut être implémenté. Il peut alors servir à analyser par post-traitement les résultats de simulation d'éprouvettes. Au delà de l'illustration de l'utilisation d'un modèle développé à partir de simulations sur cellules, ceci permet de déterminer plus finement la surface de rupture simulée, qui ne coïncide cependant pas avec la surface de rupture expérimentale.

Introduction

In the preceding chapters, ductile failure was investigated in random microstructures. Although it took into account the interaction between defects, the material was described as a cell with periodic boundary conditions, i.e. a representative volume element. This implied that the population of defects was considered statistically homogeneous and that there was a separation of scales between the defect and cell level, and the macroscopic level of the mechanical parts. This chapter now deals with the ductile failure of a material with a population of defects, coming from industry-relevant weld joints. In such a weld joint, the population is no more statistically homogeneous, and can not be easily summarized by representative volume elements.

The welding process under analysis consists in assembling two parts made of titanium alloy TA6V by pulsed (Nd:YAG) laser welding with no filler metal added. The application of this process to TA6V was studied from a metallurgical point of view by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF]. Welding has several effects for the metal in the fusion zone, such as modifying the microstructure and therefore the mechanical behavior, or creating residual stresses. Moreover, defects can be created within the weld joint. We focus on weld pores, generated by gas bubbles trapped in the fusion zone during the solidification zone. The population of such defects was studied by [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF]. These defects are responsible for a weaker mechanical resistance, for instance in ductile failure, as they represent voids and strain concentrators. Investigating the change in ductile failure behavior due to the weld pores would therefore allow to better predict mechanical behavior and resistance to ductile failure in monotonic loading.

The weld pores can be detected with X-ray computed tomography (CT), allowing to obtain 3D images of the weld joint. Non destructive control can be routinely done with laboratory tomographs, whose low resolution allows detecting only the largest voids, which are a priori the most dangerous. In order to use as precisely as possible the easily available low resolution images, the precise influence of defects on ductile failure should be investigated. This can be done by monitoring the evolution of the population of the defects throughout a mechanical test. However this study requires higher resolution CT-images, such as that produced by synchrotron tomography.

For ductile fracture, CT is generally used to analyze the local phenomena involved in failure, for instance by [START_REF] Shen | Threedimensional quantitative in situ study of crack initiation and propagation in AA6061 aluminum alloy sheets via synchrotron laminography and finite-element simulations[END_REF], even obtaining grain-scale information [START_REF] Naragani | Void coalescence and ductile failure in IN718 investigated via high-energy synchrotron X-ray tomography and diffraction[END_REF]. Simulations of ductile fracture based on CT-images have also been performed. For instance O' Keeffe et al. (2015) or [START_REF] Doroszko | A new numerical modelling method for deformation behaviour of metallic porous materials using X-ray computed microtomography[END_REF] have simulated the propagation of a ductile failure crack. [START_REF] Buljac | Numerical validation framework for micromechanical simulations based on synchrotron 3D imaging[END_REF], [START_REF] Buljac | On the calibration of elastoplastic parameters at the microscale via X-ray microtomography and digital volume correlation for the simulation of ductile damage[END_REF], [START_REF] Buljac | Effect of void arrangement on ductile damage mechanisms in nodular graphite cast iron: In situ 3D measurements[END_REF][START_REF] Navas | An examination of local strain fields evolution in ductile cast iron through micromechanical simulations based on 3D imaging[END_REF] used digital volume correlations from laminography images (similar to laterally extended tomography images) to obtain local strain fields, which can be compared with finite element simulations of the observed area. However all these simulations requiring an explicit meshing of the defects are computationally expensive and focus on small regions. Marvi-Mashhadi et al. ( 2021) identified a statistical model for pore population from CT-images and performed finite elements simulations of ductile failure (in a dynamic loading context) for synthetic realizations of the population model. CT-images are therefore rarely used for predictions of the mechanical behavior of a given specimen. For instance, in the Sandia challenge [START_REF] Boyce | The second Sandia Fracture Challenge: predictions of ductile failure under quasi-static and moderate-rate dynamic loading[END_REF], CT-images were available for the specimens whose resistance was to be predicted. However such rich images were only used to determine porosity levels to be integrated in the ductile failure models. Similarly [START_REF] Madison | Corroborating tomographic defect metrics with mechanical response in an additively manufactured precipitation-hardened stainless steel[END_REF] only correlated image statistics to mechanical behaviors. The CT-images can be used in a predictive way with reasonable computational cost, by representing explicitly the largest voids, and representing the rest of the observed voids by a homogenized model, such as that developed in the preceding chapter.

In this chapter, in-situ tensile tests are performed for the first time on TA6V weld joints, allowing to follow the evolution of weld pores. Moreover the whole specimen is simulated by finite elements, and the effect of explicitly representing the largest defects is investigated. This allows comparing simulated and experimental results. We also illustrate the application of the ductile failure model from previous chapters to represent the influence of smaller voids.

Section 5.3 describes the material under analysis and the experimental methods for the in-situ tensile tests. In section 5.4, experimental results are presented and the ductile failure of specimens is evidenced from macroscopic stress-strain curves, fracture surface analysis, and monitoring of the defect population. Finally the experiments are simulated in section 5.5, allowing to compare numerical results and experimental ones. A ductile failure model is explicitly built from the results from previous chapters, and its application is illustrated.

Materials and Methods

This section describes the material under analysis and the experimental methods to investigate its behavior in ductile fracture. The titanium alloy and its known behavior as part of a weld CHAPTER 5. APPLICATION TO WELD DEFECTS joint are first recalled. The specimens used for tensile testing and the tensile setup installed on a tomography line are then described. Finally the methodology to analyze the resulting tomography images is defined.

Material and weld joint

The material used and its known behavior, both in the base metal and as part of the weld joint are first described. [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF] already performed a characterization of the alloy in this laser welding process, so the main information will be here summarized. The specimens described below are made from the batch investigated by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF] and [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF], so the characteristics of the alloy are identical.

The TA6V alloy (or Ti6Al4V or Grade 5) is a widely used titanium alloy. Its composition, as determined by the supplier, is given in Table 5.1. Thanks to its two main alloying elements, aluminum and vanadium, TA6V can have microstructures composed of two phases: a compact hexagonal phase α and a body-centered cubic phase β . The β phase has a higher yield strength than the α phase. The typical microstructure in the base metal is shown in Fig. 5.1. The equiaxial α grains of typical size 10 µm represent the major part of the microstructure. They are surrounded by a finer contour of β phase.

As studied by [START_REF] Helbert | The influence of stress triaxiality on the damage mechanisms in an equiaxed Ti-6AI-4V alloy[END_REF] and [START_REF] Helbert | Effects of microstructural parameters and back stress on damage mechanisms in α/β titanium alloys[END_REF], TA6V is prone to ductile failure. The relative influence of void growth and nucleation depends on the stress triaxiality. For T < 1, nucleation is known to dominate, whereas void growth prevails for higher stress triaxiality levels. Therefore, for a homogeneous tensile stress state (T = 1/3), nucleation should dominate, but void growth can play a larger role for notched specimens. Recall that the model developed in the previous chapters for random microstructures can represent void growth, but contains no ingredient for void nucleation.

Two TA6V parts are then assembled by pulsed laser welding. For reasons linked to the industrial process, a half lap joint is used (Fig. 5.2). The fusion zone does not cross the whole thickness, and a recess is left on one edge of the lap. This necessary recess acts as a crack in the assembly, and thus as a weak point.

The welding process is responsible for microstructural changes. Metallography images from the weld joint (Fig. 5.3) were obtained by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF] and show that an acicular microstructure, very different from the microstructure of the base metal, is created after solidification of the fusion zone. The grains are typically 100 µm in diameter, which is considerably larger than in the base metal. Between the fusion zone and the base metal, a heat-affected zone is present with two distinct microstructures. It stayed solid throughout the welding process, but was brought to a temperature high enough to modify the microstructure. The microstructural changes are responsible for a change in mechanical behavior, without even considering welding defects such as pores. [START_REF] Huez | Characterization of the Mechanical Behaviour of Both Fusion Zone and Base Metal of Electron Beam Welded TA6V Titanium Alloy[END_REF] showed for instance that the fusion zone has a higher yield strength than the base metal, although it has a slightly lower ductility. This situation, known as overmatch, entails that the plastic strain may concentrate in the base metal (depending on the loading conditions), and the base metal is more likely to fail. In such conditions, the resistance of the whole structure is more easily predicted since the behavior of the base metal can be easily investigated. However the defects' weakening effect remains to be determined.

Symmetry plane

Recess

Molten zone

Heat affected zone

Description of the specimens and experimental methods

The nominal geometry of the specimens for tensile testing is shown in diagrams 5.4. Except for the rectangular heads, the specimens are axially symmetric. The imaging by X-ray tomography requires thin specimens, so the radius in the most reduced section is less than 1 mm. The specimens used in the test were actually for the most part already imaged by [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF]. They were initially cylindrical but a central notch was added to favor failure in the central area, as the (2018) whole specimen was not observed by tomography but only a zone of length 2.6 mm. Moreover the reduced section was chosen according to the load capacity of the tensile testing setup. Finally the notch also helps in reaching higher stress triaxiality levels. The specimens are cut from welded TA6V plates. Depending on their orientation in the plate, specimens can contain a portion of weld joint along their main axis (thereafter referred to as longitudinal welds), or orthogonal to the main axis, in the central area (transverse welds). A schematic representation of possible weld joints is shown on the diagrams 5.4. Base metal specimens do not contain any weld joint portion. For some transversal welds, the geometry II with a second notch was used (Fig. 5.4b). This second notch was designed to promote fracture in the central weld zone, although the weld joint has a higher yield strength. For the present study, weld joints were machined to remove recesses, but the machining was eventually found insufficient, complicating the analysis especially for transverse welds. The description of all the specimens tested, with their geometry and the possible position of a weld joint, is presented in Table 5.2.

The tensile testing machine, nicknamed Bulky, was developed by [START_REF] Pelerin | Development of a Versatile Mechanical Testing Device for In Situ Synchrotron Tomography and Diffraction Experiments[END_REF] (see Fig. 5.5 for diagrams). The specimen is not anchored to the machine but its lower head is held by a hook connected to the stage driven by the stepper motor. The top head lays on a metallic part with an opening letting the central axisymmetric part of the specimen through. During the test, the motor moves the lower side of the specimen, whereas the top side is kept stationary. The tensile line is completed by the transparent tube of polymer (PMMA) surrounding the specimen, which is loaded in compression. The motor is controlled in displacement. A load cell is used to measure the force applied on the specimen.

The setup is installed on the PSICHE beam line of the SOLEIL synchrotron [START_REF] King | Tomography and imaging at the PSICHE beam line of the SOLEIL synchrotron[END_REF]. A high energy X-ray pink (polychromatic) beam is used so that high resolution images can be obtained with a low exposure time. The setup is placed on a rotating stage so that the specimen remains on the path of the X-ray beam to the detector (Fig. 5.6). The whole stage rotates with the tensile testing machine and the specimen. A CT-scan is performed by producing 2D X-ray images for various angles of rotation. A three-dimensional image of the specimen can then be reconstructed by Radon transform. The optical components of the beam line and the detector allow reconstructed images of size 2048 × 2048 × 2048 voxels. The resolution is 1.3µm/vx, so that the image represents a physical cube of side 2.6 mm. As the image is centered on the central notch of the specimen, the whole cross-section is visible, on a height of 2.6 mm. CT-scans are performed during the tensile experiment, without any interruption. As strain keeps increasing and is not kept at a fixed value, there can be no stress relaxation due to viscosity effects. The different 2D images for a scan are however performed at slightly different strain values. If the strain rate is low enough so that the difference in strain values is negligible over the time needed for a scan, the reconstruction errors are negligible. The number of scans performed for each specimen is indicated in Table 5.2.

Image analysis

The procedure for analyzing the reconstructed 3D CT-images is now outlined. The reconstruction from 2D X-ray images produces 8-bit grayscale images, where high values (in white) represent the zones with high absorption (and therefore the metal), whereas low values (in black) represent the voids (weld pores within the specimen, and the exterior). The image analysis begins by applying a Gaussian filter to reduce noise.

Reconstruction artifacts, and especially rings may be present in the reconstructions, due to deficient pixels in the X-ray detector. To remove these artifacts, a method analogous to that of [START_REF] Raven | Numerical removal of ring artifacts in microtomography[END_REF] was used. The axis of rotation is first located in the image, then higher frequency components are filtered out in the Fourier space. This process is carried out slice by slice, that is for each 2D section of the 3D images orthogonal to the tensile axis and the stage rotation axis. An example of filtered image can be found in figure 5.7. In order to isolate voids, generally darker than their background, the 3D image is segmented with the method by [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF]. A (black) top hat transform is applied: Closing operations are carried out to fill voids, then the initial image is subtracted (see [START_REF] Serra | Introduction to mathematical morphology[END_REF] for mathematical morphology terminology). Finally a closing operation with a small structuring cube of size 3 allows removing very small objects. A binarized image is thus obtained. By analyzing the connectivity of the binarized images, the different voids can be identified. For each void, physical quantities, such as the volume or the inertia matrix (see subsection 5.4.3), can then be directly computed from the image. A voxelized mesh of the defects can also be obtained with a Marching Cube algorithm [START_REF] Cline | System and method for the display of surface structures contained within the interior region of a solid body[END_REF].

The above analysis considers each CT-image as independent. However the in-situ tests produce several images for the same specimen, at several strain levels. There is therefore a correspondence between defects in two different images. Determining which defects in two different images are homologous, and determining the geometric transformation joining the two defect populations, is a problem known as image registration (see Fig. 5.8 for an illustration). This problem can be solved with an Iterative Closest Point (ICP) algorithm [START_REF] Besl | A method for registration of 3-D shapes[END_REF]. This algorithm was reimplemented and particularized to the ICP images, so its functioning is here described.

The objective of the ICP is to compare two point clouds A and B, and find the corresponding point in B of each point in A. As described in the pseudocode 1, the algorithm alternates two steps: the nearest neighbour of each (transformed) point in A is found in B, then the best affine function transforming the point A to the clouds of associated nearest neighbors in B is found by least squares optimization. This process is repeated until the iterations do not change anymore the association between points in A and points in B. The algorithm has then reached convergence, and the correspondance between a point in A and a point in B can not change anymore.

In our application to CT-images and to defect population monitoring, the first image is matched to each of the other successive images by applying the ICP with the cloud of defect

A B

Figure 5.8: Objective of image registration: finding the correspondance between defects in two images A and B, and the average deformation gradient between these two images.

Algorithm 1 Iterative Closest Point

Require:

Two point clouds A = (a i ), B. n ← 0 a (0) i ← a i for all i b (0) i ← the point of B, closest to a 0 i (nearest neighbor) repeat n ← n + 1 (F ∼ ICP ,t ICP ) ← the affine transformation f : a → F ∼ ICP .a + t ICP , minimizing the distance ∑ i f (a (0) i ) -b (n) i 2 b (n) i ← the point of B, closest to a (n) i (nearest neighbor) until a (n) i = a (n-1) i for all i return F ∼ ICP , t ICP , (b (n) i )
barycenters. For a defect in an image, the barycenter is naturally defined as:

x defect = 1 V defect ∑ v∈void x v with V defect = ∑ v∈void 1 (5.1)
where the sums are taken over all voxels constituting a defect. The ICP works with point clouds, rather than directly with CT-images. Digital volume correlation could be used on 3D images, as in [START_REF] Navas | An examination of local strain fields evolution in ductile cast iron through micromechanical simulations based on 3D imaging[END_REF], but there is not enough contrast in the bulk for accurate correlation (and the defects are too scarce to act as speckle). Moreover the ICP is computationnally simple, as there is only a reduced number of points (one barycenter per void) to consider. The point clouds A and B do not play a symmetric role as each defect in A is required to have a homologous defect in B, but not reciprocally. As defects can move out from the observed area during the test, the analysis is focused on the defects sufficiently far from the edges of the tomography images. Fig. 5.9 illustrates the application of the ICP algorithm to two defect clouds. The reference cloud A in blue does not coincide with the target point cloud B in yellow. However the optimal transformation of A in red, found by ICP, agrees well with the target point cloud. The affine transformation a → F ∼ ICP .a + t ICP can be physically interpreted. The difference between two images comes from two distinct phenomena. Firstly, due to experimental issues, the specimen is not necessarily in the same position for two different images, i.e. the position of the defects is expressed in a different frame for each image. For instance, the scan does not necessarily start when the rotating stage is at the same angular position (so there is a rotation) and the lower part of the sample moves, while the upper part is stationary. This corresponds to a translation and a rotation of the cloud of the barycenters. The second phenomenon has a mechanical interpretation and corresponds to the stretching of the specimen during the tensile test. Therefore, if it is assumed that the difference in imaging conditions only entails a rotation around the tensile axis (z) and a translation, a measure of strain along the tensile axis can be computed from the ICP results :

ε ICP = F ICP zz -1 (5.2)
However strain is not homogeneous in the specimen, and the macroscopic ε ICP measure does not represent local strain fluctuations. This also means that, for the ICP algorithm, the transformed point cloud does not exactly coincide with the target cloud.

Experimental results

The present section presents a direct analysis of experimental results. Ductile failure is evidenced from macroscopic strain-stress responses, and the evolution of defect populations is observed.

Stress-strain curves

The macroscopic mechanical behavior of specimens is first analyzed from experimental force displacement curves. The displacement corresponds to the displacement of the stage of the tensile testing machine, while the force is measured by the load cell. The force can be converted to a conventional stress measure by dividing the force by the area of the initial minimal cross section. This area is measured on the initial CT-image for each specimen and not directly from the nominal geometry (Fig. 5.4a) because of machining tolerances.

Stress displacement curves are represented on Fig. 5.10, and are grouped by similar behavior. All experiments begin by an initial linear elastic loading but the overall stiffness is low. A linear regression on the initial elastic part yields stiffness values of around 1.5 • 10 6 N/m whereas an elastic simulation of the specimen (as in section 5.5) shows that the specimens have a stiffness of around 2.3 • 10 7 N/m. The difference between these two values is due to the high compliance of the tensile testing machine. For the transverse welds (except A10), the failure is brutal, immediately after the end of the elastic part of the curve. For the others, a plateau that can correspond to significant plastic deformation is found. To understand the evolution of defects, this ductility plateau is a priori the most relevant. Maximal tensile stress values around 1100 MPa are compatible with the experiments by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF].

Some serrations in the force can be found on the plastic plateau (Fig. 5.11). They are systematically associated with the beginning and end of tomographic scans. However, they are not due to viscosity effects, as the machine keeps moving throughout the test, and there can be no stress relaxation. These serrations appear due to a softening when the X-ray beam is active, which can damage and heat materials. The effect of a small increase in temperature should be limited for the metallic sample, but not necessarily for the PMMA tube of the machine surrounding it. In Fig. 5.6, the PMMA tube has yellowed where it was traversed by the X-ray beam, a sign of aging. Therefore the heat softens the PMMA tube which transmits load throughout the tensile testing machine, and can be responsible for the serrations. Those are not desired, and are generally not seen when the tensile testing machine is used with polymer samples, allowing a lower energy X-ray beam.

Using the displacement of defects between CT-images, the ICP algorithm provides a measure of the strain in the observed area. Stress-strain curves for specimens with longitudinal welds are shown in Fig. 5.12. The analysis focused on the images of these specimens, which mechanical behavior showed significant ductility. Strains at failure of 15 to 20 % are obtained for the specimens. Due to different specimen geometries, and the presence of the weld joint, these values can not be directly compared to [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF] experiments. As on figure 5.10, a rather good repeatability between similar samples is found.

Fracture surface analysis

To ascertain the physical phenomenon responsible for failure, the fracture surface of each specimen was observed by scanning electron microscopy (SEM). These fracture surfaces were used without undergoing any specific preparation.

The failure surface of the specimen 07P3, consisting only of base metal (without a weld joint) is shown in Fig. 5.13a. A typical ductile failure pattern is found [START_REF] Benzerga | Anisotropic ductile fracture: Part I: experiments[END_REF], with a dimpled zone in the center, and smoother zones on the outside. For the sample A10, with a transverse weld, the exact same type of failure surface is obtained (Fig. 5.13b). Failure therefore took place in the base metal zone (and unfortunately not in the weld, as intended). This explains the similarity in mechanical behavior between the specimen A10 with 07P3 (Fig. 5.10a).

The failure surfaces for transversal weld joints are presented in Fig. 5.14. The difference of lighting between the images is irrelevant and due to the use of two different microscopes. All failure surfaces show the presence of a recess, that can be identified by a smooth uniform surface, where parallel lines corresponding to machining can even be found. This recess has probably caused the failure, as it acts as a weak point. The failure crack then propagated through the fusion zone, where clusters of weld pores can be found. The end of the failure surface crosses a zone with no void, i.e the base metal or at least the heat-affected zone. This propagation is consistent with Sarre's (2018)'s findings that the ductile failure favors propagating through the base metal rather than through the fusion zone, which displays higher yield strength. Fig. 5.15 shows more local aspects of ductile failure (for specimen A6). In the fusion zone, the welding pores are clearly visible. At their bottom, there is a geometrical crystallographic pattern, which may be associated with plastic slip lines. The pores are also surrounded by ductile failure dimples. The smoother terminal zone of defect-free base metal is also observed, and shows typical dimples of terminal ductile failure. Finally, the failure surfaces for longitudinal weld joints (Fig. 5.16) show two distinct zones. In the center of the surfaces, the fusion zone is visible with weld pores. Compared to the transversal weld joints (Fig. 5.14), pores are less numerous and larger ones are visible. These larger pores may have been the cause of the fracture initiation. The role of the pores in fracture will be analyzed more deeply with finite element simulations in section 5.5. On the edges of the failure surface, the zone devoid of pores corresponds to the base metal.

Analysis of defect populations

The defect populations and their evolution can be directly studied from the tomography images. The initial distribution of pores in a TA6V weld joint was studied by [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF], who also determined a statistical model to generate representative distribution of pores. This subsection will therefore focus on the evolution of the defects throughout the tensile test. The specimens EPF1, EPF3 and EPF4 with transverse weld joints were prioritized as they had shown significant ductility.

Populations of the defects from CT-images were extracted with the methodology described in subsection 5.3.3. As an example, the pore population for the specimen EPF4 is represented in figures 5.17 and 5.18, for the initial image (with a negligible load) and the last image (just before failure). As the ductile failure in the TA6V alloy is quite brutal and unstable, it is not possible to have images during the final failure process, but only just before. Fig. 5.17 is a 3D representation of the voxelized meshed defects. Fig. 5.18 is constructed by superposing all the images slices transverse to the tensile axis. More precisely, the intensity level of a pixel (x 0 , y 0 ) on the superposition image is given by:

I(x 0 , y 0 ) = min z∈[[0,2048]] I(x 0 , y 0 , z) (5.3)
where I(x, y, z) is the intensity level of the voxel (x, y, z) in the CT-image. As voids appear darker than metal, this operation allows combining all the voids in a single image, and therefore showing As already pointed out by [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF], weld pores mostly have a rather simple ellipsoidal geometry. During the experiment, they are seen to change shape, mostly by a lengthening in the tensile direction. The deformation appears more limited in the transverse directions (Fig. 5.18). The defects are not seen to coalesce and newly nucleated defects are not seen, although the resolution of CT-images is not sufficient to determine accurately the behavior of small defects or defects very close one to another. The change in the shape of defects is also exemplified in Fig. 5.19, which shows the evolution of the four largest defects in the specimen EPF4 for different macroscopic strain values. The pores lengthen along the tensile axis, and keep their mainly ellipsoidal shape throughout the test. The last defect of the image possesses however several lobes. The CT-images and the defect population they display can be studied from a more quantitative point of view. The evolution of the porosity can be easily extracted from the images. It is defined as the volume fraction of voids compared to the volume of the specimen in the observed area (excluding the exterior of the specimen). For voxelized images, the volume of the voids is directly given by the number of voxels occupied by this void.

The evolution of the porosity for three specimens is shown for the three specimens in Fig. 5.20. The initial void volume fraction is around 0.1% (on the whole observed area combining the fusion zone and the base metal). [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF] found an initial void volume fraction of 0.45%, computed on only the fusion zone. Both values are compatible, but the precise distinction between the fusion zone and the base metal is difficult. During the experiment, porosity increases by 20 to 40 % depending on the specimen. This increase is appreciable, but TA6V remains a material displaying low void growth. In steels for instance, void growth ratios of 2 to 6 were found by [START_REF] Benzerga | Anisotropic ductile fracture: Part I: experiments[END_REF].

The evolution of the shape of defects can also be quantified. [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF]

I ∼ defect = ∑ v∈ defect   (y -ȳ) 2 + (z -z) 2 -(x -x)(y -ȳ) -(x -x)(z -z) -(x -x)(y -ȳ) (x -x) 2 + (z -z) 2 -(y -ȳ)(z -z) -(x -x)(z -z) -(y -ȳ)(z -z) (x -x) 2 + (y -ȳ) 2   (5.4)
where the sum is taken over all voxels constituting the void, and ( x, ȳ, z) is the barycenter of the defect (see Eq. (5.1)). The "mass" of a voxel is not included in the above equation, as it is not physically relevant and all voxels are similar. The inertia matrix is real-symmetric so it possesses three real eigenvalues I 1 ≥ I 2 ≥ I 3 . From these principal moments of inertia, morphology indicators independent of the total volume of the inertia matrix can be obtained:

λ 1 = I 1 I 1 + I 2 + I 3 λ 2 = I 2 I 1 + I 2 + I 3 λ 3 = I 3 I 1 + I 2 + I 3 (5.5)
which verify the following properties (see [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF]):

1 2 ≥ λ 1 ≥ λ 2 ≥ λ 3 ≥ 0 (5.6) λ 1 + λ 2 + λ 3 = 1 (5.7) λ 2 ≥ 1 2 (1 -λ 1 ) (5.8)
These morphology indicators completely describe the shape of an ellipsoid (which is characterized by its equivalent ellipsoid). Note that the largest inertia indicator λ 1 is associated with the shortest principal axis of the ellipsoid. The evolution of the inertia indicators for the three specimens EPF1, EPF3, EPF4 is presented on Fig. 5.21. For the initial and last CT-images of each experiment, all defects with a sufficient size are represented in the (λ 1 , λ 2 ) space. λ 3 is not represented because its value can be deduced from λ 1 and λ 2 . All experiments show a similar evolution of the morphology of defects. Initially, the defect morphology is mostly characterized by λ 1 λ 2 1/3. This corresponds to almost spherical defects. Lacourt (2019) also found that the weld pores could be well approximated by spheres. However the pores change shape during the deformation, and their morphology indicators move mostly along the line λ 1 = λ 2 . This line describes prolate spheroids, obtained by elongating spheres along one axis. This elongation is coherent with the figures 5.17 and 5.19 which showed a lengthening of defects along the tensile axis. However, for both initial and final images, defects show a significant scattering. Defects show diversity in their geometry, although most of them can be described as spheres elongating along the tensile axis.

Simulation of experiments

The experimental results are now compared to finite element simulations. After describing the methodology of these simulations, the comparison is made both from a macroscopic (overall mechanical behavior) and a microscopic point of view (at the scale of voids). The role of voids in the ductile failure and their explicit representation is explored. Finally the ductile failure criterion associated to the cell simulations of earlier chapters is implemented and its use is illustrated.

Methodology

Firstly the methodology for the finite element simulation of experiments is explained. This implies precising the meshing strategy for the specimen containing defects, the material behavior laws, and the boundary conditions.

Mesh generation

In accordance with the previous section, only specimens with longitudinal welds were simulated, and especially the specimen EPF4. The general geometry was given by the diagrams in Fig. 5.4a. However due to machining tolerances, the real geometry was different from the nominal geometry. More precise values for the minimal cross-section of the specimen and the notch radius were thus found by examining the CT-images, and especially images slices of constant z (cross sections orthogonal to the tensile axis and the weld joint). The minimal radius of the cross-section was then extracted and the notch radius was found by studying the evolution of the cross radius with z and fitting a second order polynomial. The geometry could then be described as the union of solids of revolution.

As done by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF], the specimen is supposed to be made of two materials with different behaviors: the base metal and the fusion zone. As seen in Fig. 5.3, there is actually also a heat-affected zone (with possibly two sub-zones) with a possibly different mechanical behavior. However this last zone is not taken here into account. The distinction between the fusion zone and the base metal can be inferred from CT-images. To simplify, the fusion zone is supposed to have a prismatic shape, with the axis corresponding to the tensile axis (for longitudinal welds). Its section at constant z, is mainly delimited by two planes. Moreover the fusion zone should In order to represent explicitly the identified defects in the finite element simulations, the surface of the pores were meshed. A direct mesh was obtained from CT-images by the marching cube algorithm, but such a mesh is not adapted to simulations. It is voxelized, so a large number of elements is used and the surface of the pores is very rough. This complex geometry is not physical, but due to discretization effects in the CT-images. The mesh was therefore improved, by simplificating and smoothing it. With Vedo software (Musy et al., 2022), the mesh of the defect was first smoothed with a Laplacian algorithm, then decimated. The decimation operation strongly reduces the number of nodes, while attempting to keep the overall geometry. This process is illustrated on Fig. 5.23.

The representation of the overall geometry of the specimens through solids of revolution, the separation of the fusion zone and the defect surfaces can be combined with the hierarchical mesher Netgen [START_REF] Schöberl | NETGEN -An advancing front 2D/3D-mesh generator based on abstract rules[END_REF] in a mesh of the specimen with defects. Each of the above components participate in creating a surface mesh, then the volume of the specimen is constructed.

Even if the mesh of defects is simplified compared to the CT-image, it is still composed of approximately a hundred nodes. Moreover, defects have different sizes and small defects lead to small mesh sizes. Therefore, meshing explicitly all visible defects would lead to an intractably high number of finite elements. Only the largest defects (with highest volume) were therefore meshed. The volume of the defects acts here as an indicator of sensitivity to ductile fracture, but a defect criticality indicator could also be developed, as was done by [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF] in the Figure 5.23: Meshing of weld defects: initial voxelized mesh, smoothing and decimation context of fatigue. We chose therefore to mesh only the 10 and 30 largest defects. These numbers are admittedly arbitrary, but they allow representing respectively 72 % and 81 % of the porosity visible on CT-images (Fig. 5.24). Most of the porosity is therefore concentrated in a few large voids. Three meshes with 0, 10 and 30 defects were therefore created. As a damaging material behavior is used (see Eq. (5.16) below), the mechanical behavior is sensitive to mesh size. So as to compare the three meshes, they were built with identical sizes. However, mesh sizes far from and close to the defects are not the same. To ensure identical mesh sizes, for all meshes, the 30 defects were explicitly represented. For the mesh with no defects, the defects were then filled, i.e their interior was meshed too. The mesh with no defects does not contain any void, but it is locally refined where the defects should be. For the mesh with 10 defects, the same process was carried out for the defects numbered 11 to 30 (in the decreasing volume order). All three meshes are therefore identical outside the voids.

The central zone of the mesh with 30 defects is shown as an example in Fig. 5.25. The complete mesh can be seen in Fig. 5.26 below. Note that the separation between the fusion zone and the base metal is cleanly represented, and does not show any staircase effect, unlike Sarre (2018).

Material behavior

The base metal and the fusion zone have different mechanical behaviors. Their constitutive laws for simulation were identified by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF], from tensile and notched tension tests. Unlike the base metal, the intrinsic behavior of the fusion zone cannot be directly accessed by experiment. The behavior was thus identified on a synthetic equivalent of fusion zone: The base metal is heat treated in such conditions that the obtained microstructure is similar to that of a weld joint.

For both materials, the identified behavior is an isotropic elastoviscoplastic behavior with isotropic saturating exponential hardening. Ductile failure is modeled with a Gurson-Tvergaard-Needleman (GTN) criterion with nucleation. More precisely, the same corotational formulation and notations as for the random microstructures of earlier chapters (see subsection 2.3.2) are used (except that the stress should not be corrected by J = det F ∼ for porous materials) and the behavior is given by: ė ∼ = ė ∼e + ė ∼vp (5.9)

e ∼e = 1 + ν E s ∼ - ν E (tr s ∼ )1 ∼ (5.10) ė ∼vp = ṗvp ∂ ψ ∂ s ∼ (5.11) ṗvp = max 0, ψ(s ∼ ) K n
(5.12)

(5.13)

The decomposition in an elastic and a (visco)-plastic part is similar to the formulation of previous chapters. However, the viscoplastic strain rate is no more a plastic multiplier but given by Norton's (1929) law (Eq. (5.13)). This law allows representing the viscous behavior of TA6V, even at room temperature.

Contrary to the material behavior used for random microstructures, the von Mises yield criterion is not used but a GTN, adapted to ductile failure, is favored [START_REF] Gurson | Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I-Yield Criteria and Flow Rules for Porous Ductile Media[END_REF][START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]. The yield criterion ψ is thus given by:

ψ(s ∼ ) = σ * (s ∼ ) -R(p vp ) (5.14) 0 = s vm σ * (s ∼ ) 2 + 2q 1 f cosh q 2 tr(s ∼ ) 2σ * (s ∼ ) -1 -(q 1 f ) 2
(5.15)

s vm = 3 2 s ∼ dev : s ∼ dev
(5.16)

(5.17)

where f is the porosity, σ * is an implicitly defined effective stress measure and q 1 , q 2 are two parameters characterizing the yield surface shape. The isotropic hardening R(p) is given by

R(p vp ) = R 0 + R 1 (1 -e -b 1 p vp ) + R 2 (1 -e -b 2 p vp ) (5.18)
This double saturating exponential was already used for random cells in chapter 4. Ductile damage is modeled by the porosity f increasing. Two damaging mechanisms are considered: void growth and nucleation. Coalescence could also be added by accelerating damage for high porosity values with an effective porosity [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF], but this mechanism was not used by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF]. To represent void growth and nucleation, the porosity f is decomposed into a growth part and a nucleation part with prescribed evolution laws:

f = f g + f n (5.19) ḟg = (1 -f ) tr (ė ∼vp ) (5.20) ḟn = F n ς n √ 2π exp - (p vp -ε n ) 2 2ς 2 n ṗvp (5.21)
The nucleation rate is given by [START_REF] Chu | Void Nucleation Effects in Biaxially Stretched Sheets[END_REF] who suppose a Gaussian rate characterized by parameters F n (final nucleation porosity), ε n (strain of maximal nucleation) and ς n (width of the Gaussian). The values of the parameters for the fusion zone and the base metal are given in Table 5.3. In the nucleation behavior component, the parameter f 0 , not already introduced, corresponds to the initial porosity. For the fusion zone, Sarre (2018) identified f 0 = 4.5 • 10 -4 on heat-treated base metal, but used a value ten times higher in simulations, on the account of CT-images. As we intend to explicitly mesh weld pores, their porosity should not be included in f 0 : the value f 0 = 4.5 • 10 -4 is therefore used. The porosity of the unmeshed defects (especially in the meshes with 0 and 10 defects) could be accounted for in f 0 , but for simplicity, the value f 0 = 4.5 • 10 -4 is kept. The two constitutive laws are assigned to the elements of the relevant zones of the mesh in Fig. 5.26. As in the earlier chapters, quadratic elements with reduced integration are used. An implicit Euler method is used for the local integration of the constitutive laws.

Boundary conditions

Finally, boundary conditions are applied on both ends of the samples. A schematic representation of the specimen and its boundary conditions is shown in Fig. 5.26. A displacement is applied on the right side. Due to viscous effects, the displacement can not be chosen arbitrarily, and the same value as in the experiment is therefore imposed. In order to model the low stiffness of the tensile testing machine, the left side of the specimen is linked to a linear spring with a fixed end. The stiffness of this spring can be determined as follows.

As mentioned in the subsection 5.4.1, simulating the specimen with a defect-free mesh and a fully linearly elastic behavior showed that the specimen has a stiffness K sp = 2.3 • 10 7 N/m whereas an apparent stiffness of K app = 1.5 • 10 7 N/m was found from force-displacement curves. As the apparent stiffness can be modeled as the combination in series of the specimen stiffness and machine stiffness, the latter can be computed as:

K Bulky = K app K sp K sp -K app = 1.6 • 10 6 N/m (5.22)
As the low machine stiffness can lead to instabilities, some simulations were also performed by fixing the left side, instead of using a spring. In this case, the displacement rate was reduced so that the specimen undergoes the same strain rate as previously, at least during the elastic phase. The definition of the apparent stiffness allows determining elongation rate of the specimen from the total displacement rate (machine and specimen):

Ltot K app = Lsp K sp =⇒ Lsp = K app K sp Ltot (5.23) k u

Macroscopic comparison between simulation and experiment

The finite element simulations are first analyzed from the purely macroscopic point of view, by focusing on the overall stress-strain response.

The force through the specimen (or equivalently a measure of conventional stress) can be directly extracted from the simulations, by computing the nodal forces on the nodes on one side of the specimen. In the simulation, the machine displacement can be represented by the displacement of the nodes on one side of the specimen, the other side being fixed or held by the spring. However the strain value corresponding to the strain measure given by ICP is less straightforward, as it is not simply the average strain in the visible area. From the initial tomography image of a specimen, the position of the defects were extracted and were located in the corresponding mesh of the specimen. Only the barycenters of the defects were used; they did not necessarily correspond to nodes, but their displacement was computed by interpolation within elements. From the evolution of the point cloud of defect barycenters, an ICP process was performed. As the association between defects over different increments is already known, only the optimization of the affine transformation was actually necessary. This gave a strain value directly comparable to the experimental results. Fig. 5.27 shows the comparison between the experimental results and the simulation of the EPF4 test. The simulation was done with the defect-free mesh and the spring modeling the machine stiffness. A good agreement is found for the strain stress curve, although the final stress is slightly lower in the simulation. The final failure of the specimen could not be simulated, as the force drops suddenly (visible in the force-displacement curve). The simulation was therefore stopped when the Newton algorithm required too small timesteps for convergence. However the last simulated point and the last experimental images are found for similar strain values, so the simulation appears to predict reasonably well the onset of failure. However, simulated and experimental force-machine displacement curves do not agree as well. This might be due to the difficulty to simulate machine play in the initial phase and X-ray-induced force serrations (see subsection 5.4.1). 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 This sudden force drop for a given machine displacement is due to the low machine stiffness. When necking and damage start appearing in the specimen, the force in the tensile line decreases, so the spring modeling the machine starts shrinking even if the total displacement is constant. The specimen therefore continues to elongate, worsening its damage. This is the source of an instability leading to a sudden catastrophic failure of the specimen.

This instability is demonstrated by Fig. 5.28 which compares total displacement and strain within the specimen. The latter is seen to brutally accelerate. Only the viscous effects temper this acceleration, as a higher strain rate results in a higher viscous stress which helps reduce the force drop. The role of the machine stiffness in plastic instabilities was for instance studied by [START_REF] Petit | Impact of machine stiffness on "popin" crack propagation instabilities[END_REF]. The behavior of the meshes with 0, 10 and 30 defects is compared in Fig. 5.29, in simulations with and without machine stiffness. The strain stress curves coincide almost perfectly until the failure. Although simulations representing defects show an earlier failure, as expected, the difference in strain values is negligible. Therefore explicitly meshing defects does not appear necessary to represent the macroscopic mechanical behavior.

Local comparison: evolution of defects and failure surface

The macroscopic mechanical behavior was found to be very similar between the different meshes of the specimen EPF4, even for the one not representing explicitly defects. The comparison can also be made from a more local point of view, by studying the evolution of defects, which is known from the CT-images.

Firstly, the evolution of the morphology of the four largest defects in the simulation with 30

Table 5.4: Some antidivergences useful for computing volume (first row), barycenter (second row) and inertia matrix (last rows)

h H 1 (x, y, z)/3 x (x 2 , 0, 0)/2 y 2 + z 2 (x(y 2 + z 2 ), 0, 0) -xy -(x 2 y, y 2 x, 0)/4
defects is shown in Fig. 5.30. These four defects correspond to those labeled #1 to #4 in Fig. 5.19 and visible in Fig. 5.25. Their evolution is qualitatively very similar to the experimental one, shown in Fig. 5.19. The meshes in the simulation are simplified when compared to those from the CT-images, according to the procedure described in Fig. 5.23. More generally, the whole specimen mesh and its defects (Fig. 5.17b) can be compared to the last CT-image before failure: a good qualitative agreement is reached.

To study the evolution of defect morphology from a more quantitative point of view, the volume and the morphology indicators (λ 1 and λ 2 ) of defects can be extracted from simulations and compared to experimental results. Unlike CT-images where defects are represented by voxels, the surface mesh of each defect is known, and the displacement of each node is computed during the simulation. Volume information such as volume defect can then be extracted from the surface mesh by the divergence theorem. If h is any scalar field defined in the defect, and H is a vector field satisfying divH = h, the following relation holds:

V de f ect h dV = S de f ect H • dS ∑ tri∈S de f ect S tri H ( x tri ) • n tri (5.24)
where the surface integral is evaluated by summing the contribution of each triangular element of barycenter x tri , area S tri and outwards-pointing normal n tri . The Table 5.4 lists some fields h with a possible "antidivergence" field H . With these functions, the volume, the barycenter and the inertia of each defect can be computed from Eq. (5.1) and (5.4). The evolution of the volume and morphology indicators of the four largest defects is shown in Fig. 5.32. Results from the simulations with 10 or 30 defects are compared to those from experimental CT-images. For the defects #1, #2 and #4, the evolution of the volume and the shape is accurately described in the simulation. However, the defect #3 appears to stop evolving for strain values ε zz > 7%, contrary to what is observed in CT-images. As will be shown below, the difference is due to the location of defect #3. Moreover the difference between the simulation with 10 and 30 defects is not significant, so meshing more defects does not help representing more accurately the deformation of the larger voids.

Fracture surfaces and ductile failure criteria

From the simulations, the fracture surfaces can also be extracted. The Gurson porosity f is the natural indicator in the simulation to detect the position of the fracture surface. However another ductile damage criterion can be constructed from the simulations on random microstructures in As machine stiffness destabilizes simulations and hinders the accurate representation of ductile failure, simulations without springs are used from now on to study fracture surfaces. The presence of a spring may influence macroscopic quantities such as force through the specimen, but is not supposed to affect the location of failure. In Fig. 5.33, the fracture zone is defined by the elements where a local damage indicator is higher than some arbitrary threshold. For the moment the damage indicator is the total porosity f coming from Gurson formulation (Eq. (5.16)), which accounts for void growth and void nucleation. The results corresponds to strain values of ε zz ∼ 14%, corresponding to the beginning of the final softening of the specimen (see Fig. 5.29b). For the simulation without meshed defects, the fracture zone is simply orthogonal to the tensile axis, in the central notched zone of the specimen. For the specimen with 30 defects, the fracture zone is however more complex, as it tries and joins defects. It was not possible to set the threshold on f in such a way that the fracture surface be thin and completely cross the section of the specimen, so a value of f ≥ 0.2 was chosen. Both defects #1 and #2 (the largest of the specimen) appear to participate in this failure zone, so that this zone stays mostly orthogonal to the tensile axis. The true fracture surface is clearer in Fig. 5.33c, which focuses on the equivalent plastic strain (not a damage indicator): The failure surface goes through defect #1 and is therefore clearly oblique, rather than orthogonal to the tensile axis. Note that the defect #3 is excluded from the zone of maximum damage: after the onset of necking, it cannot deform further, in accordance with the results from Fig. 5.31. In the previous fracture surfaces, the value of f = 1, theoretically signifying the failure of the element, is never reached. This late failure is a known feature of the above Gurson criterion, therefore the coalescence is generally modeled with an effective porosity [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF].

Another damage indicator can be created by adapting the response surface determined in chapter 4 from simulations on random microstructures. [START_REF] Defaisse | Ductile fracture of an ultra-high strength steel under low to moderate stress triaxiality[END_REF] describe the procedure to convert a failure strain response surface in an uncoupled failure criterion. The In accordance with the previous chapter, the strain measure E is computed by projecting F ∼ -1 ∼ along the first stress principal direction (whereas [START_REF] Defaisse | Ductile fracture of an ultra-high strength steel under low to moderate stress triaxiality[END_REF] used for instance a von Mises equivalent strain). The element fails when it reaches the value D = 1. If the loading path is proportional, the denominator is constant so the value D = 1 is reached when E = E * loc . The elements are not eroded away when the value D = 1 is reached, so values D > 1 can be attained, although they correspond to broken elements.

The indicator D is an uncoupled failure criterion inasmuch as it is completely separate from plastic deformation and can thus be computed by post-processing after the finite element simulation. On the contrary the Gurson criterion in (5.16) couples plasticity and damage since the porosity f directly affects the mechanical behavior of the material, and its evolution needs to be computed as part of the simulation. The damage criterion of Eq. (5.25) was implemented in Z-set with the simplified expression of the localization response surface (see section 4.D).

The metamodel E * loc (T, L) created in previous chapters represents the localization taking place in random microstructures with initial porosity f 0 = 6%. Using this response surface as a model for the damage around weld pores implies several hypotheses, which do not strictly apply, but the application of this model can still provide insights on ductile failure mechanisms. Firstly, the weld pores are not randomly dispersed, and there is no clear separation of scales between the weld pores and the geometry of the specimen to define a representative volume element. Secondly, the initial porosity in the weld joint is not 6 %, but lower. However, the failure strain in random microstructures increases with a lower initial porosity, so failure strain values for f 0 = 6% can act as lower bounds for the real value. Finally the material used in chapter 4 for the matrix of random microstructures used the same elasticity and plasticity behavior as TA6V in a weld joint, but viscous effects and damage, either by nucleation or void growth, were not taken into account. The viscous effects should not be critical in the ductile failure behavior but the new failure criterion is unable to represent damage by nucleation, and partly able to represent damage by void growth: It represents the growth of initially present voids but not those newly nucleated. The new failure criterion therefore represents damage mechanisms complementary to the Gurson damage identified by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF], and the two failure criteria can be used in association. The main advantages of the failure criterion derived from random cell studies over the Gurson criterion are its ability to represent the interaction between voids, and its accounting for the effect of the Lode parameter.

The simulation of the specimen EPF4 with 30 voids and without machine stiffness was post-processed and the damage indicator D was computed. This indicator is only defined in the fusion zone, as the damage criterion was developed to represent unmeshed weld pores. In the simulation, the stress state (see Fig. 5.34) shows strong local variations near defects, allowing to test the effect of T and L on the ductile damage criterion.

The fracture zone corresponding to elements reaching D > 1 is shown in for the strain value ε zz ∼ 14%. As in Fig. 5.33, this corresponds to the failure of the specimen according to the macroscopic softening. At approximately the same time, the failure zone with D > 1 starts crossing completely through the thickness of the specimen, giving another sign of the fracture of the specimen. The fracture surface is finer and clearer than in Fig. 5.33 which used the Gurson porosity. Moreover this new fracture surface agrees with the one determined from equivalent plastic strain (Fig. 5.33c).

With the more clearly defined surface thanks to the new damage criterion, the simulated fracture surface can be compared to the experimental one (Fig. 5.36). However the simulation indicates a different fracture surface than the experimental results. Firstly, the macrophotography of the specimen shows that, due to ductile fracture, the experimental surface is wedge-shaped whereas the simulation finds a mostly flat surface. Moreover, both fracture surfaces contain a large weld pore, which is likely responsible for failure. In the simulation, this large pore corresponds to the defect #2. However in the experiment the defect #1 is more likely to have caused failure. By comparing the SEM micrography and the macrophotography (their orientation by the y-axis arrow), the central valley of the fracture surface descends from the large defect #b to the defect #c, whereas the simulation indicates that the surface goes up from the large defect #2 to the defect #a. This discrepancy makes it unlikely for the defect #2 to be responsible for failure. The experiment and the simulation both indicated that one of the largest pores initiated failure but the simulation did not replicate the correct fracture mechanism. The incorrect surface may be due to too large finite elements far from the defects, which do not allow the development of the correct mechanism during the simulation. Better results may be obtained with finer and more homogeneous meshes.

Conclusion

Tensile tests were performed on TA6V specimens containing a weld joint. In such weld joints, the fusion zone has a different mechanical behavior than the base metal, due to microstructural • The specimens failed by ductile fracture and a significant ductility was found for some specimens (with longitudinal welds). A fusion zone rich in weld pores is visible on the fracture surfaces.

• Weld defects were detected and followed throughout the tensile test. Defects are initially mostly spherical, and elongate along the tensile axis during the experiment. They also increase moderately in volume.

• The defects identified in CT-images can be meshed and embedded in finite element simulations of the experiment. Sarre's model was found satisfactory in reproducing the macroscopic mechanical behavior, even without explicitly considering the defects. Meshing the defects is more useful when studying local failure mechanisms.

• A good agreement was also reached between the observed and simulated evolutions of defect morphology, except that the experimental necking band may be wider than the simulated one.

• The localization strain response surface developed in earlier chapters was converted to a ductile failure criterion, and this criterion was implemented as a post-processing operation for finite element simulations. This criterion can be used in synergy with Sarre's model of coupled damage to analyze ductile failure. It aims to represent ductile behavior of the weld joint around weld pores and take into account the effect of the Lode parameter. Using this new damage criterion allowed determining more precisely the fracture surface in the simulations.

The above results should be confirmed on other experiments, as only one specimen was up to now simulated. Moreover, further simulations with finer and more homogeneous meshes should be performed to try and replicate the experimental fracture surface.

The experimental results can also be extended in several ways. Firstly, due to the limited ductility and the unforeseen presence of the recess, test on specimens with transverse weld joints were difficult to interpret New experiments taking into account the above issues may help in better comparing the mechanical behaviors of the base metal and the fusion zone. Moreover, the weld joints appear to have low influence on the mechanical behavior due to the low overall porosity. The ductile failure behavior of weld joints and the domain on which pores should be explicitly considered could be better studying by varying weld joint porosity, even if it means artificially introducing defects. Finally the ductile damage criterion constructed from response surfaces on random microstructures was here applied only to (notched) tensile tests. It should be tested on more general loading conditions and especially in shear conditions, so as to assess the influence of the Lode parameter.

Chapter 6

Conclusions and prospects

Conclusions

This work focused on the ductile behavior of materials with randomly dispersed voids, especially in the context of pores within weld joints. The original objective was to better understand the influence of preexisting voids on ductile fracture. Moreover it aimed to provide engineers with a method to efficiently use non destructive control information so as to predict the mechanical resistance under monotonic loading.

To this aim, the study focused on better understanding ductile fracture in the model micromechanical setting of a cell with a population of randomly dispersed voids, which is able to capture the key mechanisms of void interaction during failure. From finite element simulations on such cells, a ductile fracture criterion was also built by metamodelling. Ductile fracture was also studied experimentally on real weld joints, using tomography to observe the deformation of pores in situ. Simulating these experiments at the scale of the specimen allowed to make the connection with the micromechanical cell simulations.

Computational homogenization of random microstructures

A better understanding of the influence of the defect population on ductile behavior was reached by a micromechanical computational homogenization approach. Microstructures consisting of a cubic elastoplastic matrix with randomly distributed voids were created, and these microstructures were simulated by finite elements up to ductile fracture. Various loading conditions were investigated; they were characterized by proportional loading with constant stress triaxiality and Lode parameter and any loading orientation. This work thus provides the first ductile fracture simulations of random microstructures, subjected to any prescribed constant (T, L) loading, and especially accounts for the effect of loading orientation and cell anisotropy.

The behavior of unit cells (with a single void, commonly used in similar studies) and random cells was compared. For both types of cells, failure was found in extension mode, in shear mode, or a combination of these two modes. Strain localization bands are formed during failure, and complex patterns are obtained in random cells. Whereas the dependence of failure strain on stress triaxiality was similar for unit and random cells, their behavior with respect to the Lode parameter 167 was found to be quite different: Random cells are more likely to fail early for shear-dominated conditions. This difference is due to the random cells' greater freedom in expressing shear localization. Using unit cells therefore restricts this freedom and oversimplifies the localization mechanisms and possible void interactions.

Two different criteria allowed to detect the macroscopic failure of the cell, and were compared. A new coalescence criterion was constructed from a simple geometric interpretation of the form of the macroscopic gradient deformation rate. A localization criterion, derived from Rice's mathematical analysis of ellipticity loss, was reimplemented and used for the first time to study random microstructures subjected to arbitrary loading. This macroscopic localization criterion was found more suited to complex loading conditions and also allowed to determine quantitatively the localization mode: shear, extension or combined. It should thus be preferred in future studies.

The effect of loading orientation with respect to the frame naturally associated with the cubic cell was also carefully studied. A strong anisotropy was found for localization strain, whereas ductile behavior should only depend on stress invariants for an isotropic material. This anisotropy is directly related to the boundary conditions and the cell's symmetry. Moreover the unit cell displays higher anisotropy than random cells. In order to obtain an isotropic response, the minimal failure strain should be considered. Only considering diagonal loading then results in a strong overestimation of failure strain.

The statistical dispersion arising from realizations of random microstructures was also evaluated. This dispersion is significant, and remains high even with increasing defect population size.

Multifidelity metamodel for ductile failure

Simulations on random microstructures do not only give qualitative insight on ductile fracture, they can also form the basis of a ductile damage criterion. Such a criterion can be constructed from the response surface of localization strain depending on loading conditions. Damage can then be computed by post processing finite element simulation results, and integrating elementary damage increases throughout the loading path.

An accurate description of the response surface for random microstructures is then necessary. A Gaussian process regression is used to construct the metamodel from simulations at various loading conditions. However due to the high dimensionality of the problem (5D due to loading orientation), and the high computational cost of simulations on random microstructures, a multifidelity approach was developed to accelerate the construction of the response surface. This multifidelity approach takes advantage of the information brought by numerous simulations on the unit cell which acts as a coarse surrogate, easier to simulate but unable to represent properly void interaction. Moreover this multifidelity approach is able to negate the anisotropy effect of cell simulations and obtain an isotropic metamodel thanks to a two-stage strategy. Minimum failure strains are extracted from a 5D response surface, then the isotropic model is constructed by cokriging the resulting minimum failure strains for the unit and random cells.

Several information transfer strategies were investigated. Transferring parameters trained on numerous unit cell simulations was found more efficient than using a linear model of coregionalization (or cokriging). The predictions of the metamodel were compared to test random cell simulations, and a better agreement was found with the multifidelity approach than with a metamodel only using random cell training examples. Finally methods for adaptive design of experiment were investigated. Whereas general trends were more quickly recovered than with naive sampling, the adaptive design of experiment did not allow converging faster to the absolute strain localization levels.

Experimental investigation of ductile fracture in weld joints

Whereas the previous results focused on simulated artificial cells, the ductile behavior of real weld joints was also investigated with an experimental approach. In-situ tensile tests were performed on TA6V welded specimens. The defect population and its evolution were observed with synchrotron tomography. These experiments constitute the first in-situ observation of the pore deformation in TA6V weld joints.

The macroscopic mechanical behavior and the analysis of fracture surfaces clearly indicated ductile fracture. On successive tomography images, defects were seen to increase in volume and elongate along the tensile axis.

The experiments were also simulated by finite elements. The largest defects were explicitly considered in the specimen meshes. The mechanical behavior laws for the fusion zone and the base metal were given by [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF] previously identified model. These are the first micromechanically-informed simulations of weld joints allowing the prediction of specific specimens' failure.

The experimental macroscopic stress-strain response was well reproduced by the experiments. Explicitly considering defects had little influence on this macroscopic response. A good agreement between experiment and simulation was also reached for the more local evolution of defect morphology. The uncoupled ductile damage model derived from simulations on random cells was implemented and helped analyze the simulation of tensile experiments. It can represent different damage mechanisms than Sarre's (2018) coupled damage model. However these results should be verified on more specimens, and finer meshes could be used to better replicate the experimental fracture surface.

Innovative general-purpose numerical tools developed during this work

Over the course of the study, several numerical tools were developed and implemented, and might be of interest for future research.

For the computational homogenization of random microstructures, a procedure to create and mesh a population of randomly dispersed voids in a matrix was developed. Moreover, Ling et al.'s (2016)'s strategy for constant stress triaxiality and Lode parameter loading conditions was extended to allow arbitrary loading orientations. Considering sufficiently numerous orientations is necessary to account for the anisotropy induced by the shape of the simulation cell. Finally, Rice's criterion for macroscopic strain localization analysis was implemented in the Zset software. This criterion, which considers the macroscopic tangent operator, can be used to detect ellipticity loss at the cell-level, and also to stop simulations early on. It can be used for any cell simulation, and not necessarily for ductile fracture of voided cells.

A Gaussian Process Regression framework was used to obtain a metamodel of failure strain for random cell simulations. The kriging methodology and the multifidelity approach are general and may be used in many contexts. In this work, an original strategy was proposed to obtain a metamodel of a minimum quantity over some input dimensions (whereas kriging is generally used to interpolate directly between examples). Furthermore, strategies for adaptive designs of experience were implemented.

Finally a workflow was developed for the analysis of tensile tests observed with X-ray tomography. CT-images were filtered then segmented in order to identify voids. The correspondence between defect populations on consecutive images was recovered with an image registration algorithm. Moreover, defects were meshed and integrated in specimen meshes for finite element simulations. Methods to quantify the evolution of defect morphology in experiments and simulations were also implemented.

Prospects

Although this work provided some insights on the influence of voids in ductile fracture and especially in weld joints, the results may be extended in many ways.

Extension of computational homogenization to other situations

Other void populations The computational homogenization considered a population of identical spheres. A large diversity of population characteristics could be also studied. Several defect sizes and ellipsoidal voids (the orientation becomes important) could be used for example, as was done by [START_REF] Vincent | Effective flow surface of porous materials with two populations of voids under internal pressure: II. Full-field simulations[END_REF], [START_REF] Khdir | Computational homogenization of plastic porous media with two populations of voids[END_REF][START_REF] Khdir | A computational homogenization of random porous media: Effect of void shape and void content on the overall yield surface[END_REF] for yield surfaces. The behavior of clusters of pores is also of interest. Micromechanical analyses could also be performed on realistic defect populations, observed by tomography images.

Larger void populations Larger defect populations should be investigated in cell simulations, so as to limit dispersion, and ascertain whether a limit behavior with negligible dispersion can be obtained in the infinite defect population limit. However, simulations with large defect populations suffer from prohibitive computational costs. New computational methods should therefore be implemented. Finite element simulations could be parallelized, but periodic boundary conditions can be an obstacle. Using Fast Fourier Transform methods, as Hure (2021) did, may provide an alternative path for large computations, but the computation of the macroscopic tangent operator for strain localization analysis is not directly possible.

Other material behaviors

The macroscopic strain localization analysis can be used to investigate different matrix material constitutive laws. For instance, a yield criterion whose surface hows vertices, such as Tresca's (whereas the von Mises criterion has a smooth surface), may be more prone to localization [START_REF] Needleman | Analyses of plastic flow localization in metals[END_REF]. Moreover, the matrix could be sensitive to damage, but the accurate description of the behavior past the local (but not macroscopic) ellipticity loss may require regularization. Furthermore, crystal plasticity could be used (as was done by [START_REF] Guo | A data-driven approach to predicting the anisotropic mechanical behaviour of voided single crystals[END_REF]), especially when considering the large grain size in weld joints. However this would require a fine distinction between crystal plasticity and simulation anisotropy. Finally, viscoplastic behaviors should also be investigated but Rice's localization analysis is not directly applicable for such rate-dependent materials.

Ductile failure models based on micromechanical simulations

Truly isotropic cells Periodic boundary conditions and cubic cells entail anisotropy in failure results. This anisotropy should be considered, lest the cell's resistance be overestimated. It can be reduced but not eliminated with larger defect populations. The present work proposed to consider all loading orientations, but this is computationally expensive and hinders the development of models. Simulation cells or boundary conditions inducing lower anisotropy should be developed. Porous balls were tried in this work with limited results, but the search could continue on other cell types. These cells should be highly symmetric and as close as possible to the sphere to limit anisotropy, while still allowing identification of opposite faces. Although still of cubic symmetry, Kelvin's space-filling tetrakaidecahedron (or truncated octahedron) could for instance be tried [START_REF] Jang | On the compressive strength of open-cell metal foams with Kelvin and random cell structures[END_REF].

Enriched uncoupled damage model An uncoupled ductile damage model was obtained from the localization results of cell simulations. Only the localization strain values were used. A richer model could integrate, still with a Gaussian Process Regression, not only the failure strain but also the failure mode, so as to provide several failure mechanisms. Moreover, the statistical dispersion originating from multiple realizations of a random microstructure could also be modeled, so as to provide a probabilistic failure criterion.

Full failure model from random microstructures The proposed uncoupled failure model focuses only on the cell's final failure, and does not represent the full behavior of the cell before failure. Creating a ductile fracture model coupling plasticity and damage, so as to reproduce exactly the macroscopic behavior of cells, remains an open question. Existing yield criteria, such as Keralavarma et al.'s (2020), which also interprets failure as strain localization, could be compared to random microstructure simulations and fitted. However the ductile damage model could also be constructed in a more data-driven way, as [START_REF] Guo | A data-driven approach to predicting the anisotropic mechanical behaviour of voided single crystals[END_REF] did with neural networks for unit cell simulations with crystal plasticity.

Experimental study of ductile fracture

More diverse experimental loading conditions The experiments showed that [START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF] model was able to reproduce correctly the behavior of TA6V in ductile fracture. However this model was identified and tested on predominantly tensile loading conditions, and may not be able to represent behavior in shear loading conditions. Yet, for the simulations on random microstructures, the Lode parameter had a strong influence. The behavior of TA6V both as the base metal, and in the fusion zone could therefore be experimentally investigated for a larger diversity of loading conditions.

Using a more ductile material Thanks to the tensile tests with tomography imaging, the evolution of defect morphology was observed. However, the instant of failure could not be captured, and defects had a limited influence on the ductile fracture behavior. If the understanding of the influence of defects is prioritized, rather than the realistic description of TA6V weld joints, other experiments could be performed to better observe ductile fracture. For instance, larger pores with higher impact could be obtained by artificially introducing defects. Moreover, fracture could be more easily investigated with more ductile metals, such as cast iron [START_REF] Buljac | Numerical validation framework for micromechanical simulations based on synchrotron 3D imaging[END_REF].

Prediction of mechanical performance and justification of exemptions

Transition between homogenized model and explicit meshing In the proposed approach, the largest defects were explicitly meshed, whereas the smallest were only represented by the homogenized behavior. Further investigations should be performed to formalize the threshold between the two representations. For instance, explicit meshed defects in a structure could be enlarged until the macroscopic response differs significantly from the homogeneous material's response. This distinction is of practical interest as defects smaller than this threshold can be more easily ignored, whereas larger defects require specific attention from engineers.

Defect criticality indicator Due to the limited influence of the voids in the tensile tests, the choice of the voids to explicitly represent in simulations had little importance. The larger voids were thus meshed for simplicity, but meshing the more dangerous defects would be more sensible. Determining the criticality of a given defect for ductile damage is however not straightforward, for it can be influenced by its shape, or by its neighborhood of defects (especially within clusters). Criticality could for instance be quantified by a cell simulation of a defect and its neighborhood (as for the random microstructures). Performing such simulations for each defect would however be difficult in a production setting. A machine learning approach could therefore be developed to try and predict the criticality of a cluster of defects. Krokos et al.'s (2022) methodology, initially developed for maximum von Mises stress, could for instance be adapted.

Use of radiography In this work, defects were assumed to be identified and located in space thanks to X-ray tomography. However, in practice, non-destructive controls often resort to radiography only, which identifies the largest defects but only provides two-dimensional information on their location. The preceding methodology can be used provided that the three-dimensional position of defects is estimated. As a model of the defect population was constructed by [START_REF] Lacourt | Étude numérique de la nocivité des défauts[END_REF], possible locations of the defect can be estimated by inverse modelling: among all defect populations yielding the same radiography, the one agreeing best with the generative model is chosen.

1 ,

 1 e 2 , e 3 ) = (e x , e y , e z ) Canonical orthonormal frame a = a i e i Vector A ∼ = A i j e i ⊗ e j Second

Figure 1 . 1 :

 11 Figure 1.1: Metallographic observation of a TA6V weld joint by pulsed laser welding. The fusion zone and base metal display different microstructures. Adapted from Sarre (2018).

Figure 1 .

 1 Figure 1.2: Schematic representation of the main phases of ductile fracture: void nucleation, growth and coalescence. Taken from Tanguy (1999).

  Microstructure lattice with 27 voids on a cubic lattice.

Figure 2 . 1 :

 21 Figure 2.1: Meshes of some microstructures repeatedly used in the study.

  Cumulative plastic field after coalescence for T = 1, L = -1.

Figure 2 . 2 :

 22 Figure 2.2: Typical computation results for the R1 microstructure.

Figure 2 . 4 :

 24 Figure 2.4: Identification of three ductility zones, with respect to L (R1 microstructure, constant triaxiality T = 1).

Figure 2 . 5 :

 25 Figure 2.5: Strain at coalescence in the T -L space, for the R1 microstructure. Left: Coalescence surface interpolated by Gaussian process regression for multiple loading cases (simulations shown as red points). Right: Projection of the simulation results on the L -E c plane.

  ) with ||[P]|| the standard euclidean 2-norm of [P]. If [P] and [P ] are proportional, cos(θ PP ) = 1.This quantity is extracted from Z-set computations using tools developed by[START_REF] Lacourt | Hyper-reduced direct numerical simulation of voids in welded joints via image-based modeling[END_REF].

  Figure2.6: Link between E c evolution with respect to L and cumulative plastic strain field shortly after coalescence (at strain E 11 = 1.1E c , which depends on the simulation). Fixed triaxiality T = 1

Figure 2

 2 Figure 2.8: Dispersion of strain, porosity and stress at coalescence for different loading conditions, when considering multiple (N=20) random populations of 27 defects ( : comparison with the results for unit cell)

Figure 2

 2 Figure 2.10: Comparison of the coalescence onset, as determined by the δ -indicator and the energy-based criterion. All computations at triaxiality T = 1. The hatched zones correspond to simulations for which no minimum of Ẇe / Ẇp was observed, and therefore no coalescence was identified by the energy-based criterion.

  Figure 2.11: Comparison of the evolution of power ratio Ẇe / Ẇp (left) and of porosity f (right), for the unit and R1 cells in two loading cases: T = 1, L = -1 and T = 1, L = 0.25

Figure 2

 2 Figure 2.12: Cumulative plastic strain fields after coalescence for a microstructure with 125 voids

  500 MPa, K = 343.5 MPa, n = 0.58, b = 10 or b = 200. The different yield functions are shown in figure 2.14a.

  Figure 2.14: Comparison of different hardening behaviors. (a): yield function for each hardening type. (b): Comparison of the E c -L curves for each hardening type on the microstructure R1. All computations at fixed triaxiality T = 1.

Figure 2

 2 Figure 2.15: Detection of failure through simple extension criterion (stabilization of the transverse strain) or vanishing of δ function. Microstructure R1 under the loading condition T = 1, L = -1

  Figure 2.16: Effect of varying threshold conditions for the failure indicator on the E c -L curve. All simulations on the R1 microstructure, at T = 1.

Figure 2 .

 2 Figure 2.17: Cauchy stress during the simulation for several meshes of the R1 random microstructure with different meshing parameters. Loading condition: T = 1, L = -1

  Figure 2.19: Influence of boundary conditions on the response of the cell. All computations at triaxiality T = 1.

Figure 2

 2 Figure2.20: Evolution of the strain at coalescence with respect to T for the cubic and 2D axisymmetric unit cells for porosity values f 0 = 6% and 1% (constant Lode parameter L = -1).

  Figure 2.21: Determination of the failure band and the surface porosity for a realization of a 27 void microstructure

Figure 2

 2 Figure 2.22: Surface porosity for planes orthogonal to the main loading axis, for three other realizations

  Figure 3.1: Meshes of the cells repeatedly used in the study. All cells have an initial void volume fraction f 0 = 6%

∼

  dev the deviatoric part of the rotated Cauchy stress tensor s ∼ , s vm the equivalent von Mises stress and ṗ = 2 3 ė ∼ p : ė ∼ p playing the role of the plastic multiplier. The Young modulus, the Poisson ratio and the yield strength are respectively chosen as E = 200 GPa, ν = 0.3 and R 0 = 500 MPa.

Figure 3

 3 Figure 3.2: Evolution of transverse strain, stress component σ11 and Rice criterion for a simulation of a unit cell at T = 1, L = -1. The stabilization of the transverse strain and the localization onset are almost simultaneous. This instant also corresponds to an inflection point of the stress-strain curve.

  Figure 3.3: Deformation at localization for unit and random cells (diagonal loading) and slices of the (T, L) space. The triangles indicate a slope of -3/2 in the log-lin plot: E loc ∝ exp(-3T /2)

1 (

 1 Figure 3.4: Characteristics of the localization band for the unit cell, at (T = 1, L = -1)

Figure 3 . 9 :

 39 Figure3.9: Link between the deformation at localization and the localization mode (T = 1). Same as Fig.3.3c, but with localization modes highlighted. The extension mode (resp. shear mode) is shown in blue (resp. green) and corresponds to g• n ≥ √ 2/2 (resp. g • n < √ 2/2 ) at localization

Figure 3

 3 Figure 3.11: Dependence of the deformation at localization for different rotations around the axes (three cells, T = 1, L = -1)

1 Figure 3

 13 Figure 3.12: Deformation at localization for 150 different loading orientations, represented in the disorientation space

Figure 3

 3 Figure3.14: ṗ/ Ė11 just after localization, for the orientations maximizing (left) and minimizing (right) E loc , for the two random cells and two loading conditions (rows). To better visualize the failure mode, the deformed mesh is represented, and a factor is applied to the displacements.

Figure 3

 3 Figure3.17: Anisotropy of the maximum principal stress, for the unit cell and two random ones. Computations at T = 1, varying L, for 30 different loading orientations.

  Figure 3.19: Validation of the computation of macroscopic tangent operator on homogeneous volume elements

Figure 3

 3 Figure 3.20: Validation of the computation of macroscopic tangent operator on an elastic unit cell

Figure 3

 3 Figure 3.21: Verification of the -E /4 theoretical localization result on a homogeneous cube: for a von Mises material with softening R(p), the localization takes place when dR/dp = -E /4

Figure 3

 3 Figure 3.22: Sensitivity of the determination of localization by Rice criterion to the mesh size and to the time step size

Figure 3 Figure 3

 33 Figure 3.23: Deformation at localization for the homogeneous cube of matrix and comparison with results from the unit and random cells (diagonal loading, T = 1, varying L)

  Figure 3.25: Localization strain for 30 loading orientations on the hollow sphere with periodic boundary conditions. The results are displayed only if a localization was found by Rice criterion.

Figure 3 .Figure 3

 33 Figure 3.26: Hollow ball with periodic boundary conditions. Left: mesh of an hemisphere. Right: ṗ/ Ė11 , at E 11 = 1 (last computed time step) for (T = 1, L = -1) and diagonal loading orientation

Figure 4 . 1 :

 41 Figure 4.1: Meshes of the unit cell and the random microstructures used in the study. The void volume fraction is 6%

4. 3 . 2

 32 Figure 4.2: Strain at localization for 30 different orientations, at T = 1 and L = -1. Representation in the cubic-orthotropic misorientation space.

  Figure 4.3: Dispersion of strain at localization for the unit and the random cells. Fixed triaxiality T = 1, varying L and 30 different loading orientations.

  Figure 4.4: Training data for the unit and random cells. Each point is projected in the T, L space, ignoring the three components of the loading orientation.

  Figure 4.5: Estimation of the minimal failure strain for the random and unit cells

  Fig. 4.3 depicting results from simulations at T = 1 for different values of L, and 30 different loading orientations, can give an estimation of E * R (T, L) in the T = 1 space. It is emphasized that the data from Fig. 4.3 were not used to train L

Figure 4 . 6 :Figure 4

 464 Figure 4.6: Application of strategy 1-SK and 1-CK with different numbers of examples on the random cell. Gaussian process predictions in the T = 1 plane.

Figure 4

 4 Figure 4.8: Comparison of the predictions from the metamodel to results from 300 experiments

Figure 4

 4 Figure 4.9: RMSE and KL errors for varying number of training examples for the different information strategies. All the strategies use the same dataset ordering

Figure 4 .

 4 Figure 4.10: Dispersion of RMSE errors due to dataset ordering. Ten random orderings were used; minimum, maximum and average errors are represented for each strategy.

Figure 4

 4 Figure 4.11: Comparison of the adaptive experimental design strategies. Evolution of the error with the number of examples, from direct (left column) and scaled results (right column).The comparison is performed for three choices of first ten points (each choice of initial points corresponds to a row). The results are compared with the results from ten random orderings (with fixed 10 first points): the average is shown with a red line, while the minimum and maximum values delimit a red area.

  Figure4.12: Sensitivity analysis for the kernel and the number of inducing points of the SVGP. A 3-fold cross validation was performed. The RMSE score compares the predictions of the SVGP trained on two folds when compared with the reference values from the the third fold. For test score and fit time, the average and the standard deviation over the three folds are reported.

Figure 4 .

 4 Figure 4.13: Influence of initial porosity on failure strain. Each datum corresponds to the minimum failure strain over the same 30 loading orientations as in Fig. 4.2.

Figure 4

 4 Figure 4.14: Response surface obtained by fitting Lou and Huh's (2013) model

Figure 5 . 1 :

 51 Figure 5.1: SEM micrography of base metal, taken from Sarre (2018). The microstructure is made of equiaxial α grains surrounded by a fine contour made of β phase (in light grey).

Figure 5

 5 Figure 5.2: Schematic representation of the weld joint, adapted from Sarre (2018)

  Figure 5.3: Optical micrographies of the microstructure in the weld joint. Taken from Sarre(2018)

  Figure 5.4: Nominal geometries of the specimens for tensile testing, and schematic representation of the weld joint orientation

  Figure 5.5: Diagram of the tensile testing machine. Taken from Pelerin et al., 2019.

  Figure 5.6: Experimental setup on the PSICHE line

Figure 5

 5 Figure 5.7: Example of image analysis: ring artifacts are removed, and defects are segmented.

Figure 5 . 9 :

 59 Figure 5.9: Illustration of the ICP on the point cloud of defect barycenters for two different images: the transformed point cloud A' is optimally close to the target point cloud B

  Figure 5.10: Experimental displacement-stress curves for the tested specimens

Figure 5

 5 Figure5.12: Stress strain curves for the specimens EPF1, EPF3 and EPF4. The stress is computed from the displacement of defects visible in the tomography images.

Figure 5

 5 Figure 5.13: Fracture surfaces showing ductile failure in the base metal (SEM)

Figure 5 Figure 5

 55 Figure 5.14: Fracture surfaces for specimens with transverse weld joints. The failure starts from the recess, continues in the molten zone with welding pores, and finishes in the defect-free base metal.

  Figure 5.18: Evolution of the population of defects from specimen EPF4. All transversal slices of the X-ray tomography are superposed, and the minimum of all intensity levels is taken, showing the whole population of defects.

(Figure 5

 5 Figure 5.19: Evolution of the four largest defects in samples EPF4. Each row corresponds to one defect. The different columns correspond to strain levels of approximately 0% (initial image), 1.5%, 6.5%, 11.5% and 14.4% (shortly before failure).

Figure 5 .Figure 5

 55 Figure 5.21: Comparison of inertia-based morphology indicators λ 1 and λ 2 between the first and last tomography images, for three specimens with longitudinal weld

Figure 5

 5 Figure5.24: Evolution of the cumulative relative porosity (as a fraction of the total porosity) with the number of defects, for the EPF4 specimen. Defects are ordered by decreasing volume.

Figure 5 .

 5 Figure 5.25: Central zone of the meshed sample, with 30 meshed defects. A wireframe visualization is used to show the voids, which are hidden inside the mesh of the sample. The four largest defects (shown in Fig. 5.19) are labeled.

Figure 5 .

 5 Figure 5.26: Mesh of the specimen with its boundary conditions. The fusion zone (in red) and the base metal (in white) are shown in different colors.

  Figure5.27: Comparison of experimental and simulation results for the EPF4 specimen. The void-free mesh with a spring is used for the simulation.

  Figure5.28: Acceleration of the strain increase in the specimen, due to the machine stiffnessinduced instability

  Figure 5.29: Comparison of the response of three meshes of the EPF4 specimen: with no defects, with the 10 or 30 largest defects explicitly meshed. Responses are almost indistinguishable.

  Figure5.30: Evolution of the meshes of the four largest defects in the simulation of specimen EPF4, with 30 meshed defects and machine stiffness. Each row corresponds to one defect. The different columns correspond to strain levels of approximately 0% (initial image), 1.5%, 6.5%, 11.5% and 14.4% (shortly before failure).

Figure 5 .

 5 Figure 5.31: Mesh of the specimen with defects at the end of the simulation. A good agreement in the shape of defects is found when compared to Fig. 5.18b.

Figure 5

 5 Figure 5.32: Comparison of the evolution of the volume (first column) and morphology indicators (second and third columns) in the simulation and in the EPF4 experiment Each row corresponds to one of the four largest defects. All results are normalized by the initial value.

  Figure 5.33: Fracture zones in the simulations. The elements with a total porosity f ≥ 0.2 or equivalent plastic strain p ≥ 0.3 are highlighted. The outline of the mesh and the defects is shown as a wireframe

  Figure 5.34: Local variations of the stress state in the specimen during ductile failure. Representation in a plane cutting vertically the specimen in half. For the Lode parameter, due to implementation choices, the opposite of our convention is used in the definition of L: uniaxial tension corresponds here to L = 1.

Figure 5 .

 5 Figure 5.35: Fracture zones in the simulations according to the damage criterion. The elements reaching a damage D = 1 are highlighted

Figure 5 .

 5 Figure 5.36: Comparison of the experimental and simulated fracture surfaces

Table 4 . 1 :

 41 Coefficients of the material model

	Young modulus Poisson Ratio Yield Strength		Hardening
	E	ν	R 0	R 1	b 1	R 2	b 2
	115 GPa	0.3	753 MPa	222 MPa 290 285 MPa 10

Table 5 .

 5 1: Chemical composition of base metal TA6V (from supplier's statement, cited by Sarre (2018))

	Element	Ti	Al	V	O	H	N	C	Fe
	% (wt) Base 6.42 4.02 0.19 0.003 0.0055 0.06 0.15

Table 5

 5

		.2: Description of the specimens
	Specimen Geometry	Welding	Number of Images
	07_P3	I	Base metal	30
	A3	I	Transverse	11
	A6	II	Transverse	13
	A7	II	Transverse	31
	A8	I	Transverse	14
	A9	I	Transverse	13
	A10	I	Transverse	38
	EPF1	I	Longitudinal	27
	EPF2	I	Longitudinal	67
	EPF3	I	Longitudinal	18
	EFP4	I	Longitudinal	17

Table 5 .

 5 3: Parameters defining the constitutive laws of the base metal and the fusion zone. Taken from[START_REF] Sarre | Influence du soudage laser Nd:YAG sur les propriétés métallurgiques et mécaniques de l'alliage de titane TA6V[END_REF]

	Behavior	Parameter	Unit	Base metal Fusion zone
	Elasticity	E	GPa	115	115
		ν		0.3	0.3
	Viscosity	K	MPa /s 1/n	310	310
		n		15	12
		R 0	MPa	700	753
		R 1	MPa	285	222
	Plasticity	b 1		10	290
		R 2	MPa	0	285
		b 2		/	10
		f 0		10 -5	4.5 • 10 -4
	Nucleation	F n		0.35	0.20
		ς n		0.10	0.038
		ε n		0.275	0.107
	Void growth	q 1		1.1	1
		q 2		1.08	1.13

Remerciements

RÉSUMÉ

Lors d'un assemblage par soudage, des bulles de gaz peuvent être piégées au sein du cordon de soudure. Ces défauts réduisent la résistance des pièces à la rupture ductile, mais peuvent être détectés par des moyens de contrôle non destructifs comme la tomographie X. Cependant l'utilisation directe des défauts dans une simulation mènerait à des temps de calcul prohibitifs. L'objectif de l'étude est d'analyser l'influence des pores sur la rupture ductile de pièces soudées, et de proposer une méthode permettant de prendre en compte de façon simplifiée les images de tomographies afin de prédire efficacement et précisément la résistance en rupture ductile. L'influence sur la rupture ductile est étudiée par des simulations micromécaniques sur des cellules élastoplastiques contenant une distribution aléatoire de pores, permettant de prendre en compte l'interaction entre défauts. La détection de la rupture à l'échelle de la cellule par coalescence ou localisation, l'influence de l'orientation du chargement et la dispersion pour différentes populations sont en particulier étudiées. A partir des simulations sur microstructures aléatoires, la dépendance de la déformation à localisation aux conditions de chargement peut être représentée par un métamodèle obtenu par krigeage. Une approche multifidélité est utilisée pour combiner les résultats de cellules unitaires et aléatoires. Enfin, la rupture ductile est étudiée expérimentalement par des essais observés par tomographie synchrotron, permettant de suivre l'évolution des défauts. Les résultats expérimentaux sont comparés à la simulation, à la fois au niveau du comportement macroscopique, et à celui de la déformation locale des défauts.

MOTS CLÉS

Rupture ductile, Homogénéisation numérique, Localisation de la déformation, Métamodèle, Tomographie X ABSTRACT During welding, gas bubbles can be trapped within the weld joint. These defects reduce the resistance to ductile fracture, but can be detected by non-destructive controls such as X-ray tomography. However, the direct use of defects in a simulation would lead to prohibitive computation times. The objective of the study is to analyze the influence of pores on the ductile fracture of welded parts, and to propose a method to take into account the tomography images in a simplified way so as to predict efficiently and precisely the resistance in ductile fracture. The influence on ductile fracture is studied by micromechanical simulations on elastoplastic cells containing a random distribution of pores, allowing to represent defect interaction. The study focuses on the detection of fracture at the cell scale either by coalescence or localization, the influence of the loading orientation and the dispersion for different void populations. From simulations on random microstructures, the dependence of the localization strain on the loading conditions can be represented by a metamodel obtained by kriging. A multi-fidelity approach is used to combine results from unit and random cells. Finally, ductile fracture is studied experimentally by means of tests observed by synchrotron tomography, allowing to follow the evolution of defects. The experimental results are compared to the simulation, both at the level of macroscopic behavior, and at the local level of the deformation of defects.

KEYWORDS

Ductile fracture, Computational Homogenization, Strain localization, Metamodelling, X-ray tomography