Remerciement

En premier lieu, je tiens à remercier les rapporteurs et membres de mon jury. Je souhaite adresser ma gratitude à Eric Lorentz, pour son rapport méticuleusement détaillé, qui m'ont permis de continuer ma formation jusqu'à mon dernier jour de doctorant.

Je souhaite également remercier Yann Monnerie, d'avoir pris et dument assuré la tâche difficile de rapporteur et président de jury. Je souhaite te remercier pour ton rapport soigneusement détaillé, ta synthèse de mon sujet de thèse d'un point de vue de mécanicien théorique, ainsi que tes compliments sur mon travail, qui m'ont mis en confiance le jour de ma soutenance et m'ont permis d'aborder la discussion avec un grand enthousiasme.

Professeur Meinhard Kunna, thank you for the very passionate discussion during my defense, I was very glad to finally meet you in person in Madeira and continue our discussion.

Je souhaite également gratifier les efforts continus de Pierre-Oliver Bouchard, à travers les points réguliers que l'on a pu mener avec Hazem. Il faut avouer qu'au vu de la similitude de nos sujets, la qualité des résultats de la thèse d'Hazem fut un élément de pression intimidant, que j'ai finalement pu convertir en source de motivation, grâce à ton accompagnement, afin de donner le meilleur de moi-même et assurer le bon déroulement de ma soutenance.

Nicolas Jacques, merci de t'être déplacé de loin pour assister en présentiel à ma soutenance. Je te remercie également de la discussion scientifique qualitative durant ma soutenance.

Ensuite, je souhaite remercier mon encadrement : En premier, je souhaite adresser un remerciement particulier à la personne ayant le plus contribué à mon travail de thèse, mon père spirituel comme dirait Basile, Jacques Besson. Je te remercie pour les innombrables lec ¸ons que j'ai pu apprendre de toi autant sur le plan scientifique que sur le plan humain. Tu as su parfaitement illustrer l'exemple à suivre d'un chercheur passionné, avec beaucoup de recul, et qui ne cesse d'apprendre pour rester toujours à la frontière de la nouveauté scientifique. J'ai été impressionné par ton humilité et ta capacité de te remettre en question qui sont à mon sens la meilleure qualité qu'espère avoir un chercheur de ton calibre. Je suis très content de pouvoir encore travailler et surtout d'apprendre te toi.

Je souhaite également remercier ma directrice de thèse Sylvia Feld-Payet. Merci de m'avoir confié la continuité de ton travail. Merci de m'avoir consacré des heures innombrables pour que je puisse finalement m'approprier l'avancement de ton travail sur la transition endommagement-rupture. Je te remercie pour tes encouragements pendant la période du Covid qui fut une période compliquée pour moi sur le morale. Et finalement, un grand merci pour tes conseils, remarques ainsi que nos échanges sur mon manuscrit. En plus de m'avoir apporté des corrections sur les différents chapitres, tu me laissais des commentaires très détaillés pour que je puisse apprendre de mes erreurs.

Toujours dans mon encadrement, je souhaite remercier Nikolay Osipov et Stéphane Quilici. Je vous remercie pour le support technique vous m'avez apporté, sans lequel, je n'aurai pu apprendre à naviguer dans le monde très vaste de Z-set. Je vous remercie pour votre exigence sur la robustesse de mon code, qui m'a permis de m'améliorer en tant qu'informaticien et de produire un code propre et efficace. Je vous remercie aussi de m'avoir accueillie dans l'équipe Z-set et de m'avoir donné la chance de continuer mon aventure avec Z-Ductile.

Un grand merci à la personne qui m'a recommandé pour cette thèse et qui sans elle, je n'aurais pas eu à rédiger ces remerciements, Basile Marchand. Merci d'avoir cru en moi et d'avoir contribué à la réalisation de mon rêve de devenir un mécanicien numérique. Je te remercie pour tout le temps que tu m'as accordé et les bugs que tu m'as aidé à résoudre.

Ensuite je souhaite remercier mes amis du CDM que sans eux cette expérience de thèse aurait été beaucoup moins marrante. Je souhaite adresser un merci tout particulier à Hugo L. , mon compagnon de ce chemin de thèse. Merci pour ta disponibilité et pour tes retours toujours enrichissants. Nous avons tellement échangé sur nos deux sujets, au point que j'ai l'impression d'avoir soutenu deux thèses.

Un grand merci à Benjamin de m'avoir intégré et accueilli dans la "Dream team" et pour l'organisation parfaite de chaque petit ou grand événement, Romain de m'avoir montré le chemin quand j'étais bien perdu au fameux lac, Adrien de m'avoir montré que je ne suis pas Chapter 4 Damage crack transition using mesh adaption Contents

aussi débile que ce que je pensais, Hugo D. pour tous les repas du midi très philosophiques qui se sont transformé en Zoom philosophiques pendant le confinement, Glwadys pour les super parties de volley et cours de natation, Luc de m'avoir accueilli dans son académie, Chiraz , Mohamed A. et Robin et tous ceux que j'oublie pour les super moments passés ensemble.

Je souhaite remercier mes co-Bureau du B106 Mélanie et Théophile des discussions toujours passionnantes et des moments de convivialité qu'on a passé ensemble. Je vous remercie aussi de m'avoir fait honneur d'être le parrain de la petite terreur Ulysse.

Un grand merci tout particulier à mon meilleur ami Houssam qui a toujours été à l'écoute de mes problèmes de thèse qui pouvaient sembler très abstraits et pour ton aide à l'organisation pour la soutenance. Merci à Oussama d'avoir toujours cru en moi, Youssef pour toutes les heures au téléphone pendant le confinement qui m'ont permis de voir la lumière au bout du tunnel et Fadl de m'avoir encouragé à entamer ce chemin de thèse ainsi que de m'avoir épaulé tout au long avec tes retour d'expériences précieux.

Pour ma famille, toutes mes gratitudes aux deux personnes qui m'ont amené sur cette terre. Merci papa de m'avoir donné le gout des mathématiques et de la science depuis mon jeunes âge. Cela m'a permis de garder une soif d'apprendre qui ne diminue pas quel que soit la difficulté du problème. Maman, merci de m'avoir appris le sens de l'équilibre dans la vie et surtout de m'avoir donnée confiance en moi pour que je puisse être à la hauteur ce défi. Je souhaite remercier ma petite soeur et frères Mohamed, Malak et Ali pour tous leurs encouragements et admiration qui m'ont permis de garder la tête haute quand le défis de la thèse était à son paroxysme. Un grand merci à mon grand frère Mokhless pour tous ses conseils, soutiens et de m'avoir toujours rappeler que la grandeur se mesure par rapport à un point de départ relatif à chaque personne. Un énorme merci à ma belle famille, Amine, Mohamed, Zineb et Leila qui ont voyagé de si loin pour être présents pour moi à ma grande journée.

Enfin, Kenza, l'amour de ma vie, merci pour m'avoir supporté, encouragé, écouté et rassuré tout au long du chemin, surtout quand je n'avais plus de confiance en moi. Merci de m'avoir offert le meilleur voyage juste après la fin de rédaction sans oublier le meilleur cadeau de thèse. (1) , n (2) , n (3) ; (b) a prolate void; and (c) anoblate void ( from [START_REF] Cao | A model for ductile damage prediction at low stress triaxialities incorporating void shape change and void rotation[END_REF]). . . . . . . . . . . . . . . . . . . . . 2.4 Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for mesh orientations in the case of the axisymmetric specimen. The global responses are similar for all the orientations. However, the total porosity for each mesh size at total failure shows different a crack path for each mesh orientation, which indicates mesh dependence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Geometry of a localization band from [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] . . . . . . . . . . 2.6 Typical nonlocal weighting functions: Gaussian (a) and bell-shaped polynomial functions (b) (from [START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF]). . . . . . . . . . . . . 2.7 Distribution of ductile damage inside the truncated side-pressed cylinder from [START_REF] Saanouni | On the numerical prediction of the ductile fracture in metal forming[END_REF]. The crack is represented in the simulation using the classical "remove elements" method. . . . . . . . . . . . . . . . . . . . . . 2.8 Evolution of damage and mesh at the center of the specimen from [START_REF] El | An anisotropic mesh adaptation strategy for damage and failure in ductile materials[END_REF]. The crack is represented in the simulation using kill elements coupled with remeshing. . . . . . . . . . . . . . . . . . 2.9 Crack surface and damaged elements: (a) crack surface; (b) damaged elements from [START_REF] Yun | A computational methodology for simulating quasi-brittle fracture problems[END_REF]. Controlled kill element was used to represent the crack in a fracture simulation of a quasi-brittle material. . . . . 2.10 Scheme displaying a) intra-element cracking and b) inter-element cracking. . 2.11 Damage distribution reported on the deformed configuration from [START_REF] Cuvilliez | A finite element approach coupling a continuous gradient damage model and a cohesive zone model within the framework of quasi-brittle failure[END_REF]. Intra-element cracking coupled with a cohesive band model inserted on a pre-defined planar crack path. . . . . . . . . . . . . . . . . . . 2.12 Damage distribution reported on the deformed configuration from [START_REF] Feld-Payet | Amorc ¸age et propagation de fissures dans les milieux ductiles non locaux[END_REF]. Intra-element cracking with adaptive remeshing for crack insertion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Crack on a uniform mesh (left) and on a non-uniform mesh (right) from [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF]. The circled nodes are enriched by the jump function whereas the squared nodes are enriched. . . . . . . . . . . . . . . . . . . . . 2.14 Scheme of the global method algorithm from [START_REF] Oliver | Continuum approach to the numerical simulation of material failure in concrete[END_REF] . . . . . . 2.15 Schematic distribution of the function µ in a 2D problem, where the crack path Γ is defined as the zero level set of µ [START_REF] Dias | Strain-injection and crack-path field techniques for 3d crack-propagation modelling in quasi-brittle materials[END_REF] . . . . . . . . . 2.16 polar grid search around the crack front in 2D from [START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF] 2.17 Schematic representation from [START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF], of the damage profile plotted for two different evaluation radii (left) and the corresponding circles centered at the point of maximum damage, thus already on the ridge (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 Schematic representation for of the multiple polar grid search for the maximum from [START_REF] Brokken | Numerical modelling of the metal blanking process[END_REF]] . . . . . . . . . . . . . . . . . . . . . 2.19 Scheme of the search by projection on the orthogonal line to the crack path procedure from [START_REF] Brokken | Numerical modelling of the metal blanking process[END_REF]] . . . . . . . . . . . . . . . . . . . . . 2.20 Scheme of the search by projection on a circle around the the crack path procedure from [START_REF] Seabra | Damage driven crack initiation and propagation in ductile metals using xfem[END_REF] [START_REF] Javani | Three-dimensional finite element modeling of ductile crack initiation and propagation[END_REF] . . . . . . . . . . . . 51 2.23 3D crack initiation scheme using the marching ridges in 3 orthogonal planes from [START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF] . . . . . . . . . . . . . . . . . . . . . . . . . 51

List of Figures

2.24 (a) Visualization of the optimal sampling points of a function (b) and its gradient for 1D (quadratic ) (from [START_REF] Zienkiewicz | The finite element method[END_REF]). . . . . . 57

2.25 Order of super-converging sampling points for different C 0 elements, (from [START_REF] Zienkiewicz | The finite element method[END_REF] 3.1 Dimensions and boundary conditions for the tensile test simulations of the axisymmetric specimen 3.1a and the plane strain specimen 3.1b. . . . . . . . 74 3.2 Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for the axisymmetric specimen (left) and the thickness reduction -∆e/e 0 for the plane strain specimen (right). The comparison with global curves from the literature (respectively from [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] and [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]) indicates a correct parameter fit. Total porosity is displayed at total failure for both simulations: a cup-cone failure (left) and a slant failure (right) can be respectively observed. . . . . . . . . . 75 3.3 Evolution, for the nonlocal model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for the axisymmetric specimen with a central spatial discretization of N h = 90. The confrontation with global curves from the litterature [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] and the local model indicates a correct parameter fit. Total porosity is displayed at various times of the simulation to illustrate the formation of the expected cup-cone crack path. . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.4 Evolution, for the nonlocal model, of the engineering stress as a function of the thickness reduction -∆e/e 0 for the plane strain specimen with a central spatial discretization of N h = 90. The confrontation with global curves from the literature [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF][START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] and the local model indicates a correct parameter fit. Total porosity is displayed at various steps of the simulation to illustrate the propagation of a flat crack instead of a slanted one. This illustrates the need to properly calibrate the nonlocal characteristic lengths. . . . . . . . . . . . . . . . . . . 77 3.5 Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for different spatial discretization in the case of the axisymmetric specimen. The similarity of the global responses does not enable to observe a dependence on the mesh size. The total porosity for each mesh size at total failure shows a cup-cone failure with a slight deviation near the axis of symmetry. . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.6 Evolution, for the nonlocal model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for different spatial discretization in the case of the axisymmetric specimen. The global responses are similar for the mesh sizes N h > 20. The total porosity for each mesh size at total failure shows a cup-cone failure with a straight first horizontal part except for the coarsest mesh size N h = 20, which exhibits a flat crack path. The finest mesh size N h = 180 enables the crack path to bifurcate twice, which can also be experimentally observed. These results indicate mesh independence. . . . . . 79 3.7 Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for various mesh orientations in the case of the axisymmetric specimen. The global responses are similar for all the orientations. However, the total porosity for each mesh size at total failure shows a different crack path for each mesh orientation, which indicates mesh dependence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Evolution

, for the nonlocal model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for various mesh orientations in the case of the axisymmetric specimen. The global responses are similar for all orientations. Besides, the total porosity for each mesh size at total failure shows an identical cup-cone crack path for each mesh orientation. These observations indicate mesh independence. . . . . . . . . . . . . . . . . . . . 3.11 Evolution of the measured band width for both l κ = l ω = 40µm and l κ = l ω = 80µm on both the axisymmetric and plane strain tensile tests as a function of the element height in the current configuration H r1 . There appears to be a linear relation between the band width and the element height. This means that it is possible to estimate a value of the band width l ∞ b free from any discretization error for H r1 = 0. . . . . . . . . . . . . . . . . . . .

3.

12 Evolution of the estimated band width free from discretization error for both the axisymmetric and plane strain tensile tests as a function of the characteristic length (l κ = l ω ). There appears to be a proportional relation between the band width and the characteristic length, with slightly different coefficients for each geometry. This means that it is possible to estimate the size the band width for a given characteristic length and a given geometry. . .

3.13

Contours of total porosity at total fracture for different characteristic lengths l κ and l ω for the axisymmetric specimen. Cup-cone crack paths are observed only for l κ < 80µm. The characteristic length l κ should then be relatively small compared to the specimen size to observe a cup-cone fracture. The larger the characteristic lengths, the wider the localization band (which is linked to the displacement at failure) and the further from the symmetry axis the bifurcation.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.14 Contours of total porosity at total fracture for different characteristic lengths l κ and l ω for the plane strain case. Slant fracture is observed only for l κ = 20µm. The effect of l ω appears to be limited in this case. . . . . . . . . . . .

4.1

Flowchart describing the proposed continuous-discontinuous approach with the different involved components and their inputs/outputs. After resolution of the mechanical problem, the considered scalar field related to material degradation f is smoothed thanks to the superconvergent patch recovery method. Then a crack path tracking algorithm is used to determine where the next crack increment should be inserted and its geometry. The crack increment is not inserted until the insertion criterion is verified. Only then is a new mesh generated thanks to the mesh generator meshGems based, on the one hand, on the intersection of the structure's mesh and the crack increment mesh, and on the other hand on the tracking of the active process zones which are meant to be refined. The variables are then transferred onto the new mesh. To facilitate the return to equilibrium and ensure a satisfactory convergence rate, a reequilibrium step is proposed. Computation of the next load increment can then be performed. . . . . . . . . . . . . . . . . . . . .

4.2

Ellipse shape representing the first crack increment build from a cluster of points having reached the critical value f crit . The ellipse (in red) is defined, from its center G corresponding to the barycenter of the broken points, so that it contains all the projections on its surface of the considered points (in blue). . 4.6 Evolution of the engineering stress as a function of the diameter reduction -∆φ 0 /φ for different minimal mesh sizes obtained when simulating a tensile test on an axisymmetric specimen. With minimal mesh sizes h min ∈ {10, 20, 40}µm, the proposed continuous-discontinuous method has been used whereas the simulation with a very fine mesh is without any mesh adaption or crack insertion. For all the considered meshes, the global responses are similar and the total porosity at complete failure shows a cupcone failure. This similarity indicates that the proposed procedure is rather robust to the choice of the mesh size. . . . . . . . . . . . . . . . . . . . . . 98 4.7 Evolution of the engineering stress as a function of the diameter reduction -∆φ 0 /φ for different crack increment lengths ∆a, using the proposed continuous-discontinuous method, in the case of the axisymmetric specimen. It appears that ∆a has no effect on the global response. However, the largest increment length, ∆a = 180µm, leads to a cup-cone crack path which is not completely smooth compared to the other crack increment lengths. . . . . . . 99

4.8 Evolution of the number of nodes as a function of the number of mesh adaption procedures for the proposed procedure with refinement of only the active process zones (in blue), mesh adaption with only refinement without coarsening of the process zones (in orange) and a fixed mesh (in green). The firt approache consistently leads to a smaller and rather constant number of elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.9 Evolution of the cumulated number of time increments as a function of imposed vertical displacement for the simulation of the axisymmetric tensile test on a fixed mesh (in black), with mesh adaption but without reequilibrium (in red) and with mesh adaption and the proposed reequilibrium procedure (in dark red). The higher number of time increments is obtained for the simulation with mesh adaption without reequilibrium. With the reequilibrium procedure, the total number of increment is approximately divided by two.

The lowest number of time increments is obtained for the fixed mesh for which there is no discrete crack represented, so no transfer needed. In all cases, the crack propagation phase is accompanied by an increase in the number of increments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.10 Specimens' geometry specifications and boundary conditions. . . . . . . . . 102 4.11 Evolution of the engineering stress as a function of diameter reduction -∆φ 0 /φ for all the considered axisymmetric specimens, using the proposed continuous-discontinuous approach. The final diameter reduction decreases with the notch radius. For each global response, three particular instants of the simulation are highlighted with a numbered arrow in order to track damage evolution and crack path in figure 4.12. . . . . . . . . . . . . . . . . 102 4.12 Total porosity f t contours at the three instants of the simulation highlighted in figure 4.11, for all the considered axysimmetric specimens. For all the specimens, there is a cup-cone crack path at stage 3. The diameter corresponding to the bifurcation increases as the notch radius decreases. Let us note that, for the round bar specimen, striction causes the boundary to appear as curved at these advanced stages of material degradation. . . . . . . 4.13 Evolution of the engineering stress as a function of the thickness reduction -∆e 0 /e 0 of the plane strain specimen for both a simulation with the proposed continuous-discontinuous procedure and a simulation with a fixed mesh. The global curves are similar up to the latest stages of the simulation where the simulation with the proposed continuous-discontinuous procedure ends up with a larger final thickness reduction. Total porosity is displayed at various times of the simulation using the proposed method, and at total failure for the fixed mesh simulation. The proposed method predicts an axisymmetric slant crack path with a bifurcation whereas the fixed mesh simulation predicts a symmetric slant crack path. The bifurcation coincide with the deviation of the global responses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Evolution of the engineering stress as a function of Crack MOuth Displacement (CMOD) for two meshes of the CT specimen with different minimum element sizes: h min ∈ {20, 40}µm. The global response is identical for both meshes which tends to indicate mesh convergence. Contours of the total porosity f t are displayed at different times of the simulation. Before crack initiation (marked with a green cross), crack tip blunting is well described thanks to mesh adaption which enables to preserve good quality elements near the crack tip. Then the crack propagates over a relatively long distance with still well shaped fine elements near the crack tip and a coarsened mesh along the crack lips. . . . . . . . . . . . . . . . . . . 4.15 Evolution of the engineering stress as a function of diameter reduction -∆φ 0 /φ for the simulation of a bar with round section with the proposed approach in 3D (in red), the proposed approach with 2D axisymmetric elements and a smaller minimum element size (in green) and a fixed mesh with 2D axisymmetric elements (in blue). All the global responses are similar, which thend to indicate the validity of the used 3D mesh. Contours of the total porosity are displayed at different times for the 3D simulation to illustrate the cup-cone crack propagation. . . . . . . . . . . . . . . . . . . . 4.16 a: initial mesh of quarter of the specimen for the 3D simulation. b: fracture surfaces obtained from the 3D simulation. The predicted cup-cone crack surfaces are smooth, well defined and qualitatively in agreement with what can be expected from the experimental results obtained by [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.

1 Microstructure of DP450 steel after nital etching from [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF]. 5.2 Specimen used in the experimental campaign(measurement points are marked with •). Specific dimensions used at high rates are outlined in red (from [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Experimental setups and measure devices for quasi-static tensile tests. (from [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Illustration of shadow tracking measure on NT1 specimen. (from [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Cylinder subjected to internal pressure. Error estimator using the recovered stress σ * 22 (from [START_REF] Ródenas | Improvement of the superconvergent patch recovery technique by the use of constraint equations: the spr-c technique[END_REF]). . . . . . . . . . . . . . . . . . . . . 5.6 Scalar incremental strategy adapted for thin 3D meshes. . . . . . . . . . . . . 5.7 For the standard flat specimen, the load-Engineering transverse strain obtained from simulation and experiments are compared. Also, the map of total porosity f t is displayed at complete failure. . . . . . . . . . . . . . . . . 5.8 For the NT1 specimen, the load-notch opening is obtained from simulation and experiment are compared. Also, the map of total porosity f t is displayed at complete failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 For the NT2 specimen, the load-notch opening is obtained from simulation and experiment are compared. Also, the map of total porosity f t is displayed at complete failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.10 For the V45 specimen, the load-notch opening is obtained from simulation and experiment are compared. Also, the map of total porosity f t is displayed at complete failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.11 For the M-shape specimen, the load-relative displacement curve ∆L is obtained from simulation and experiment are compared. Also, the map of total porosity f t is displayed at complete failure. . . . . . . . . . . . . . . . . 122 5.12 For the FN specimen, the load-displacement obtained from simulation and experiment are compared. Also, the map of total porosity f t is displayed at complete failure. . The numerical developments of this work are implemented in the finite element software Z-set, dedicated to general structural analysis, with a particular emphasis on advanced material representation. The software is composed of different modules that allow performing pre/post-processing operations in addition to the finite element analysis of a structure.

Industrial context

Predicting the development of damage and tearing is a major challenge when trying to prevent failure. This is particularly important to prevent crack propagation in critical parts of an aircraft or a pipeline, or to control the propagation of cracks, for example, in the case of the metal part forming. More broadly, the control of damage phenomena is part of most design phases, when a technical and economic trade-off between structural integrity and efficiency is required.

These predictions increasingly rely on numerical simulation tools. These numerical simulations are generally based on the finite element method in the context of fracture mechanics. In particular, these methods are essential in the case of complex industrial parts for which full-scale experiments are too expensive or impractical.

Developing a simulation tool capable of predicting the structural failure of complex industrial parts is one of the objectives of Transvalor to broaden the use of numerical simulations to structural failure, and thus propose a complete simulation device ready to be used by technicians and engineers.

Nowadays, the modeling and numerical simulation of ductile fracture for industrial structures remain complex due to finite deformations, softening phenomena, and the need for more accurate characterization of materials. In this industrial context, this thesis aims at producing a framework for the simulation of ductile fracture compatible with actual industrial simulations.

Scientific context

Several approaches exist in the literature to describe the development of damage. The first approach, referred to as "global approach", is based on macroscopic parameters such as the energy release rate (G) or stress intensity factors (K I ). Additional parameters describing the level of constraint at the crack tip, such as the T and Q parameters, can also be used [START_REF] O'dowd | Family of crack-tip fields characterized by a triaxiality parameter-I. Structure of Fields[END_REF][START_REF] Thaulow | Constraint correction of high strength steel. Selection of test specimens and application of direct calculations[END_REF]. The global approaches introduce a reduced number of parameters, which explains their popularity for industrial applications. However, they have some limitations: (i) they are limited to cases where an existing crack is present and are unable to predict crack initiation, (ii) they are limited to cases where plasticity is limited, (iii) material parameters determined using this approach (e.g. toughness (K J ) or J-∆a curves) may strongly depend on the geometry of the investigated test specimens [START_REF] Sumpter | Constraint based analysis of shallow cracks in mild steel[END_REF][START_REF] Sumpter | An experimental investigation of the T stress approach[END_REF][START_REF] O'dowd | The role of geometry can crack growth on constraint and implications for ductile/brittle fracture[END_REF] so that transferability to structures remains difficult, (iv) the approach does not explicitely accounts for physical degradation mechanisms.

A second approach, referred to as "local approach to failure", is based on a micromechanical analysis of the material degradation. In the case of ductile failure of metals, micrographic investigations have evidenced different stages of damage [Garrison and Moody, 1987b]: chronologically, the nucleation of voids (creation of porosity), the growth of voids, and the coalescence of voids. Different models, based on micromechanical consideration [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][START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF] or purely phenomenological [START_REF] Lemaitre | A continuous damage mechanics model for ductile fracture[END_REF], are already available in the literature. All these damage models share a common characteristic nevertheless: softening or negative strain hardening. Without special treatment, this behavior leads to a mathematically ill-posed and non-physical problem since the dissipated energy at failure cannot be determined [START_REF] Benallal | Bifurcation and Localization in Rate-Independent Materials. Some General Considerations[END_REF]. In practice, this leads to pathologically dependent results on spatial discretization. The modeling of the fracture phenomena must be enriched by coupling the responses between neighboring material points; this is called nonlocal modeling [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF][START_REF] Lorentz | A variational formulation for nonlocal damage models[END_REF][START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]. The main advantage of continuous damage models based on the local approach is that they allow for crack initiation without the need to introduce a pre-crack, unlike the global approach. However, this continuous description of fracture is only valid until the initiation of a crack. After crack initiation, the continuous model cannot represent displacement jumps across the crack.

A third approach, called the continuous-discontinuous approach to failure, consists of using a continuous damage model up to crack initiation. Then, when a zone is sufficiently damaged, a discrete crack representation is introduced into the model to represent the displacement jumps across the crack lips accurately. Crack propagation is represented similarly. In practice, this approach requires the choices of several elements. How to represent the crack? Where to represent the crack? When to switch from a continuous representation of the crack (smeared manner) to a discrete representation of the crack? A multitude of combinations is then possible, giving rise to as many different procedures.

Objectives

In this work, this third approach is chosen to simulate ductile fracture. The main objective of this work is to develop a framework for the continuous-discontinuous approach, based on a micromechanical model adapted to finite strains for the continuous part and on remeshing techniques for the discontinuous part. The developed approach must be numerically robust and efficient. Its applicability is validated by making comparisons between simulations and experimental tests.

Structure

This work is structured as follows:

In chapter 2, a review of the literature on continuous damage models and micromechanical models for ductile fracture is first presented. The localization problem, which is the major drawback of these coupled damage models in their local form, is explained and then illustrated by a simulation using the local GTN model. At the end of this first part, the GTN model is chosen as the damage model to be used for this work. Then, the regularization approaches of the damage models are presented with their strengths and weaknesses. After this, the implicit gradient method is chosen as the regularization method for the GTN model. Finally, the different elements of damage-tocrack transition approaches are detailed. The proposed continuous-discontinuous approach based on remeshing techniques is chosen to simulate ductile failure. At the end of the first chapter, the framework for the simulation of ductile fracture is detailed. The existing components of the framework are presented, and the components that have been developed during this work are associated with reference to the corresponding chapter detailing them.

In chapter 2, a continuous nonlocal damage model based on the implicit gradient methodology is developed. The model uses two characteristic lengths in order to regularize each damage mechanism independently: (i) growth of primary voids, (ii) nucleation of secondary voids on strengthening particles. The model is written within a finite strain framework, with the updated Lagrangian formulation in order to be compatible with the chosen remeshing strategy. The implementations of the nonlocal finite element and the GTN damage model are detailed. A detailed mesh convergence study is conducted to verify mesh insensitivity to size and orientation, which is essential in remeshing. In addition, with a convergence study on the localization bandwidth, a relationship between the characteristic length and the localization bandwidth is identified. This relationship is used in order to determine the minimal element size to achieve mesh convergence. At the end of this chapter, applications to cup-cone and slant failure are presented using the developed model. The effect of the specimen size on the crack paths is demonstrated.

In chapter 3, the continuous-discontinuous approach is presented with a pragmatic remeshing strategy to optimize mesh size and, therefore, computation time. At first, new components to improve the simulation of ductile crack propagation using 3D meshes are developed. In particular, a novel crack initiation algorithm is developed in order to initiate correct shaped cracks. In addition, a new criterion for the continuous-discontinuous transition is developed in order to determine the crack increment size to be inserted. The problem of the deterioration of the convergence rate after remeshing is addressed, and a reequilibrium method is then proposed to overcome this issue. The efficiency of the developed continuous-discontinuous approach is tested with respect to spatial discretization and time discretization. Finally, applications using the developed method using different 2D and 3D meshes of specimens used in experimental tests in the literature are presented.

In chapter 4, the main objective is to evaluate the predictive capability of the proposed approach. For that purpose, an experimental database from a previous thesis [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF] is used to perform comparisons between simulations and experimental tests. At first, the material of the experimental database is briefly recalled. Then, the constitutive model is adapted in order to reproduce the nonlinear phenomena observed in the experiments for quasistatic loading. The nonlocal framework with two characteristic lengths from chapter 2 is used. Finally, simulations of the experimental tests using the developed continuous-discontinuous approach are compared to the experimental result. A good general agreement is found.

Finally, chapter 5 summarizes this work and proposes some avenues for future research.

Chapter 2

Toward an efficient continuous-discontinuous strategy for the simulation of ductile failure: Review In this work, the main focus is ductile fracture, which is a mode of material failure that occurs generally at high strain levels. The fracture surface is characterize by the presence of dimples ; inside these dimples second-phase particles or inclusions are often observed [Garrison and Moody, 1987a].

Résumé en franc ¸ais

Mechanisms of ductile failure

The micrographic analysis of metals and the observation of the surface or of the bulk of the deformed material containing inclusions has led to identify the process of ductile failure as a three stages process [START_REF] Besson | Damage of ductile materials deforming under multiple plastic or viscoplastic mechanisms[END_REF]:

1. Voids nucleate at inclusion and second-phase particles by cohesion from the metallic matrix or by fracture of these inclusion. Voids may also preexist in the material, 2. Due to large plastic strain and triaxial stresses, voids expand and their shape is determined by material properties and test conditions, 3. The final step of the three stage process of ductile fracture is coalescence. During this stage, large enough voids link together, either by internal necking of the ligament or by void sheeting.

In general, the ductile fracture process is accompanied by large plastic deformations. It is therefore mandatory, in the numerical simulation of this process by the finite element method (FEM), to use a large deformation framework.

Continuum damage models

Continuum damage model aim to characterize, represent and model at the macroscopic scale the effects of distributed defects and their growth on the material behavior. In order to model the evolution of defects, a damage variable, generally called D as described in the work of [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF], is introduced. In figure 2.1, a Representative Volume Element (RVE) of section S on which is applied a normal load F is illustrated: where the stress in the undamaged volume is σ = F S . It can be observed that the effective area of the RVE subjected to uniaxial tension is S -S D , where S D represents the defects area in the considered plane. The damage variable D is defined as:

D = S D S (2.1)
σ = F S -S D = F/S 1 -S/S D = σ 1 -D > σ (2.2)
The strain equivalent principle ε = ε as introduced in [START_REF] Lemaitre | Mécanique des Matériaux Solides[END_REF]:

σ Ẽ = σ E (2.3)
Where Ẽ = (1 -D)Eis the effective stiffness. This effective stiffness is used numerically to describe the loss of load crying capacity due to damage. The continuum damage models are mostly defined in the framework of thermodynamic of irreversible processes. In this framwork, state variable are introduced to define the thermodynamic state of the RVE. These variables are split into "measurable" state or external variables (as total strain tensor ε and temperature T ), and internal state variables that are not directly "measurable" (as the plastic strain p, damage D, ...).

For each stat variable an thermodynamic force is associated: σ is associated to the elastic strain ε, the hardening-law radius R is associated to the cumulative plastic strain p and the elastic release rate Y is associated with the damage variable D. Then, to link these force to the state variable, two states function are defined:

1. The free energy density per unit volume ψ which is expressed as a function of the state variables. The function ψ is divided in case of an elastoplastic material wit damage, into an elastic and plastic contributions as:

ψ(ε e , D, p) = ψ e (ε e , D) + ψ p (p) (2.4)
where ε e is the elastic strain tensor.

The thermodynamic forces are then obtain by differentiating ψ with respect to the stat variables:

σ = ∂ψ e ∂ε , R = ∂ψ p ∂p , Y = ∂ψ e ∂D (2.5)
2. The dissipation potential D that considers the contribution plasticity and damage irreversible processes. In case of standard generalized elastoplastic materials, the dissipation potential is expressed using the yield function Φ(σ) and the damage dissipation potential F D as:

D = Φ(σ, R) + F D (Y, D) (2.6)
The evolution of the internal state variables is then obtained by differentiating D with respect to thermodynamic forces:

εp = ε -ε e = λ ∂Φ ∂σ , ṗ = - λ ∂Φ ∂R , Ḋ = - λ∂F D ∂Y (2.7)
where λ is the plastic multiplier.

Damage evolution depends of the choice of the damage potential F D . The most used potential is the one introduced in [START_REF] Lemaitre | A course on damage mechanics[END_REF]:

F D (Y ) = S 0 (b + 1)(1 -D) -Y S 0 b+1 (2.8)
However, other formulations have been presented in literature:

1. [START_REF] Tai | A new microvoid-damage model for ductile fracture[END_REF] :

F D (Y ) = S 0 D 2 -Y S 0 2
(2.9)

2. [START_REF] Chandrakanth | An isotropic damage model for ductile material[END_REF] :

F D (Y ) = 1 2 S 0 (Dp) α/n -Y S 0 2 (2.10)
3. [START_REF] Bonora | A nonlinear CDM model for ductile failure[END_REF] :

F D (Y ) = 1 2 S 0 (D c -D) (α-1)/α (1 -D)p (2+n)/n -Y S 0 2 (2.11)
4. [START_REF] Bouchard | An enhanced lemaitre model formulation for materials processing damage computation[END_REF]]

F D (Y ) = S 0 b + 1 1 1 -D -Y S 0 b+1 1 p α (2.12)
where S 0 , D c , b, α and n are damage materials parameters. The damage evolution Ḋ for each potential can be deduced using eq.2.7.

Micromechanical models for ductile failure

Another type of models exist in literature called the micromechanical damage models. These type of models represent ductile failure based the direct description of the ductile damage process at micro scale: the nucleation, growth and coalescence of cavities. A brief review of these models are presented in this section.

Rice & Tracey Model

The first attempts to model voids growth are the theoretical analysis by McClintock [Mc Clintock, 1968] and [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF]. The first model by McClintock computes the growth of cylindrical voids in a plastic material. The second model by Rice & Tracey studies the growth of a spherical void in a perfectly plastic material. From both analysis, quantitatively similar formulations of voids growth are found:

Ṙ R = A + B exp C σ ∞ m σ ∞ eq ε∞ eq (2.13)
where R is the radius of the cavity, ε∞ eq is the von Mises equivalent strain imposed at a boundary considered infinitively far from the void. σ ∞ m and σ ∞ eq are the remote mean and von Mises stresses. A, B, etC are the constants of the model. A is close to one and depend of loading direction, C equal to 1.5 in the analysis of Rice & Tracey, in the McClintock model (κ as the hardening coefficient). At last B is a numerical coefficient not so well defined in their analyses, for Rice & Tracey is equal to 0.283 and To adapt the Rice&Tracey model to the elastoplastic case, [START_REF] Mudry | Etude de la rupture ductile et de la rupture par clivage d'aciers faiblement alliés[END_REF] replaces ε∞ eq by the rate of the equivalent plastic strain ( ṗ∞ ). To account for the hardening of the metallic matrix, σ 0 is substituted by the flow stress σ F which is a function of p ∞ and which is equal to the von Mises stress:

Ṙ R = B exp C σ m σ eq ṗ∞ (2.14)
After integrating the expression in 2.14:

log R R 0 = p pc B exp C σ m σ eq dp (2.15)
where R 0 is the initial void radius and p c is the critical strain at which nucleation from particles starts. Using the model define in 2.15, a fracture criterion can be derived: the fracture of the material begins if the value of R/R 0 exceeds a critical value R/R 0 | c . This value is supposed to be a material property.

The major finding of this model is that the void growth process is a function that depends exponentially on the stress triaxiality which is the third of the trace of the stress tensor ( 13 trace(σ)) divided by the equivalent von Mises stress (σ eq ).

The Rice and Tracey model and its variants can easily be implemented to post-process FEM computations to evaluate the damage state in a structure. This approach was generalized by introducing a damage indicator D. Its evolution can be expressed in a generic manner as (see [START_REF] Defaisse | Ductile fracture of an ultra-high strength steel under low to moderate stress triaxiality[END_REF]):

Ḋ = FUNCTION(stress state, p) ṗ (2.16)
In particular the effect of the third invariant of the stress tensor (i.e. the Lode parameter) is nowadays accounted for [START_REF] Bai | A comparative study of three groups of ductile fracture loci in the 3d space[END_REF]. Although this model type is still used in computations of ductile damage to this day, it should be restricted only to preliminary studies or to model failure initiation only (preferably in uncracked structures). This is because this type of model lacks the ability to describe the interaction between voids and the effect of void growth (i.e. softening) on the mechanical behavior of the material.

The Gurson model

A second analytical model was introduced by [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criteria and flow rules for porous ductile media[END_REF]. This model is a result of an upper bound analysis of a perfectly plastic finite sphere (matrix) containing a spherical void.

Damage is represented in this model by the porosity (f ) which is the ratio of void volume to total volume. The Gurson analysis leads to the definition of a yield plastic surface that depends on the macroscopic stress and the porosity:

Φ = σ 2 eq σ 2 0 + 2f cosh 3 2 σ m σ 0 -1 -f 2 (2.17)
Where σ m is the mean stress, σ 0 is the yield stress For an undamaged material without porosity, the von Mises criterion is retrieved. Failure of the material occurs when the stress tensor is null (i.e no force transmission by the material). In the classical Gurson model, complete failure occurs when the porosity f is equal to 1, meaning the material consists only of voids! The Gurson model has been modified in order to obtain a more realistic failure condition.

The evolution law of porosity can be deduced from volume conservation, with a constant matrix volume:

f = V -V m V ⇔ ḟ = V m V V V = (1 -f )trace( εp ) , Vm = 0 (2.18)
where V m is the volume of the matrix and V is the total volume. From the definition of the yield surface 2.17, the elastic domain shrinks while the porosity increases (see figure 2.2):

By deriving the yield criterion with respect to the stress and using the normality rule, the evolution law of the plastic deformation can be written as: where κ is a plastic multiplier.

εp = (1 -f ) κ ∂Φ ∂σ (2.19)
The model can be also used for material with hardening by substituting in the yield surface the yield stress σ 0 by the flow stress (σ F ) that takes into account work-hardening. A cumulated plastic strain can be defined by the following equation:

(1 -f )σ F ṗ = σ : εp (2.20)
The Gurson model is considered one of the very sound models, because of its strong micro-mechanical background. Nevertheless, this model only describes the growth of cavities, but ignores the other two mechanisms of ductile fracture. Other models, based on the Gurson model, are developed in order to be more consistent with the experimental reality.

The Gurson-Tvergaard-Needleman (GTN) model

In order to fit properly experimental test using the Gurson model, the nucleation and coalescing of porosity needed be accounted for in the model. A porosity of nucleation f n was added to the model following the work of [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF], and briefly recalled in the first part of this section. Also, in the work of [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF], the Gurson yield surface was modified to account for the coalescing of porosity, and the new yield surface is named after Gurson-Tvergaard-Needleman (GTN). This modification is detailed in the second part of this section.

(i) Nucleation To introduce nucleation in the GTN model, the porosity is split into two contributions, the classical growth porosity deduced from the mass conservation and a new contribution that depend on the cumulated plasticity rate κ, as the following:

ḟ = ḟg + ḟn = (1 -f )trace( εp ) + A n κ (2.21)
A Gaussian form for the nucleation rate introduced by [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF], is usually associated with the GTN model:

A n = f N s N √ 2π exp 1 2 ( κ -κ N s N 2 ) (2.22)
Where f N is the fraction of inclusions at which porosities could be nucleated; κ N is the cumulated plastic deformation that correspond to the maximum of nucleation rate; s N the standard deviation on the nucleation strain.

(ii) Coalescence In order to take into account the effect of coalescence in the Gurson model, [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] introduced the following modified yield surface:

Φ = σ 2 eq σ 2 0 + 2q 1 f * cosh( 3q 2 2 σ m σ 0 ) -1 -(q 1 f * ) 2 (2.23)
Where q 1 and q 2 are two material parameters to be identified, and f * is the effective porosity defined by:

f * = f if f < f c f c + δ(f -f c ) if f ≥ f c (2.24)
where f c and δ are respectively the threshold and speed of the coalescence mechanism. Although this particular form for f * is very often used, other functions could be suitable.

Using low values of f c can lead to convergence problems using the FEM. f c can be computed from a unit cell computations (see [START_REF] Zhang | A complete Gurson model approach for ductile fracture[END_REF]).

(iii) Failure Using the GTN model, failure is achieved when the porosity f or effective porosity f * reaches it maximal value:

f = f r -f c δ + f c (2.25)
where f r = 1/q 1 and is the critical porosity. If this porosity is reached and using 2.23, a null stress tensor is obtained σ = 0. These two equation resume the most known model for ductile failure using local approach: the GTN model (iv) An effective stress from the GTN model Micro-mechanics based models can have a generic expression. For the GTN model, this generic presentation can be expressed by defining an effective stress σ * as the following (see [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF]):

   F (σ * , R) = σ * -R G(σ, σ * , f * ) = σ 2 eq σ 2 * + 2q 1 f * cosh( 3q 2 2 σm σ * ) -1 -(q 1 f * ) 2 def. σ * = 0 (2.26)
a stress measure σ * is defined implicitly to satisfy the equation G(σ, σ * , f * ) = 0.

GTN model extensions for low stress triaxiality

One of the major drawbacks of the GTN model concerns the fact that it underestimates damage under shear or, more generally, under low triaxiality loading [START_REF] Li | Ductile fracture: Experiments and computations[END_REF]. In order to improve the performance of the GTN model outside of high triaxiality loading, some extensions take into account the rotation of the cavities in addition to the change of their shapes while others phenomenologically associate the evolution of the damage to the deviatoric part of the shear. The common principle of all these improvements is to involve the third invariant or Lode parameter in the evolution of porosity. In this section, we describe these extensions that incorporate the effect of shear in the GTN model.

(i) Nahshon& Hutchison extension In the work of [START_REF] Nahshon | Modification of the Gurson model for shear failure[END_REF], in addition to the classical growth porosity of the GTN model, a shear porosity is introduced where its evolution is expressed as:

ḟs = kf (1 -L(σ) 2 ) σ dev : εp σ e q (2.27)
where k in a constant between 0 and 3. σ dev is the deviatoric part of the stress tensor. L is the Lode parameter and can be expressed as:

L(σ) = 27det(σ) ( 3 2 σ dev : σ dev ) 3/2
(2.28)

The shear porosity f s is not directly related to the plastic strain as is the case for the growth and nucleation porosity. Thus this f s variable can be seen only as a damage variable. This shear modification was used to simulate successfully ductile failure under low stress triaxlality. However, it effect is non negligible in case of high stress striaxiality loading [START_REF] Nielsen | Effect of a shear modified Gurson model on damage development in a FSW tensile specimen[END_REF]. In order to limit this effect, a slit modification to this shear porosity was proposed in [START_REF] Nielsen | Effect of a shear modified Gurson model on damage development in a FSW tensile specimen[END_REF]:

ḟs = kf (1 -L(σ) 2 )Ω(η)
σ dev : εp σ eq (2.29) where Ω(η) is a function of the stress triaxiality η = σ H /σ eq :

Ω(η) =      0 if η > 1/2 1 -2η if 0 ≤ η ≤ 1/2 1 if η < 0 (2.30) (ii) Xue Extension
In the work of [START_REF] Xue | Constitutive modeling of void shearing effect in ductile fracture of porous materials[END_REF], a similar approach to the previous model where a new shear damage variable was introduced. Based on the analysis of a hollow plate under shear loading, a new damage variable D was introduced where it evolution depends on both the total porosity f and the shear loading:

Ḋ = δ D (q1 ḟ + Ḋs ) (2.31)
where the coefficient δ D is defined as:

δ D = 1 if D ≤ q 1 f c 1/q 1 -fc f -fc if q 1 f c < D ≤ 1 (2.32)
and Ḋs is the evolution of the damage under shear loading: Ḋs = q 3 f q 4 g(θ L )ε p eq εp eq (2.33) where

θ L = arctan 2-σ II -σ I -σ III √ 3(σ I -σ III)
is the Lode angle and g(θ L ) is a function that verifies g(0) = 0 and g(1) = 1. The following expression for this function was proposed in [START_REF] Xue | Constitutive modeling of void shearing effect in ductile fracture of porous materials[END_REF]:

g(θ L ) = 1 - 6 π |θ L | (2.34)
And from the analytic solution of the hollow plate, q 3 = 6 (π), q 4 = 1/2 for 2D simulations and q 3 = 3(6/π) 1 /3, q 4 = 1/3 for 3D simulations. For the implementation of this model, two approach were proposed in [START_REF] Xue | Constitutive modeling of void shearing effect in ductile fracture of porous materials[END_REF]: (i) to consider the damage variable D as an internal variable that influence the integration of the behavior; (ii) as a parameter computed based on the porosity and stress tensor after the integration of the behavior.

(iii) Shape and rotation extensions [START_REF] Kailasam | A general constitutive theory for linear and nonlinear particulate media with microstructure evolution[END_REF]] have performed analyses on the behavior of ellipsoidal cavities in a composite and have shown by second order homogenization that the shape and orientation of the cavities that evolve with the deformation play an important role in the macroscopic behavior of the material. In order to take into account these microscopic phenomena, they introduced new variables describing microscopically the deformation and the rotation of the cavities. The model of [START_REF] Kailasam | A general constitutive theory for linear and nonlinear particulate media with microstructure evolution[END_REF]] was formalized and later completed by [Danas and Ponte-Castaneda, 2009a,b] by the latter work, [START_REF] Cao | A model for ductile damage prediction at low stress triaxialities incorporating void shape change and void rotation[END_REF] recently performed an application to the GTN model and developed a model capable of simulating damage at low triaxiality, by making an ad hoc modification to the GTN law :

Φ = σeq 2 σ 2 0 + 2q 1 f * cosh( 3q 2 2 σm σ 0 ) -1 -(q 1 f * ) 2 (2.35)
where

   σeq = (σ eq (1 -αf ) + 2αf (σ : ((I -1 3 I ⊗ I) : Q) : σ)) 1/2 σm = σ:( 1 3 I⊗I:Q):σ (3|σ:( 1 3 I⊗I:Q):σ|) 1/2 (2.36)
Where Q is a fourth order tensor where its components are the semi-axes, a 1 , a 2 , a 3 , of the general ellipsoidal void and the three orientation vectors n (1) , n (2) , n (3) as illustrated in 2.3. The expression of the tensor Q is detailed in the appendix of [START_REF] Cao | A model for ductile damage prediction at low stress triaxialities incorporating void shape change and void rotation[END_REF].

Unlike the previous two extensions, the evolution of porosity in this model keeps the same expression as in the GTN model. The shear effect is introduced implicitly by the modified threshold function due to the consideration of anisotropy in the stress tensor. From the simulations in [START_REF] Cao | A model for ductile damage prediction at low stress triaxialities incorporating void shape change and void rotation[END_REF], this model leads to physically relevant results under hydrostatic and shear loading. The authors also showed differences in failure behavior under prolate and oblate ellipsoids. But in this method, the expressions resulting from the second order homogenization seem very heavy. Its numerical performance also remains to be compared with the two previous models which are more pragmatic. 1) , n (2) , n (3) ; (b) a prolate void; and (c) anoblate void ( from [START_REF] Cao | A model for ductile damage prediction at low stress triaxialities incorporating void shape change and void rotation[END_REF]).

Damage localization

In the process of ductile failure of a material, a significant accumulation of plastic strain in a thin zone near the failure point is observed in experimental tests. The size of this zone, commonly called the localization band, depends on the material properties, the geometry of the structure and the nature of the loading.

Nevertheless, the numerical simulation by finite elements, with local models, shows a strong dependence on the mesh: the damage tends to be localized in a band formed by a single layer of elements, whatever its size. Moreover, the orientation of the damage band also depends on the orientation given to the mesh. If an infinitely refined mesh is used, the thickness of the localization will tend towards zero, leading to a zero dissipated energy, totally unrealistic. In order to illustrate this mesh dependency, we conducted simulations of an axisymmetric tensile test, for different mesh orientations, using the local GTN model. The load-diameter reduction curves are different for each mesh orientation. Also the crack path follows the mesh orientation. Thus the crack path is also mesh sensitive This problem is mainly due to the instability of the macroscopic behavior of the material (softening) induced by the increase of damage. During the softening of the material, the system of equations describing the mechanical problem becomes non-elliptic, thus, the continuum damage models that are motivated by micromechanical models lose part of their validity.

Localization is assumed to occur when it becomes possible to form a strain rate discontinuity in a planar band. This band is characterized by its unit normal n and the displacement jump across the band whose direction is denoted g as shown in figure 2.5. In [START_REF] Hadamard | Lec ¸ons sur la propagation des ondes et les équations de l'hydrodynamique[END_REF], it was shown that the strain rate tensor depends on n and g as:

∆ ε = 1 2 ( n ⊗ g + g ⊗ n) (2.37)
The normal vector n can be obtained as the vector that minimizes the Rice criterion [START_REF] Rice | The localisation of plastic deformation[END_REF][START_REF] Rice | The mechanics of earthquake rupture[END_REF] which is expressed in case of small strain hypothesis as:

det( n.D. n) = 0 (2.38)
where D in elasto-plastic materials is the elastoplastic tangent matrix ( σ = D : ε). Additional terms to D are needed in case finite strain (in case of hypo-elasticity see [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF], and hyper-elasticity see [START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF]). In [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF], in case of visco-plastic material, the bifurcation analysis is conducted using a perturbation analysis. It was also proposed to use in general the consistent tangent matrix as an approximation of the elastoplastic tangent.

Conclusion

In this section, several local damage models have been presented. These models describe the degradation of the material by means of a spatial field, generally called damage field. The damage field can either: affect the material behavior, which is the case of coupled damage models; or not, which is the case of uncoupled damage models. The use of coupled damage The global responses are similar for all the orientations. However, the total porosity for each mesh size at total failure shows different a crack path for each mesh orientation, which indicates mesh dependence.

Figure 2.5: Geometry of a localization band from [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] models leads to softening of the material at critical damage values. The high damage field is often localized in narrow bands that depend on the discretization, in other words, local coupled models suffer from pathological mesh dependency.

From a physical point of view, this localization phenomenon is explained by a lack of information on the scale in the modeling. Indeed, the construction of a macroscopic behavior law in local models is based on the assumption that the variation of macroscopic mechanical fields is negligible below a minimal scale characterizing the Representative Volume Element (RVE). This assumption is no longer respected if localization bands develop below this scale. It is therefore necessary, in order to remain in the spirit of the continuous medium, to introduce a notion of characteristic length corresponding to the scale of REV.

Regularized damage models

Different strategies exist in the literature to overcome the pathological mesh dependency when using continuum damage mechanics models (CDMs). There is cohesive zone models, thick level set method and finally the nonlocal methods. In this section only a brief review of nonlocal models will be presented. This type of models consist in introducing a material internal length. With the use of nonlocal models, the constitutive equations are preserved with small adaptations. This is beneficial especially in case of constitutive models based on a strong micro-mechanical description of the damage processes.

Integral methods

The underlying idea of integral methods is to take into account a spatial vicinity effect to describe the behavior of a material point. This is achieved by computing the weighted average of a local variable around the material point under consideration [Pijaudier-Cabot and[START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Bazant | Non local continuum damage. localization, instability and convergence[END_REF]. Consider D a damage variables causing localization and D its nonlocal counterpart averaged at a space position x in the body Ω:

D(x) = Ω D(y)w(y -x)dΩ Ω w(y -x)dΩ (2.39)
Where w is a weight function, verifying the following condition:

V w(x)dΩ = 1 (2.40)
An example of this weight function can be expressed as Gaussian:

w(x) = 1 (2π) 3/2 l 3 exp - ||x|| 2 2l 2 (2.41)
Another example of the weight function can be expressed as bell-shaped polynomial:

w(x) = 105 32π 2 l 3   1 - ||x|| l 2   2 (2.42)
Where l is the characteristic length. The bell polynomial function 2.42 is often preferred in [START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF]).

numerical computations, because it ensures null interactions between any pair of material points such that their distance is larger than l (i.e., w(x) = 0 / ||x||/l ≥ 1). And subsequently, the nonlocal field is locally averaged over a small sub domain. Implementing this integral method in the finite element method can be challenging. In particular, the calculation of the integral 2.39 requires the consideration of a patch of elements around the node under consideration. This can be very intrusive in the FEM code unless a fully explicit scheme is used. In that case the average is performed over the increments of the damage variable The integral method was used to regularize GTN model in [START_REF] Enakoutsa | Numerical implementation and assessment of a phenomenological nonlocal model of ductile rupture[END_REF]. The integral method considers that any sufficiently close node must interact. However, nodes that are sufficiently close but separated by a crack, a geometric hole or even by a localization damage band should interact less. One solution consist in modeling the damage evolution as a wave propagation in continuous media based on the Wentzel-Kramers-Brillouin (WKB) approximation. In this formulation, the interaction distances are defined as the solution of a stationary Eikonal equation dependent on the damage variable. It allows to model nonlocal interactions that progressively disappear in the damaged areas as the wave is blocked by cracks or geometric holes [START_REF] Desmorat | Nonlocal models with damage-dependent interactions motivated by internal time[END_REF][START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF].

Explicit/Implicit gradient methods

In order to overcome the complex and intrusive implementation of integral methods, gradient methods were developed [START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF]. This method was initially proposed to regularize the local equivalent deformation ε eq but was later extended to deal with any type of local variable D which can be a scalar or a tensor. The Taylor series expansion of a local variable D around a spatial position x i is as follows:

D(y) = D(x) + ∂D ∂x i (y i -x i ) + 1 2! ∂ 2 D ∂x i ∂x j (y i -x i )(y j -x j ) (2.43) + 1 3! ∂ 3 D ∂x i ∂x j ∂x k (y i -x i )(y j -x j )(y k -x k ) + 1 4! ∂ 4 D ∂x i ∂x j ∂x k ∂x l (y i -x i )(y j -x j )(y k -x k )(y l -x l ) ...
The substitution of 2.43 into 2.39, and considering the integral over an infinite body (R 3 ) leads to the following:

D(x) = D(x) R 3 w(y -x)dΩ + ∂D ∂x i R 3 w(y -x)(y i -x i )dΩ (2.44) + 1 2! ∂ 2 D ∂x i ∂x j R 3 w(y -x)(y i -x i )(y j -x j )dΩ + 1 3! ∂ 3 D ∂x i ∂x j ∂x k R 3 w(y -x)(y i -x i )(y j -x j )(y k -x k )dΩ ...
Taking into account 2.40 and the fact that the odd integral terms are null and that the remaining integrals are constants (respectively c2, c4, . . . ), the nonlocal variable can be expressed as:

D(x) = D(x) + c 2 ∇ 2 D + c 4 ∇ 4 D . . . (2.45)
Thus, neglecting higher terms in expression 2.45 leads to the definition of the explicit nonlocal operator:

D(x) = D(x) + c 2 ∇ 2 D = D(x) + l 2 c ∆D (2.46)
The parameter c 2 is a squared length, so an internal length l c appears in the explicit formulation of the gradient as in the integral method. This explicit formulation is not adapted to the finite element method. Indeed, the computation of the Laplacian of the local variable D, which is defined on the integration points is problematic for two reasons: (i) It requires the extrapolation to the nodes which create a source of error. (ii) It equires computing the Laplacian of the extrapolated variable which is not necessarily C 1 . A solution to this limitation is to construct an implicit formulation [START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF]. To begin with, the Laplacian operator must be applied to 2.46:

∆D(x) = ∆D(x) -c 2 ∇ 4 D(x) (2.47)
Once substituted into 2.45, the following equation is obtained:

D(x) = D(x) + c 2 ∇ 2 D(x) -(c 4 -c 2 2 )∇ 4 D(x) . . . (2.48)
Finally, the implicit gradient formulation is than obtained after neglecting higher order terms:

D(x) = D(x) + c 2 ∇ 2 D(x) = D(x) + l 2 c ∆D(x) (2.49)
In the case of ductile failure, it was proposed in [Mediavilla et al., 2006b,a;[START_REF] Feld-Payet | Finite element analysis of damage in ductile structures using a nonlocal model combined with a three-field formulation[END_REF], to replace the local effective cumulative plastic strain p with its nonlocal counterpart p in the damage evolution variables. In particular, for the GTN model, the nonlocal cumulated plastic strain was used for the regularization of both the growth porosity (f g ) and nucleated porosity (f n ). The nonlocal plastic strain can also be used to regularize hardening as outlined in [START_REF] Peerlings | An implicit gradient plasticity-damage theory for predicting size effects in hardening and softening[END_REF]. However the implicit gradient formulation does not provide any guidelines to express the constitutive equations with respect to the nonlocal variable(s). Using p allows regularizing three local variables (p, f g and f n ) at the computational cost of one additional unknown field (p) [START_REF] Linse | Simulation of crack propagation using a gradient-enriched ductile damage model based on dilatational strain[END_REF]. However, in some cases the quantity trace(n) can be very localized and thus the growth porosity f g = (1 -f g )trace(n)p remains localized. A local tensor variable can also be used as shown is [START_REF] Hütter | Simulation of ductile crack initiation and propagation by means of a non-local Gurson-model[END_REF]. In this work, based on the plastic strain tensor ε p , a nonlocal counterpartε p is computed and then used in calculation of the evolution of the growth and nucleation porosity.

The idea of using multiple internal length scales was introduced in [Nguyen et al., 2020] for ductile failure. In the proposed model, an internal length was associated to each damage mechanism: (i) void growth based on the classical GTN model, (ii) coalescence by internal necking governed by a heuristic extension of the Thomason yield surface based on the maximum principal stress, (iii) shear-dominated coalescence mechanisms triggered by the maximum shear stress. Although the model proposes different internal length scales, the same characteristic length was assigned to the three nonlocal variables in the presented simulations.

Micromorphic models

Another approach, to regularize the ill-posed mechanical problem, is the micromorphic approach. This approach consist in enriching the kinematic description of the material behavior, by introducing an additional micromophic variable and its first gradient in the Helmholtz energy [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF]. The local internal variables have a micromorphic counterpart which takes into account the interactions at the microscopic level. In case of an elasto-plastic material, the effective cumulative plastic strain p and its micromorphic counterpart p χ are expressed in the Helmholtz energy [START_REF] Mazière | Strain gradient plasticity modeling and finite element simulation of lüders band formation and propagation[END_REF] as follows:

ψ = 1 2 (ε -ε p ) : E : (ε -ε p ) + ψ p (p) + 1 2 H χ (p -p χ ) 2 + 1 2 A χ ∇p χ .∇p χ (2.50)
Where A χ and H χ are two model parameters. E and ε e are respectively the elasticity and the elastic deformation tensors. ψ p is the contribution to the Helmholtz energy governing the hardening of the material. Thus, the thermodynamical associated forces are:

σ = ∂ψ ∂ε = E : (ε -ε e ) (2.51) R χ = ∂ψ ∂p = R + H χ (p -p χ ) (2.52) a χ = ∂ψ ∂p χ = -H χ (p -p χ ) (2.53) b χ = ∂ψ ∂∇p χ = A χ ∇p χ (2.54)
The yield surface is now expressed as:

φ = σ eq -(σ 0 + R χ ) = σ eq -(σ 0 + R + H χ (p -p χ ))
Unlike the gradient methods where no rule exists for the formulation of the constitutive law with the nonlocal variables, the micromorphic approach defines the coupling between the micromorphic variable p χ and the constitutive equation naturally in the writing of the thermodynamic force representing the hardening R χ .

The dissipation on the body Ω is expressed in the following equation:

D = Ω σ : εp -R χ ṗ -a χ ṗχ -b χ .∇ ṗχ dΩ (2.55)
Using the Green formula on eq 2.55, the dissipation potential can be rewritten as:

D = Ω (σ : εp -R χ ṗ -(A χ ∆p χ + H χ (p -p χ )) ṗχ ) dΩ (2.56) - ∂Ω A χ ∇p χ . n ṗχ dS
Based on the hypothesis that the dissipation potential related to the micromorphic variable p χ is null, the following equation is obtained:

A χ ∆p χ + H χ (p -pχ) = 0 ⇒ p χ = p + A χ H χ ∆p χ (2.57)
This equation shows an equivalence between the micromorphic approach and the gradient implicit approach where l 2 c = A χ /H χ . Thus the implementation of the micromorphic approach is very similar to the implicit gradient method.

Based on the general micromorphic approach presented in [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF], an extensions of the model to ductile fracture was introduced by [START_REF] Brepols | Gradient-extended two-surface damage-plasticity: Micromorphic formulation and numerical aspects[END_REF], a fully coupled damage-plasticity model where damage exhibits gradient effects. This proposed extension relies upon a 'two-surface' formulation in which plasticity and damage are treated with independent physical mechanisms, considering both a yield and damage function as well as appropriate loading/unloading conditions. The numerical implementation of the model in FEM was detailed and the results presented are independent of mesh size. Moreover, this model can well capture the initiation and propagation of cracks.

Another fully coupled damage model taking into account the strong coupling between all fields (Cauchy stress, kinematic hardening and isotropic hardening) and micromorphic damage was used in [START_REF] Diamantopoulou | Micromorphic constitutive equations with damage applied to metal forming[END_REF]. The implementation in the commercial FEM code ABAQUS/Explicit have been detailed. The implemented model has been successfully validated by simulation-test comparisons of a classical tensile test as well as by a bending test.

Another class of the micromorphic approach based on the improvement of Gurson's original homogenization procedure was initially introduced in [START_REF] Gologanu | Continuum Micromechanics, chapter Recent extensions of Gurson's model for porous ductile metals[END_REF]. The GLPD model is of micromorphic nature since it involves the second gradient of the macroscopic velocity and generalized macroscopic stresses of moment type, together with some characteristic microstructural length. The implementation of the GLPD method was detailed in [START_REF] Enakoutsa | Numerical implementation and assessment of the glpd micromorphic model of ductile rupture[END_REF], numerical 2D simulation were presented. The independence of the GLPD on mesh size was verified. In addition, a good agreement between the simulation and some experimental tests was reported. However, using the developed algorithm, a slow convergence rate of the model was obtained: the reported CPU time is 2 to 3 times longer than the standard local Gurson's model. In the work of [START_REF] Huetter | A micromechanical gradient extension of gurson's model of ductile damage within the theory of microdilatational media[END_REF], a gradient extension of Gurson's model which combines computational efficiency with a sound micromechanical basis was proposed. However no simulation was shown to support the numerical efficiency of the developed extension.

The micromorphic models use a local variable and its micromorphic counterpart. In the case where the two variables are needed to be very close so that to use the gradient of the local variable in practice, high values of the penalty term (Hχ in the case presented above) are employed. The use of high values of H χ results in poor convergence [START_REF] Scherer | Lagrange multiplier based vs micromorphic gradient-enhanced rate-(in)dependent crystal plasticity modelling and simulation[END_REF]. This problem is partially resolved with the next nonlocal approach which uses an energy enriched framework.

Enriched energy models

The enriched energy approach consists in regularizing the ill-posed mechanical problem by adding a positive function of the gradient of a local variable to the Helmholtz energy. Thus, minimizing the Helmholtz energy results in constraining the gradient of the local variable to have small values. Controlling the gradient of the local variable involves controlling the localization bandwidth by a parameter that can be related to the internal length of the material. According to [START_REF] Lorentz | A variational formulation for nonlocal damage models[END_REF][START_REF] Lorentz | Ill-posed boundary conditions encountered in 3d and plate finite element simulations[END_REF][START_REF] Lorentz | Gradient damage models: Toward full-scale computations[END_REF][START_REF] Zhang | Ductile damage modelling with locking-free regularised gtn model[END_REF][START_REF] Chen | Properties of a nonlocal gtn model within the context of small-scale yielding[END_REF], the Helmholtz energy is expressed on a Ω body as :

ψ = 1 2 (ε -ε p ) : E : (ε -ε p ) + ψ p (p) + 1 2 c∇p.∇p (2.58)
where c is a model model parameter (unit N) which can be expressed as c = σ 0 l 2 c so as to introduce a material characteristic length l c .

This approach has been applied in the case of a GTN model [START_REF] Zhang | Ductile damage modelling with locking-free regularised gtn model[END_REF][START_REF] Chen | Properties of a nonlocal gtn model within the context of small-scale yielding[END_REF]. In general, any localization variable of the GTN model can be chosen in the case of this approach. In order to regularize both the nucleation porosity f n and the growth porosity f g , the effective cumulative plastic strain p was chosen since it is involved in the evolution both porosities. Numerically, the gradient of the variable p, as p is defined at the integration point, cannot be calculated correctly as in the case of the explicit gradient model. Therefore, a relaxed formulation is proposed in [START_REF] Zhang | Ductile damage modelling with locking-free regularised gtn model[END_REF], introducing a nodal variable a and ensuring the equality of p and a in a weak sense. This equality is enforced by an ed Lagrangian using a Lagrange multiplier (l) used to enforce the weak equality between p and a ; it is expressed as follows:

L(ε, ε p , p, a, l) = Ω ψ(ε, ε p , p, a, l)dΩ (2.59) = Ω ψ loc (ε, ε p , p)dΩ + Ω ψ gdt (a)dΩ + Ω ψ rlx (p, a, l)dΩ
where

ψ loc (ε, ε p , p) = 1 2 (ε -ε p ) : E : (ε -ε p ) + ψ p (p) (2.60) ψ gdt (a) = 1 2 c∇a.∇a (2.61) ψ rlx (p, a, l) = l(a -p) (2.62)
The associated forces to the state variables are then expressed in the following:

σ = ∂ψ ∂ε = E : (ε -ε p ) (2.63) R nl = ∂ψ ∂p = R(p) -l (2.64) R a = ∂ψ ∂a = l (2.65) R l = ∂ψ ∂l = a -p (2.66) R ∇a = c∇a (2.67)
The yield surface is now expressed as:

φ = σ eq -(σ 0 + R nl ) = σ eq -(σ 0 + R -l)
Similar nodal forces are obtained as the ones in obtained with micromorphic models, except for the Lagrangian term (R l ). Note that a is introduced to numerically solve the problem. Therefore it should not be considered as the micromorphic counterpart of p. The equality between both variables is weakly enforced. As explained above the near equality between p and p χ is obtained for the micromorphic model using a very high value for H χ which then acts as a penalty factor. Both models then become close but CPU time becomes very large for the micromorphic model as H χ increases [START_REF] Scherer | Lagrange multiplier based vs micromorphic gradient-enhanced rate-(in)dependent crystal plasticity modelling and simulation[END_REF].

The dissipation retains its usual definition:

D = Ω (σ : ε)dΩ - L = Ω σ : εp -(R -l) ṗ -l ȧ -c∇a.∇ ȧ -(a -p) l dΩ (2.68)
Using Green's formula one obtains:

D = Ω σ : εp -(R -l) ṗ -(a -p) l -(l -c∆a) ȧ dΩ - ∂Ω c∇a. n ȧdS (2.69)
Considering that the new variables a and l do not contribute to the dissipation, their evolution is obtained in eq 2.70. Thus, as in the case of the micromorphic model, the strong form is obtained:

a = p l -c∆a = 0
(2.70) ∆a. n = 0 on ∂Ω Using this result, the yield surface can be rewritten as :

φ = σ eq -(σ 0 + R -l) = σ eq -(σ 0 + R -c∆a) (2.71)
which is very similar to the yield surface obtained using the micromorphic model. A locking-free FE formulation was also proposed in [START_REF] Zhang | Ductile damage modelling with locking-free regularised gtn model[END_REF] in a logarithmic total lagrangian large strains framework. The numerical implementation is detailed for the 3 then 5 field formulation. Simulation in 2D and 3D were presented. Mesh independence is verified with respect to mesh size and mesh arrangement. Furthermore, the model was applied to an experimental test base and good correlation with the experimental tests was obtained. However, it was reported that the convergence rate of simulations starts to decrease after the onset of crack propagation: cracked elements lose their stiffness and may become extremely distorted.

A solution to convergence problems was proposed in [START_REF] Chen | Properties of a nonlocal gtn model within the context of small-scale yielding[END_REF], by adding a penalty term which brings an additional coercivity so as to avoid the potential appearance of spurious plastic strain localization.

Phase field damage

The phase-field approach to model systems with sharp interfaces consists in incorporating a continuous field variable that differentiates several physical phases in a given system by a smooth transition. In the context fracture mechanics, such order parameter describes the smooth transition between the fully broken and intact material phases, thus approximating the sharp crack discontinuity, and is, therefore, referred to as the crack field. The evolution of this field as a result of the external loading conditions models the fracture process. A complete review of phase field based approaches is presented in [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF].

The phase field model for fracture was introduced in a variational form by [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]] as a generalization of Griffith's criterion to predict the critical stress for brittle fracture. The minimization of a functional that contains the sum of total elastic strain energy and fracture energy enables predicting the initiation, propagation, merging and branching of multiple cracks under complex loading conditions. The entire (quasi-static) process of crack initiation, propagation and branching is governed by a minimization problem of the energy functional:

E( u, Γ) = Ω ψ e (ε( u))dΩ + G c Γ dS (2.72)
where ψ e is the elastic energy density function, G c the material fracture toughness and Γ the admissible crack set.

In order to minimize this energy, it is needed to determine a priori the fractured surface in order to calculate the fracture energy. This results in a computational scheme that is non-tractable and inconvenient to apply. To overcome this limitation, another variational formulation based on regularizing the energy functional which is expressed, following the work of [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], as:

E l ( u, d) = Ω ((1 -d) 2 + η)ψ e (ε( u))dΩ + G c Ω 1 2l c d 2 + l c 2 l c |∇d| 2 dΩ (2.73)
where d is the phase field variable describing the degradation of the material; d = 0 represent intact material and d = 1 represent a totally broken material. l c is the characteristic length that controls the width of the transition zone, defined by d and η is a dimensionless parameter needed to prevent numerical difficulties that occur at total failure. Initially the phase field approach was developed for quasi-brittle material. However, recently this approach was extended to ductile fracture. For most damage models describing ductile fracture, the damage variable is related to the development of plastic zones. A phase field framework for ductile failure was proposed in [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF]. the Helmholtz energy of gradient damage approaches coupled with ideal elastic-plastic is given as:

ψ(ε, ε p , p, d) = (1 -d) 2 ψ e (ε e ) + (1 -d) 2 ψ p (p) + ψ d (d, ∇d) (2.74) = (1 -d) 2 1 2 (ε -ε p ) : E : (ε -ε p ) + ψ p (p) + G c 2l c d 2 + l 2 c |∇d| 2
The associated thermodynamic forces:

σ = ∂ψ ∂ε = ((1 -d) 2 + η)E : (ε e ) (2.75) R p = ∂ψ ∂p = ∂ψ p ∂p = (1 -d) 2 R(p) (2.76) R d = ∂ψ ∂d = -2(1 -d)(ψ e + ψ p ) + G c l c d -l 2 c ∇ 2 d (2.77)
The classical dissipation function:

D = Ω σ : εdΩ - d dt Ω ψdΩ (2.78) = Ω σ : εp -R p ṗ --2(1 -d)(ψ e + ψ p ) + G c l c d ḋ -G c l c ∇d∇ ḋ dΩ
Using the Green formula, the dissipation can be rewritten as:

D = Ω σ : εp -R p ṗ --2(1 -d)ψ e + G c l c (d -l 2 c ∆d) ḋ dΩ - ∂Ω G c l c ∇d. n ḋ dS
(2.79) If the dissipation related to the phase field variable is considered null, the strong form of the evolution of the phase field variable can be written as:

G c l c (d -l 2 c ∆d) = 2(1 -d)(ψ e + ψ p ) in Ω (2.80) ∇d. n = 0 on ∂Ω (2.81)
In the work of [START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF], the evolution of the regularized crack surface has been given in it local form as a function of the crack driving force H as:

η d ḋ = (1 -d)H -(d -l 2 c ∆d) in Ω (2.82) ∇d. n = 0 on ∂Ω (2.83)
By the identification of the two equation 2.80 and 2.82, the expression of the crack driving force is found:

H = 2 G c /l c (ψ e + ψ p ) (2.84)
In order to over come the numerical difficulties of the phase field models and to account for the irreversibly of the cracking process, the elastic energy ψ e was replaced by max s∈[0,t] ψ e (s) (see [START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]).

In brittle materials, the evolution of the phase field variable d only depends on the elastic energy. However, in case of ductile fracture, the plastic strain energy is added to the crack driving force H: a plastic degradation function is added in order to weaken the material where the plastic deformation is localized as proposed in [START_REF] Michael | A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF]. A similar approach was used in [Eldahshan et al., 2021b] coupled with remeshing to simulate ductile fracture. Its was shown that a fine remesh zone should be attributed to the localization band at the appropriate time in order to accurately predict the initiation of the crack in accordance with the experimental test.

A phase field approach adapted to porous ductile materials using a modified GTN-type plasticity model to account for a temperature-dependent growth of voids on micro-scale followed by crack initiation and propagation on macro-scale was presented in [START_REF] Dittmann | Phase-field modeling of porous-ductile fracture in non-linear thermo-elasto-plastic solids[END_REF]. Also simulations in accordance with experimental results in terms of hardening, necking, crack initiation and propagation were shown.

Conclusion

In order to overcome mesh dependency, one solution is to use nonlocal models. Several nonlocal models and the phase field approach were briefly recalled. The presented models regularize the ill-posed mechanical problem during softening. The principle of each method was detailed, and their extension to ductile failure was discussed. In this thesis work, only the implicit gradient method developed in chapter 3, is used to simulate ductile failure. The same model is chosen to simulate ductile failure using the crack insertion coupled with remeshing techniques as described in 4. Note that any other regularized model can be used in the developed framework for ductile fracture, using crack insertion coupled with remeshing.

Regularized damage models allows the modeling of crack initiation without the introduction of a pre-crack. However, continuum damage models still suffer from various problems during the crack propagation phase. The first problem is related to the fact that with a continuous damage field, the description of the crack lips (line in 2D simulations and surface in 3D simulations) is not possible. The second problem is related to the high deformation in the damaged zone. Especially near the crack front, all the softened elements stretch terribly, and thus possibly degrade the convergence rate of the calculation. These problems can be solved by using the so-called continuous-discontinuous methods. A brief review of the different continuous-discontinuous methods is presented in the next section.

Damage to crack transition

Continuous-discontinuous methods aim at modeling discontinuities/cracks based on a relatively localized continuous field. Several technical challenges must be addressed to achieve this goal. First, a discrete crack model must be defined to represent the discontinuity in the simulations. Second, a method should be chosen to determine the discrete crack localization from the continuous field, hereafter referred to as the discontinuity localization methods. It should be noted that these methods are not necessarily part of every damage to crack transition framework. And finally, a criterion should be formulated to indicate when the transition from a continuous damage field to its equivalent crack model is most relevant, hereafter called the insertion criterion.

Let us note that this review is mainly based on Feld-Payet [to be published]

Discontinuity modeling methods

In this section, the most common methods used in the literature to model discontinuity are presented. These methods are compared in terms of their accuracy to model the kinematics of crack propagation, their compatibility and ease of implementation within the finite element method and their applicability to ductile failure simulations.

Remove/kill elements

The kill elements method, also known as the remove element method, consists in eliminating the contribution of a considered broken element from the global rigidity matrix. In the following three variants of the same methodology are presented.

(i) Classical kill elements Only an insertion/deletion criterion (for example, an element is deleted if the chosen variable exceeds the critical value) is needed in the classical kill elements method. Each element which validates this criterion is removed from the simulation. This method is widely used [START_REF] Tvergaard | Influence of void nucleation on ductile shear fracture at a free surface[END_REF][START_REF] Tvergaard | Three dimensional microstructural effects on plane strain ductile crack growth[END_REF][START_REF] Gao | Cell model for nonlinear fracture analysis -ii. fractureprocess calibration and verification[END_REF][START_REF] Hambli | Finite element model fracture prediction during sheet-metal blanking processes[END_REF][START_REF] Saanouni | On the numerical prediction of the ductile fracture in metal forming[END_REF][START_REF] Xue | Numerical simulation of fracture mode transition in ductile plates[END_REF] because of its simplicity of implementation and its low computational cost. However, on the one hand, the use of this technique should lead to the suppression of a large band of elements in the case of nonlocal models. Lat us not that the bandwidth of the removed elements is not controlled. Subsequently, the computation body may suffer a mass loss, which may be quite significant depending on the size and/or the bandwidth of the removed elements. On the other hand, using the classical "kill elements" method results in a mesh-dependent crack path. In addition the crack lips are not smooth due to the removed elements (see figure 2.7).

(ii) Kill elements with remeshing This variant of the kill element method is obtained when coupled with remeshing techniques [START_REF] El | An anisotropic mesh adaptation strategy for damage and failure in ductile materials[END_REF][START_REF] Borouchaki | Adaptive remeshing in large plastic strain with damage[END_REF]. This variant aims at minimizing the size of the removed elements. The method thus allows to partially solve two of the three major problems of the "remove elements" Figure 2.7: Distribution of ductile damage inside the truncated side-pressed cylinder from [START_REF] Saanouni | On the numerical prediction of the ductile fracture in metal forming[END_REF]. The crack is represented in the simulation using the classical "remove elements" method. method, i.e. (i) minimizing the mass loss by minimizing the size of the removed elements, (ii) improving the crack lip surface with remeshing before element removal (see figure 2.8).

However, the problem of uncontrolled bandwidth of the removed elements can still lead to multiple rows of removed elements, and thus to unrealistic crack path.

Figure 2.8: Evolution of damage and mesh at the center of the specimen from [START_REF] El | An anisotropic mesh adaptation strategy for damage and failure in ductile materials[END_REF]. The crack is represented in the simulation using kill elements coupled with remeshing.

(iii) Controlled kill elements This variant of the "kill elements" method is obtained if the remove element technique is used with a discontinuity localization method (detailed in the next section) as proposed in [START_REF] Yun | A computational methodology for simulating quasi-brittle fracture problems[END_REF]. The discontinuity localization method determines the localization the crack surface from the continuous damage field (fig. 2.9a). Only the elements crossed by the localized crack surface are removed (fig. 2.9-b).

The improvement brought by this variant lies in controlling the bandwidth of the removed elements which is then limited to one row of elements.

Element cracking

A more realistic way to represent a crack is to consider for 3D problems a 2D crack surface (or in 2D, a 1D crack line). The geometry of the crack surface (or crack line in 2D) is then incorporated into the mesh. This representation of the localized crack can be done either by: (i) inter-element cracking, (ii) intra-element cracking.

(i) Inter-element cracking

The inter-element cracking representation consists in inserting a discontinuity at the interfaces between the elements. However in order to avoid mesh dependence, a pre-adapted mesh to the geometry of the crack is needed, which implies an a priori knowledge of the crack path. Practically, this method can be implemented by duplication nodes on the crack path in order to create the crack lips (see figure 2.10-a).

Figure 2.11: Damage distribution reported on the deformed configuration from [START_REF] Cuvilliez | A finite element approach coupling a continuous gradient damage model and a cohesive zone model within the framework of quasi-brittle failure[END_REF]. Intra-element cracking coupled with a cohesive band model inserted on a predefined planar crack path.

(ii) Intra-element cracking using remeshing Unlike the previous method, intraelement cracking with remeshing do not require prior knowledge of the crack path. However, an adaptive strategy is necessary to incorporate the geometry of the crack in the structure's mesh. A global or local mesh adaption is performed, in order to model the initiation of a new crack or the propagation of an existing crack. In practice, an intersection operation between the crack path and the mesh is necessary to generate the crack nodes in the mesh then duplicate them in order to create the crack lips (see figure 2.10-b) (see [START_REF] Chiaruttini | Advanced remeshing techniques for complex 3d crack propagation[END_REF]). Inter-element cracking has been used for ductile failure in 2D simulations [START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF]Mediavilla et al., 2006a] and extended to 3D in [START_REF] Feld-Payet | Amorc ¸age et propagation de fissures dans les milieux ductiles non locaux[END_REF][START_REF] Areias | Effective 2d and 3d crack propagation with local mesh refinement and the screened poisson equation[END_REF][START_REF] Javani | Three-dimensional finite element modeling of ductile crack initiation and propagation[END_REF]. In order to optimize the computational cost of the adaptive remeshing procedure, an error estimator [START_REF] Askes | Remeshing strategies for adaptive ale analysis of strain localisation[END_REF][START_REF] Rodriguez-Ferran | Error estimation and adaptivity for nonlocal damage models[END_REF][START_REF] Feld-Payet | Amorc ¸age et propagation de fissures dans les milieux ductiles non locaux[END_REF]] can be used: the error estimator generates an optimized mesh, refined in the areas where the damage field is localized and coarsened mesh in the areas where the damage field is more regular (see figure 2.12 ).

Let us note that both the previous representations of cracks can be coupled with cohesive zone bands. The cohesive band have its own interface behavior that can be coupled to the CDM model or not. An application of the inter-element cracking method coupled with cohesive band is presented in [START_REF] Cuvilliez | A finite element approach coupling a continuous gradient damage model and a cohesive zone model within the framework of quasi-brittle failure[END_REF].

eXtended Finite Element (X-FEM) method for crack representation

As opposed to the inter-element cracking method which necessitate mesh adaptation in order to describe the crack lips opening/closing, the X-FEM aims at representing the discontinuity using the same mesh.

The X-FEM method relies on the use of unit partitioning [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF]] to enrich a standard finite element approximation in displacement by discontinuous functions capable of modeling the crack opening (see figure 2.13).

Figure 2.13: Crack on a uniform mesh (left) and on a non-uniform mesh (right) from [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF]. The circled nodes are enriched by the jump function whereas the squared nodes are enriched.

The X-FEM method has been used with the Thick Level Set method to simulate crack propagation in softening materials. The TLS method was then used to regularize the mechanical problem at material softening and the X-FEM is used in the transition from continuous to discontinuous modeling (see [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF]). Also the X-FEM method has been used with a cohesive zone model (CZM) to simulate ductile fracture in [START_REF] Nikolakopoulos | Progressive failure of ductile metals: Description via a three-dimensional coupled czm-xfem based approach[END_REF] with a relatively coarse mesh. However the small strain framework considered in this paper may not be completely appropriate for ductile failure which is usually accompanied by large plastic strain.

Note: Although the remeshing technique for inter-element cracking and X-FEM offer the best representation of the crack lips kinematics during the simulation, they lead to a change of the discretization which brings its own challenges. These challenges are mainly related to the transfer of fields from the old mesh to the new mesh: transfer operations causes unavoidable errors, thus the mechanical stat after the transfer is unbalanced. This unbalance if not treated may results in the degradation of the convergence rate or the complete divergence of the resolution algorithm.

Discontinuity localization methods (DLM)

In the context of the local approach to fracture, discontinuity localization methods aim at defining a crack surface from a spatial (field scalar or tensor). In the case of a mesh with a preexisting crack, the discontinuity localization method must define the next crack increment.

However, in the case of a mesh without any crack, the discontinuity localization method must define where the first crack increment should be inserted. To the best of the author's knowledge, few procedures treat this case in the literature, especially for crack initiation away from any boundary.

Direction criteria-Pre-processing of DLM

Direction criteria aims at the orientation of the crack path. They are based on indicators that can either be a vector field to the crack path, or a scalar field that is maximal where the crack is supposed to be inserted. This indicator can be a global field defined on the whole mesh so the method is qualified as global, or a field only exploited in the vicinity of the crack path so the method qualified as local. This indicator is used as an input of the so-called crack path tracking algorithm to determine the exact position of a continuous crack path.

(i) Vector field from stability or bifurcation analysis This vector field can be computed in the case of coupled local (i,e. non-regularized) damage models. As mentioned in the section (1.5), localization is assumed to occur when it becomes possible to form a strain rate discontinuity in a planar band. This band is characterized by its unit normal n and the displacement jump across the band. The normal vector n can be obtained as the vector that minimizes the Rice criterion [START_REF] Rice | The localisation of plastic deformation[END_REF][START_REF] Rice | The mechanics of earthquake rupture[END_REF] :

det( n.D. n) = 0 (2.85)
where D in elasto-plastic materials is the elastoplastic tangent matrix ( σ = D : ε). In practice, the condition det( n.D. n) = 0 is never exactly met. So the normal n is obtained as the first occurrence of det( n.D. n) < 0. In [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] the consistent tangent matrix was used and gave similar prediction for the bifurcation analysis. However, for viscoplastic materials, linear perturbation analysis can be used to obtain the acoustic tensor n.D. n. Note that additional terms are added to 2.85 in case of finite strain ( see [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] for hypo-elastic material or [START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF] for hyper-elastic material). This approach was used to model ductile failure in [START_REF] Crété | Numerical modelling of crack propagation in ductile materials combining the GTN model and X-FEM[END_REF][START_REF] Wolf | Numerical modeling of strain localization in engineering ductile materials combining cohesive models and xfem[END_REF][START_REF] Huespe | A finite thickness band method for ductile fracture analysis[END_REF] using GTN model within the small strain framework and extended to large strain framework in [START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF].

This type of analysis has very strong arguments for the choice vector field indicator of the crack path, but requires a very high computational cost for 3D simulations.

(ii) Vector field from principal component analysis (PCA) Another technique is based on PCA of a tensor field, in order to determine normals to the crack path. This technique is computationally less expensive than the stability or bifurcation analysis and can be used in case of nonlocal models. Different tensor fields have been used in literature:

1. The strain tensor was used in 2D [START_REF] Dufour | Estimation of crack opening from a two-dimensional continuum-based finite element computation[END_REF] and in 3D [START_REF] Areias | Effective 2d and 3d crack propagation with local mesh refinement and the screened poisson equation[END_REF] for quasi-brittle fracture simulation;

2. The stress tensor used in 2D for quasi-brittle fracture simulation [START_REF] Bobinski | in cons[END_REF];

3. The anisotopic damage tensor was used in for 3D quasi-brittle fracture simulation in [START_REF] Javanmardi | Extended finite element method and anisotropic damage plasticity for modelling crack propagation in concrete[END_REF].

The vector field corresponds to one of the components of the eigenvectors obtained from the principal component analysis of the considered tensor.

The bifurcation analysis and the PCA of a tensor provide several vectors for each evaluation point. But only one vector is needed to define the normal to the crack front (in the absence of branching). To choose the appropriate vector field, in the literature, most authors using these methods follow this recommendation: the direction of propagation is the direction of the output vector field that maximizes the displacement gradient [START_REF] Rabczuk | A three-dimensional large deformation meshfree method for arbitrary evolving cracks[END_REF] or the equivalent deformation [START_REF] Crété | Numerical modelling of crack propagation in ductile materials combining the GTN model and X-FEM[END_REF].

(iii) Scalar fields Various scalar fields have been used in the literature to determine the crack surface. Some authors use directly the damage field (see [START_REF] Broumand | The extended finite element method for large deformation ductile fracture problems with a non-local damage-plasticity model[END_REF][START_REF] Seabra | Damage driven crack initiation and propagation in ductile metals using xfem[END_REF][START_REF] Feld-Payet | Amorc ¸age et propagation de fissures dans les milieux ductiles non locaux[END_REF][START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF][START_REF] Bottoni | Topological search of the crack pattern from a continuum mechanical computation[END_REF]). But other scalar fields can be considered, such as the cumulative plastic strain which is used in the calculation of the damage field in the case of the GTN model.

In the case of quasi-brittle failure, the scalar field can be the L 2 norm of the eigenvectors V ε i of the deformation tensor [START_REF] Simone | From continuous to discontinuous failure in a gradient-enhanced continuum damage model[END_REF]:

ε eq = Σ i < V ε i > 2 .
Another example of a scalar field is the equivalent strain calculated on the basis of the maximum eigenvalue maxσ i of the stress tensor [START_REF] Bobinski | in cons[END_REF]: ε eq = maxσ i E .

From a direction criterion to the crack surface-crack path tracking algorithm

The crack path tracking algorithms use the vector field or scalar field mentioned in the previous section to define a continuous crack path. This definition can be either explicit, through the construction of an auxiliary crack mesh, or based on the iso-value (equal to θ c ) of a scalar field θ evaluated at the nodes of the structure's mesh.

(i) Global method This method was initially presented in [START_REF] Oliver | On strategies for tracking strong discontinuities in computational failure mechanics[END_REF] for single 2D and 3D crack propagation and in [START_REF] Oliver | Continuum approach to the numerical simulation of material failure in concrete[END_REF] for multiple 3D crack tracking with examples of quasi-brittle fracture simulations. This method considers as input data the vector field of the normals to the crack path n, to compute the tangential vector(s) t in 2D (or ( t, s) in 3D), of a scalar field θ(x) (see figure 2.14.a), so that:

∂θ ∂t = 0 for 2D cases ; ∂θ ∂t = ∂θ ∂s = 0 for 3D cases (2.86)
In order to obtain the scalar field θ (see figure 2.14.b), a heat-conduction-like problem with adiabatic heat flux boundary condition at the boundary ∂Ω and no internal heat source need to be solved:

div( q) = 0 in Ω (2.87) q = -( t ⊗ t + s ⊗ s).∇θ in Ω (2.88) q. v = 0 on ∂Ω (2.89)
As described in figure 2.14.c, from the nodal values of the monotonous field θ the crack path Γ c can then be implicitly defined:

Γ c = {x ∈ Ω/θ(x) = θ c } (2.90)
Where θ c is a threshold value set by an adiabatic heat flux boundary condition.

Figure 2.14: Scheme of the global method algorithm from [START_REF] Oliver | Continuum approach to the numerical simulation of material failure in concrete[END_REF] This problem solved thanks to the finite element method on the whole mesh, hence the name "global method". Note that the stiffness matrix K = (t ⊗ t + s ⊗ s) is singular. This is why a perturbation term is added in the form

K ε = K + ε1.
This global method was used in [START_REF] Cervera | Mesh objective tensile cracking via a local continuum damage model and a crack tracking technique[END_REF] in the case of a quasibrittle problem with strain localization in order to track cracks and to allow the crossed elements to crack in order to obtain an objective mesh solution. Also In the work of [START_REF] Dufour | Estimation of crack opening from a two-dimensional continuum-based finite element computation[END_REF], the global method was used to estimate the crack opening in quasibrittle simulation. In the case of ductile failure, the same method was used within the small strain framework [START_REF] Huespe | A finite thickness band method for ductile fracture analysis[END_REF] to track cracks and represent them using a "strong discontinuity" method, and then adapted to large deformations in [START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF]. More recently, the same method has been used for polycrystalline materials to select the elements crossed by the crack path and then enrich the surrounding nodes using X-FEM method [START_REF] Beese | 3d ductile crack propagation within a polycrystalline microstructure using xfem[END_REF].

The global method uses as input the normals to the crack paths n, a global quantity. This global quantity varies little along time in the case of quasi-fragile simulations as indicated by [START_REF] Oliver | Continuum approach to the numerical simulation of material failure in concrete[END_REF]. This means that the crack path can be almost entirely determined at the early stages of the degradation. However, in other cases where vector field n, and thus the crack path, may vary in time, it is necessary to update the n field and to consider crack propagation along the current θ iso-value that crosses the last element on the consolidated part of the crack path.

(ii) Crack-path field technique This technique was originally introduced in [START_REF] Dias | Strain injection techniques in numerical modeling of propagating material failure[END_REF] as an improvement of the global method. it requires more inputs than the global method: Indeed, in addition to the vector field of the normals to the crack path n, it uses a localization scalar field, such as the equivalent deformation field ε eq . The mathematical framework of the method has been introduced in [START_REF] Oliver | Crack-path field and strain-injection techniques in computational modeling of propagating material failure[END_REF]: the idea is to use the derivative of the scalar field ε eq in the direction n to define the crack path thanks to its zero iso-value. To do so, the first step is to consider the localized field ε eq (a field computed at the integration point) in order to obtain a smoother nodal field ψ t = ψ(t) equivalent to the local form in a weak sense by solving:

∀ψ B ψ (x) (ψ t -ε eq (x, t)d) dB = 0 (2.91)
where ψ are test functions. It is then possible to consider the derivative of this smoothed field along the assumed normal vector n given as an input data:

∂ψ t ∂ n (x) = ∇ψ t . n.
In order to determine the zero iso-value of this new scalar field, it is possible to build a smoothed version of it using the same procedure:

∀µ B µ (x) µ t - ∂ψ t ∂ n (x)d dB = 0 (2.92)
Figure 2.15: Schematic distribution of the function µ in a 2D problem, where the crack path Γ is defined as the zero level set of µ [START_REF] Dias | Strain-injection and crack-path field techniques for 3d crack-propagation modelling in quasi-brittle materials[END_REF] As described in 2.15, the crack path Γ c is then be retrieved in an element basis from the smoothed µ:

Γ c = {x ∈ B/µ t (x) = 0} (2.93)
Note that this method is qualified as a local crack path tracking algorithm because of the use of smoothing operations only on a restriction B of the domain Ω where the scalar variable is strongly localized.

The crack-path field was used for crack tracking to simulation 2D crack propagation in [START_REF] Oliver | Crack-path field and strain-injection techniques in computational modeling of propagating material failure[END_REF]] [Dias et al., 2016] and extended to 3D simulations in [START_REF] Dias | Modeling tensile crack propagation in concrete gravity dams via crack-path-field and strain injection techniques[END_REF] for quasi-brittle cases. This method is more efficient than the global method in terms of computational cost: the local smoothing operation are done locally and no global resolution of an auxiliary equation is needed. However, the computational cost associated with the computation of the normals to the crack path using a bifurcation analysis or a PCA analysis remains significant, especially is cases where the normals to the crack paths vary during the propagation process which is the case for ductile fracture.

A major limitation of this method is related to the choice of a vector field assumed to be sufficiently close to the normals to the crack path. This choice is far from obvious, especially in cases where a bifurcation analysis is not relevant. Approximations to the vector field n using the gradient of a function f of displacement field u(x, t): ∇f (u(x,t)) ||u|| ≈ n were proposed in [START_REF] Lloberas-Valls | Strain injection techniques in dynamic fracture modeling[END_REF] in case of dynamic simulation. To the author knowledge, this crack-path field technique was not used in case of nonlocal models to simulate ductile fracture.

(iii) Marching ridges algorithm This method Feld-Payet et al. [2015] shares similarities with the crack-path method field. For instance they both aim to construct an scalar function with a zero iso-value corresponding to the crack path. However, the marching ridges algorithm, by considering a finite set of vectors to test whatever the application may be instead of guessing a vector field sufficiently close to the crack normals, simplifies the inputs: its only requires the gradient vector of a localized scalar field. By the same token, this algorithm is more efficient cost wise then the crack-path method since the heavy computational step of normals to the crack is not needed.

The first step of this algorithm is a smoothing operation on the considered scalar variable f (computed at integration points). This smoothing operation can be done following the super convergent patch recovery method described in [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF]. From the smoothed field f defined at the nodes, the computation of the spacial gradient ∇f is straightforward. The idea is then to approximate locally the crack path by a straight kine over a distance R and to determine its direction defined by an angle. To do so, a polar grid centered at the last known point of the crack path is used (see figure 2.16). From here a second scalar function µ(x) = ∇f. e θ is evaluated, where e θ is the tangential vector to the polar grid at a given angle θ. The crack corresponds to the zero iso-level of the function µ:

Γ c = {x ∈ Ω/µ(x) = 0} (2.94)
In practice, the µ function is evaluated only at a finite number of discrete points, regularly distributed on a polar grid, and the zero iso-level is considered to be locally well approximated by the straight line corresponding to the angle between two evaluations of µ where its sign changes from positive to negative.

Another difference of this method is the explicit definition of the crack path after the polar grid search unlike the implicite definition with the above mentioned methods. One The adequate choice of the radius of search is important in order of obtain an accurate crack path. In particular, this algorithm is only valid if the local radius of curvature is smaller than the distance between the point of evaluation and the considered point of the crack front R (see figure 2.17). The evaluation radius should enlarged for treating quasi-plateau areas with small evaluation radii.

(iv) Local maximum search over different radii This crack path tracking algorithm was originally introduced in [START_REF] Brokken | Numerical modelling of the metal blanking process[END_REF][START_REF] Brokken | Numerical modelling of ductile fracture in blanking[END_REF] in the case of uncoupled damage models in 2D simulations. It considered a scalar field S which is evaluated on polar grids in order to determine the direction Θ i dir corresponding to the maximum of S at radius r i . Finally the selected direction Θ dir is the median of the directions Θ i dir . The input scalar field considered by [START_REF] Brokken | Numerical modelling of ductile fracture in blanking[END_REF] was the S = ∆R/R Rc , where ∆R/R is the void growth according to the Rice & Tracey criterion and R c is the critical void growth. More precisely, the authors used the nodal extrapolation S n of this field S, originally computed at the integration points (with 0 ≤ S n ≤ 1) which is evaluated at several radii r i = 1/2∆L, 3/4∆L, 5/4∆L, 3/2∆L (see figure 2.18). Figure 2.18: Schematic representation for of the multiple polar grid search for the maximum from [START_REF] Brokken | Numerical modelling of the metal blanking process[END_REF] This method was used in the case of nonlocal coupled damage models in [START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF]Mediavilla et al., 2006a,a] in 2D simulations and then extended to 3D simulation in [START_REF] Javani | Consistent remeshing and transfer for a three dimensional enriched mixed formulation of plasticity and non-local damage[END_REF][START_REF] Javani | Three-dimensional finite element modeling of ductile crack initiation and propagation[END_REF]. The input scalar variable used is the nonlocal damage variable ω p .In the case of nonlocal models, if the damage variable is used as the nonlocal variable, this crack path tracking algorithm is the most efficient with respect to computational cost for the following reasons: (i) the input nodal scalar field is obtained directly from the nonlocal computation, (ii) the polar grid search is performed locally around the crack tip (in 2D) or the crack front (in 3D), (iii) the crack path is also expressed explicitly by the x i c positions. However, since the local maximum is sought by the polar grid method, only one Θ dir direction is obtained. Therefore, in the presence of branches, only the branch where the damage is maximum will be found. Moreover, the resulting crack path depends on the correct choice of the evaluation radii ∆r i as it is the case for the Marching ridge algorithm presented earlier.

(v) Maximum of projection onto a geometrical form This approach is quite similar to the local maximum search over different radii. The principle is to search for the maximum of the projection of the considered scalar field onto a line. As with any method dealing witg a scalar field, there is a smoothing operation involved before the final evaluation. Smoothing can be performed, for example, by using Bézier curves [START_REF] Seabra | Damage driven crack initiation and propagation in ductile metals using xfem[END_REF], or a convolution product between the scalar field and a Gaussian function [START_REF] Bottoni | Topological search of the crack pattern from a continuum mechanical computation[END_REF]. This smoothing operation is done to obtain a field rather independent of the mesh. As far as the projection is concerned, it can be done on a straight line [START_REF] Bottoni | Topological search of the crack pattern from a continuum mechanical computation[END_REF] or a portion of circle [START_REF] Seabra | Damage driven crack initiation and propagation in ductile metals using xfem[END_REF]. The position of the geometric shape is determined as the projection of current crack front considering the vector ---→ x i-1 x i at a given distance ∆L (see figure 2.20). From the smoothed field projected on the geometrical shape, the maximum of the latter is the next point on crack-path. Then, the crack path is determined explicitly in a step-by-step manner (see figure 2.19).

Figure 2.19: Scheme of the search by projection on the orthogonal line to the crack path procedure from [START_REF] Brokken | Numerical modelling of the metal blanking process[END_REF] Figure 2.20: Scheme of the search by projection on a circle around the the crack path procedure from [START_REF] Seabra | Damage driven crack initiation and propagation in ductile metals using xfem[END_REF] This approach was presented and validated only in 2D simulations. The proposed approach benefits from the same computational cost efficiency as the previous method. But for the same reason, no branching detection is possible.

(vi) A medial-axis-based approach This approach was introduced in [Tamayo- [START_REF] Tamayo-Mas | A medial-axis-based model for propagating cracks in a regularised bulk[END_REF] to track cracks in a regularized medium. The algorithm takes as input a smooth scalar damage field, and considers the path of the crack to be the medial-axis (or the medial surface in 3D) of the localization band. The medial-axis is constructed as the centers of circles/spheres that are bi-tangential to the a given iso-level surface of damage (see figure 2.21). This approach can track cracks in a 2D or 3D simulation and can also capture branching as shown in [START_REF] Tamayo-Mas | A medial-axis-based model for propagating cracks in a regularised bulk[END_REF][START_REF] Tamayo-Mas | A continuous-discontinuous model for crack branching[END_REF]. However, this approach is not able to deal with crack initiation away from any boundary. Besides, it suppose that the regularized damage field is symmetric with respect to the crack path. Which is not necessarily the case if there is some change during loading in the propagation direction.

Initiation of cracks methods

In the previous section dealing with direction criteria the methods were explained assuming the existence of a pre-crack or if the structure is crack free, assuming the knowledge of at least a point on the future crack surface.

The problem of crack initiation is often simplified in the literature. Most authors consider only the simulation of crack initiation at the structure's surface [START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF][START_REF] Seabra | Damage driven crack initiation and propagation in ductile metals using xfem[END_REF]. However, only few authors have addressed crack initiation inside structure. The problem is even more challenging in 3D where the first crack increment is a surface whereas it is simply a point in 2D.

In this section, the different crack initiation methods are briefly presented keeping in mind these difficulties.

(i) Initiation of 2D crack using maximum search or the marching ridges algorithm When damage is not coupled with the constitutive behavior, it is possible to use material bifurcation as a criterion to define in which element a crack should initiate. Another option is to search for local maxima on a scalar field to define a starting point for the crack path as in [START_REF] Brokken | Numerical modelling of ductile fracture in blanking[END_REF][START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF][START_REF] Feld-Payet | Amorc ¸age et propagation de fissures dans les milieux ductiles non locaux[END_REF][START_REF] Bottoni | Topological search of the crack pattern from a continuum mechanical computation[END_REF][START_REF] Seabra | Damage driven crack initiation and propagation in ductile metals using xfem[END_REF]. Let us note that depending on the smoothing technique used to construct the input scalar field, this starting point can either be a node in the mesh or any point in the structure.

Let us note that it is usual to take some measure to avoid cracks initiation in the wake of already existing cracks. A simple solution allow initation of only a single crack in each separate localization band. To distinguish different localization bands, a threshold should be defined and the smoothed input scalar field can be constructed only locally for each cluster of points exceeding this specific threshold.

(ii) Initiation in3D with implicit definition of the crack When considering a nodal scalar field where a specific value of this field indicates the presence of the crack (2.90 2.93), the crack surface corresponding to crack initiation can be determined in quite a similar way as for crack propagation. Indeed a scalar field is always evaluated at the element level to find points on the edges of the element that coincide with the crack path. A crack is initiated if three different edges of the element are intersected by the surface representing the crack path, thus defining the crack surface of the first crack ( [START_REF] Oliver | On strategies for tracking strong discontinuities in computational failure mechanics[END_REF][START_REF] Oliver | Crack-path field and strain-injection techniques in computational modeling of propagating material failure[END_REF]).

(iii) Initiation of 3D crack using a centered rectangle in a cluster of damaged points This method was described in [START_REF] Javani | Three-dimensional finite element modeling of ductile crack initiation and propagation[END_REF] to initiate cracks in the outer surface of the mesh or in the body of the mesh. The initiated 3D crack is modeled by a rectangle placed in a group of nodes that have exceeded the critical damage value f crit . The center of the rectangle is computed as a weighted average of the damaged points {x i crit ∈ Ω / f (x i c ) > f crit } (where f is the scalar damage field). In the case of a near-surface crack, the weights are calculated as the inverse of the distance to the external surface multiplied by the damage values, so that the center will be the closest point to the external surface. The plane of the rectangle is the median plane through the cluster of damaged points. The dimensions of the rectangle are calculated to obtain the minimum rectangle containing all orthogonal projections of the damaged points.

Figure 2.22: 3D crack initiation algorithm from [START_REF] Javani | Three-dimensional finite element modeling of ductile crack initiation and propagation[END_REF] Let us note that the number of damaged points must be carefully chosen to obtain a rectangle with the appropriate normal.

(iv) Initiation of 3D crack using Marching ridges algorithm Another possible approach is to use the same initiation and propagation algorithm as proposed in [START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF]. In 3D, the Marching ridge algorithm starts the search from a discretized line. For each segment of the line, its next projection is computed using the 2D Marching ridge algorithm in the plane orthogonal to the segment. Thus, the Marching ridge algorithm in 3D only needs one segment to initiate the crack. Figure 2.23: 3D crack initiation scheme using the marching ridges in 3 orthogonal planes from [START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF] In order to determine the first segment for crack initiation, the global maximum of the smoothed scalar field is first determined. Then, three orthogonal planes centered on the corresponding point denoted x 0 and representing the global reference frame are defined(see 2.23). A search with the 2D Marching ridges algorithm is performed in each of these planes, to provide a maximum of the three intersection lines with the damage crest surface. The longest intersection is then selected and a new 2D search in an orthogonal plane can be performed with the Marching ridges algorithm. This search enables to determine the second direction to build a rectangular shaped crack increment.

Using this approach, a unified method for initiation and propagation is proposed. However, the computational cost of the three-plane search is not negligible. Especially if it is used during a computation where the algorithm is called at the end of each converged load increment.

Insertion criteria

The purpose of an insertion criterion is to determine the right time for the transition from continuous to discontinuous modeling. The insertion criterion depends on the nature of the constitutive law. In particular, in the case where a coupled local damage model is used, the insertion criterion must select the insertion time in order to avoid the loss of solution uniqueness that occurs with material softening. However, in the case of a nonlocal damage model, the uniqueness of the solution is guaranteed by the nonlocal treatment. The insertion criterion must then allow the crack to be inserted, unless a cohesive zone model is used.

Furthermore, the evaluation of the insertion criterion also depends on the chosen representation of the crack. In the case where the representation is mesh dependent, the evaluation of the insertion criterion can be performed at the integration points for the kill element method, or at the nodes of the mesh for inter-element cracking method. However, in the case of a mesh-independent representation of the calculation (e.g. intra-elemental remeshing cracking or X-FEM), the input field of the criterion should be related to points on the crack surface.

(i) Stability and bifurcation criteria In case of coupled local damage model, this criterion is usually used with the global method or the crack-path field method. This criterion select the moment of insertion of the crack in an element Ω e is usually allowed if these two conditions are met:

1. The bifurcation criterion 2.85 is validated within Ω e .

The crack surface from the discontinuity localization method crosses the considered

element Ω e .

Representation of the crack is updated whenever these conditions are me in a new element. This is why, this method is usually used with enrichment method for crack.

(ii) Critical value criterion When considering a spacial field and a critical value, the crack increment is inserted when the chosen field reaches the critical value in the corresponding area.

A quantity as the equivalent strain can be used as an input of the insertion criterion [START_REF] Jirásek | Embedded crack model. part ii: combination with smeared cracks[END_REF]. This choice is interesting because the critical equivalent strain could eventually be determined experimentally using image correlation if crack initiation and propagation is properly detected.

In case where a damage variable is available in the material model, it can be chosen as the input field of the insertion criterion.

(i) For a normalized damage variable in the sense of Kachanov, D = 1 represents a completely damaged material and D = 0 represents intact material.

1. If the discrete crack is inserted directly without a cohesive band, the natural critical value should be f crit = 1. Practically for numerical reasons, this value is approximated f crit ∈ {0.9, 0.99, 0.999}, and the higher critical value, the later the insertion is authorized by th insertion criterion as reported in [START_REF] Simone | From continuous to discontinuous failure in a gradient-enhanced continuum damage model[END_REF].

2. If a cohesive band is used in order to dissipate the remaining energy, the critical value is chosen depending on the identification parameters for the cohesive model. Therefore lower values of f crit can be used.

(ii) In case of micro-mechanical based models as the GTN, the damage variable is the porosity in the material. For these models the critical porosity is determined as the porosity corresponding to the zero effective stress. For the GTN model:

f crit ≈ 1/q 1 .
(iii) Dimensions of the crack increment For crack representation methods that proceed on an element basis, the dimension of the crack element is not controlled by the used, but directly by the verification of the insertion criterion. However, when inserting crack increments larger than one element, the user can choose the minimal size of the crack increment to be inserted. The dimensions of the crack increment should be carefully selected for the following reasons: if a large crack increment is chosen, the crack insertion will be delayed, and therefore the softened elements in the crack tips will be overstretched which may affect the convergence rate. In the opposite case, if a too small crack increment is inserted, the remeshing will be triggered too often and the calculation cost will be higher. This is why, a minimal crack-increment has been imposed in [START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF]; [START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF]; [START_REF] Javani | Three-dimensional finite element modeling of ductile crack initiation and propagation[END_REF].

Conclusion

One of the main objectives of this thesis is to allow the use of a nonlocal model to guide the initiation and the propagation automatically over long distances. In this case, a fine mesh along all possible crack paths is required, and thus the propagation distances are limited by the size of the mesh i.e. limited to the number of degrees of freedom. The selected method based on a continuous-discontinuous transition must be adapted to this problem.

First, the crack representation adopted for the rest of this thesis is inter-element cracking with mesh adaption. This choice is based on three considerations: (i) inter-element cracking provides the best representation of the crack lips while being mesh independent, (ii) mesh adaption avoids elements distortion due to large plastic deformation (especially ahead of the crack front), (iii) while ensuring a good quality/cost ratio by enabling both to refine the elements in the process zone before crack insertion to properly capture the nonlocal damage evolution and then to increase the element size after crack opening.

Secondly, for the discontinuity localization method, the Marching ridges algorithm, is selected to allow a robust detection of cracks, potentially initiating crack completly inside the structure in 2D and 3D. Finally, as the chosen continuous damage model is a nonlocal GTN model, the natural choice for the insertion criterion is the critical porosity criterion. With this choice, crack insertion is necessarily performed within zones that transmit almost no forces.

A remaining challenge for this thesis is to define the optimal dimensions of the crack increment, in particular in 3D. This choice must take into account both the remeshing frequency and the convergence rate after each remeshing.

In the following part, the different elements required for automatic mesh adaption, even without crack insertion are tackled.

4 Mesh adaption for an efficient continuous-discontinuous approach

Introduction

While nonlocal models enable to solve the mesh dependence problem for modeling ductile cracks, they however require a fine mesh along the localization band, at least on the areas where degradation is evolving. For crack propagating over long distances, using a fine mesh throughout the structure is not possible, otherwise the computational size may become too large.

In order to keep a reasonable computational cost, this work proposes to resort to an automatic mesh procedure to provide an optimal mesh characterized by: (i) a fine mesh in the active localization band (i.e. near the area where to the future crack increment should be inserted), (ii) a coarser mesh in the areas far from the crack path or in the already cracked areas. In order to maintain a good compromise between computational cost and accuracy, a possible solution is to couple the remeshing procedure with a discretization error estimator.

In the first part of this section, the most commonly used discretization error estimators are recalled, then a remeshing procedure based on these error estimators is presented. Finally, the transfer operators that are necessary whenever there is a change in discretization are also presented.

Discretization error estimators

A discretization error estimator informs on the capability of a mesh to represent the unknown solution.

This can be done prior to the computation so the error estimator is called "a priori error estimator" of after to the computation so the error estimator is called "a posteriori error estimator" [START_REF] Tinsley Oden Ainsworth | A posteriori error estimation in finite element analysis[END_REF].

Only a posteriori error estimator will be presented, because the priori error estimators cannot be implemented in FEM software, in most cases (theses estimators use incalculable constants to bound the exact error ).

A posteriori error estimators can be regrouped in three classes:

-Error estimators based on the equilibrium defaults (see ??) in the finite element solution.

-Error estimators by measuring the error between a constructed admissible field and the field computed using the constitutive equation, -Error estimators built on the irregularity defects of the finite element solution, At first, in order to introduce briefly discretization error estimator for each category, a simple static mechanical problem on an elastic linear solid will be introduced. The adaptation to non-linear problem will then be presented for the chosen discretization error estimator.

(i) Discretization error The discretization error resulting from a finite element approximated displacement solution, denoted u h , to the standard mechanical problem is defined as: (see [START_REF] Tinsley Oden Ainsworth | A posteriori error estimation in finite element analysis[END_REF])

e = u -u h (2.95)
where u stands for the exact solution. The same definition can be written in terms of the approximated stress σ h :

e σ = σ -σ h (2.96)
The discretization error is generally measured thanks to the energetic norm, which has a strong physical meaning and is equivalent to other norms like the L 2 norm:

η ex = || e|| H,Ω = [ Ω ε( e) : H : ε( e)dΩ] 1 2 η ex,σ = ||e σ || H -1 ,Ω = ||σ -σ h || H -1 ,Ω = [ Ω (σ -σ h ) : (H -1 : (σ -σ h ))dΩ] 1 2 (2.97)
where η ex and η ex,σ are respectively, the exact global measure of the error on the displacement and the stress field. H is the fourth order elasticity tensor and Ω ∈ R n is the domain of the solid structure A local form of the exact measure of error can be expressed within each element of the mesh as shown bellow. This local error is necessary to all adaptive procedures.

Within each element E K of the mesh, an elementary contribution to the global error η ex can be defined as:

η ex E k = || e|| H,E k = E k ε( e) : H : ε( e)dΩ 1 2 avec E k ∈Ω h || e|| 2 H,E k = || e|| 2 H,Ω h (2.98)
(ii) Quality indicator of discretization error estimators Since the exact solution is rarely known, the same goes for the exact error. The exact error η ex is only evaluated in some benchmark test cases, where the analytical solution is known, to quantify the efficiency of the estimation. However, η ex in most cases is out of reach. In practice, a solution on a more finer mesh u f in h can be used to approach η ex . But the aim of the different error estimators is to avoid such an expensive computation while providing the best approximation possible of the exact error. In order to compare different discretization error estimators, some indicator have been proposed in the literature. The most popular and straightforward indicator is the efficiency index, which is defined as:

γ = η η ex (2.99)
where η is an estimation of the exact error.

The closer the efficiency index to 1, the more precise the estimator is (i.e γ ≈ 1 is equivalent to η ≈ ηex). If the efficiency index γ > 1, the estimator is said to be reliable because it's provides an upper bound to the exact error.

Other properties of an error estimator are taken into account to evaluate its performance as, the cost in terms of computation, the complexity/intrusion-level in terms of its implementation in FEM software, etc.

Residual discretization error estimators

This category of error estimators introduced initially by [I. Babuška, 1978], exploits the lack of equilibrium in the FEM approached solution σ h to evaluate the error due the discretization. This type of error estimator exploits the fact that e = u -u h is a solution to the mechanical problem. If this quantity is injected in the weak form of the mechanical problem, the following equation is obtained:

Ω σ( e) : εdΩ = E∈Ω h Ω E r E v dE + Γ∈∂Ω E Ω E r Γ v dΓ (2.100)
where r E is the residual within each element E of the mesh under the volume forces F , and can be expressed by:

r E = div(σ h|E ) + F |E (2.101)
and r Γ is the residual on at the interface of two adjacent elements E and E , and can be expressed by:

r Γ =      σ h|E • n E + σ h|E • n E if Γ = ∂E ∩ ∂E σ h|E • n -F d if Γ ∈ ∂ f Ω 0 if Γ ∈ ∂ u Ω (2.102)
An upper bound of the discretization error can then be estimated as(see [START_REF] Bank | Some A Posteriori Error Estimators for Eleptic Partial differential Equation[END_REF]):

|| e|| H,Ω ≤ η = C   E∈Ω h 2 E ||r E || 2 L 2 (E) + Γ∈∂Ω h Γ ||r Γ || L 2 (Γ)   1 2 (2.103)
where h E and h Γ are respectively the minimal element size of volume elements and boundary elements.

All the terms of the second member of 2.103 are explicitly computable but the constant C. This constant is difficult to evaluate and in most cases mesh dependent. From this global estimator, local contribution of error η E can be derived for the adaptive procedure.

The estimator defined in 2.103 is called the explicit residual error estimator. Implicit versions also exist, such, as the flux free residual estimator [START_REF] Díez | Accurate upper and lower error bounds by solving flux-free local problems in "stars[END_REF] and the equilibrated residual estimator [START_REF] Ainsworth | A modified gurson-type plasticity model at finite strains: formulation, numerical analysis and phase-field coupling[END_REF]. They use the resolution of the problem 2.100 locally on elements or patches of elements.

Error estimators using the constitutive equation

This category of error estimator was introduced initially by [START_REF] Ladeveze | Error Estimate Procedure in the Finite Element Method and Applications[END_REF]]. The error is defined by the energetic norm of the difference between a statically admissible field τ and the computed kinematically admissible u h field:

e rdc ( u h , τ ) = ||τ -H : ε( u h )|| H -1 ,Ω (2.104)
This choice is due to the fact that the constitutive equation can be considered as the less reliable.

The most interesting asset brought by this error estimator is the reliability of the estimation. For two admissible fields, an upper bound is guaranteed for η ex :

η ex = || u -u h || H,Ω ≤ e rdc ( ˆ u h , σh ) (2.105)
where σh is a statically admissible field to be constructed from σ h . The construction of an admissible stress field σh can be done with different more or less complex methods [START_REF] Beckers | Pure equilibrium tetrahedral finite elements for global error estimation by dual analysis[END_REF][START_REF] Parret-Fréaud | Estimation d'erreur de discrétisation dans les calculs par décomposition de domaine[END_REF].

Overall this category of estimators is considered as the most reliable, but the high cost of the construction of an admissible stress field is not negligible and the implementation is intrusive in a FEM software.

Gradient recovery discretization error estimators

Recovery procedures in error estimation were introduced by [START_REF] Olgierd | A simple error estimator and adaptive procedure for practical engineerng analysis[END_REF]. These procedures are widely popular because of their simplicity in terms of concept and also in implementation. They consist in approximating the discretization error by the difference between the finite element solution σ h and the so-called recovered solution σ * , which is constructed to be a better approximation of the exact solution then σ h .

||e * ,σ || = ||σ * -σ h || H -1 ,Ω (2.106)
In homogeneous media, let us recall that the approximated stress field σ h presents some discontinuities from an element of the mesh to another, while the exact stress field is continuous in homogeneous media. The recover field is constructed to be continuous, thus closer to the exact stress field σ than σ h . The construction of the continuous field σ * is based on the concept of super-convergent points, and is explained in the next section. Note that in case of heterogeneous media, this type of error estimator is not accurate because the exact stress tensor can be discontinuous.1 (i) Super-convergence and optimal recovery strategy The proof of existence of particular points in the mesh where stress and displacement fields are most accurate was discussed in [START_REF] Grätsch | A posteriori error estimation techniques in practical finite element analysis[END_REF]. In general, the best sampling points for displacement are the nodes of the mesh. For the gradient (like stress) field there are optimal points inside the elements (often, the integration points).

Let consider a typical second order equation for which the exact solution is assumed to be known. In figure 2.24, a comparison between FEM solution quadratic 1D element and the exact solution is displayed. The approximate solution is exact on the nodes of the elements. And for the gradient field (which can represent the strain or the stress), the approximate solution is exact on points inside the element. In the present case, these locations are Gauss integration points for 1D quadratic element.

The concept of super-convergent point has been generalized to 2D elements. For quadrangular elements, on the super-convergent points, the error decreases quicker than any other location within the element (the error is of the order of o(h p+1 )) 2 . On the other hand, for the triangular elements, it has been shown that the super-converging points do not exist, but one can refer to the table 2.25 for optimal sampling points and their degree of convergence.

The information available in 2.25 is not valid for distorted elements, the full study is available in [START_REF] Zienkiewicz | Superconvergence and the superconvergent patch recovery[END_REF].

(ii) Superconvergent patch recovery (SPR) The SPR method and its variants use the superconvergent points to construct the smoothed field σ * . Based on the SPR method, the ZZ2 discretization error estimator was introduced in [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF].

In practice, the smoothed field σ * is constructed by patch (see 2.26). Each component of the stress field is approximated on a patch by a polynomial function with the same degree as The coefficients a i of the polynomial function are solution of the following minimization problem: where (x k , y k ) are coordinate of the superconvergent points of the patch. The coefficients a i are immediately obtained:

Π = n k=1 [σ h,i (x k , y k ) -P k a] 2 P k = P(x k , y k ) (2.108)
a = n k=1 (P t k P k ) -1 • n k=1 P t k σh,i (x k , y k ) (2.109)
The polynomial function is identified for each patch using 2.109. For each vertex node, the smoothed approximation is computed directly using the polynomial function. In the case of quadratic meshes, the value of the field on the middle nodes is an average value of contributions computed from the patches sharing this node. The final smoothed approximation is obtained by using these nodal values in combination with the displacement shape functions.

The smoothed field σ * constructed using the SPR procedure can be used to compute an error estimator, the so called ZZ2 error estimator:

||e * ,zz2 || = ||σ * SP R -σ h || H -1 ,Ω =   Ω E ||e * ,zz2 E || 2   (2.110)
(iii) Adaptation to non linear problems An extension of the ZZ2 discretization error estimator to non-linear problems was proposed by [START_REF] Boroomand | Recovery procedures in error estimation and adaptivity. part ii: Adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF]]. An incremental measure of error adapted for non linear problem has been defined:

||e * ,incr || = Ω (σ * -σ h ).(∆ε * -∆ε h )dΩ 1 2 =   Ω E ||e * ,incr E || 2   1 2
(2.111)

In order to compute this incremental error estimator, two SPR recoveries are needed: a first one for the stress field σ and a second one for the increment of strain field ∆ε. According to [START_REF] Boroomand | Recovery procedures in error estimation and adaptivity. part ii: Adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF], this error estimator is also valid in finite strain framework.

The ZZ2 incremental error estimator can be compared to the global energy, computed using the smoothed field, in order to obtain a relative version of the error estimator:

η * ,incr = ||e * ,incr || E * Where E * = Ω σ * .∆ε * dΩ 1 2 (2.112)
This relative version is essential for adaptive procedures.

Since the principle of recovery procedure can be applied to any integration point field, tensor or scalar, a ZZ2 incremental error estimation can be applied to a scalar field f was proposed by [START_REF] Andrade Pires | On the finite element prediction of damage growth and fracture initiation in finitely deforming ductile materials[END_REF][START_REF] Feld-Payet | Finite element analysis of damage in ductile structures using a nonlocal model combined with a three-field formulation[END_REF]:

||e * ,incr scalar || = Ω (f * -f h ).(∆f * -∆f h )dΩ 1 2 =   Ω E ||e * ,incr scalar E || 2   1 2 (2.113) η * , incr scalar = ||e * ,incr scalar || E * where E * = Ω f * .∆f * dΩ 1 2
(2.114)

Conclusion

The gradient recovery discretization error estimators represents a good compromise for the mesh adaption purposes, as evidenced in Feld-Payet [2010]. Indeed, they involve a relatively small computational cost. And they can easily be extended to non-linear mechanical problem and can be either used for scalar or tensor fields. In our case of the nonlocal GTN model, the porosity or cumulative plasticity variables can be used in order to detect the high gradient areas i.e., the localization bands. Besides, since the discretization error uses the increment of a scalar field and damage hardly evolves in the already cracked areas or far from the crack tip, the error estimated by the incremental estimator should be small and therefore the mesh in these areas should be relatively coarse, once the error estimation is combined with a mesh an optimality criterion.

Mesh adaption procedure

An adaptive procedure should allow to ensures the respect of the desired precision during the computation with minimal cost. For this reason driving the adaptive procedure using a discretization error was chosen. To achieve this goal, both the global and local error estimations are used to generate an optimal mesh when the precision is unsatisfactory. Different criteria are presented in this section. For a more complete review see [START_REF] Díez | A unified approach to remeshing strategies for finite element h-adaptivity[END_REF].

Global criterion for mesh adaption

The global criterion consists in verifying that the global relative estimated error η do not exceed the threshold defined by the user ε 0 , other wise a mesh adaption is triggered. The criterion can be expressed by the following inequality:

η ≤ ε 0 (2.115)

Construction of a map size using the discretization error

The common motivation behind the local criteria is: how to build an "optimal mesh" based on a distribution of error which should be, in some sens, equally distributed and/or should enable to obtain the minimal elements number for the new mesh.

The different criteria presented bellow where constructed based on these hypotheses: 1-The mechanical problem is linear.

2-The field used to measure the error is regular enough so that the convergence of the error only depends on the type of the finite element:

η ≈ O(h q ) η E k ≈ O(h q )mes(Ω E ) 1 2 = O(h q+d/2 ) (2.116)
where q is the order of the shape function used in the mesh, and d is the dimension of the space.

For non-linear simulation, the convergence of the incremental ZZ2 error estimator is different from the convergence rate for linear problems (in 2.116). According to [START_REF] Boroomand | Recovery procedures in error estimation and adaptivity. part ii: Adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF], the actual rate of convergence is:

η = O(h q ) + O(h q/2 ) (2.117)
In practice, this modification is not taken into account for the construction of the local criterion.

All the local criteria try to determine the elementary size factor between the current and the new mesh sizes:

r E = h * E h E (2.118)
where h E and h * E are respectively the element size on the current and the next mesh, for a considered elementΩ E . Let us note that in the following all the quantities on the new mesh are all annotated with an asterisk (i.e η * , η * E , h * E ).

(i) Uniform repartition of the error on the current mesh There are different ways if enforcing the equal distribution of the global error, or more exactly of the square of the global error. For example, [START_REF] Oñate | Adaptive Mesh Refinement Techniques for Structural Problems[END_REF] have proposed to formulate the local condition on the current mesh:

η 2 E = η 2 N
where N is the number of elements in the current mesh (2.119) Writing equation 2.116 for the current and the new meshes, and using the hypothesis in 2.119, the elementary size factor can be written as:

r E = η 2 2q+d -1 q ε 1/q 0 η E • √ N 2 2q+d
(2.120) A simplified version is introduced in by [START_REF] Boroomand | Recovery procedures in error estimation and adaptivity. part ii: Adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF] coupled with the incremental ZZ2 error estimator:

r E = ε 0 η E √ N 1/q
(2.121)

(ii) Uniform repartition of the error on the new mesh [START_REF] Ladevèze | La maîtrise du calcul en mécanique linéaire et non linéaire[END_REF] chooses to express the equal distribution of error on the new mesh:

(η * ) 2 E = (η * ) 2 N *
where N * is the number of elements in the new mesh (2.122) In order to use the previous equation, an estimate for the new mesh size N * is needed:

N * = Ω E 1 r d E (2.123)
Following the same method to construct the previous criterion 2.120, the elementary size factor can be written as:

r E = ε 1/q 0 η 2 2q+d E Ω E η 2d d+2q E 1/2q
(2.124) (iii) Uniform repartition of the specific error on the current mesh [START_REF] Bugeda | A general methodology for structural shape optimization problems using automatic adaptive remeshing[END_REF] proposed to impose a uniform distribution of error divided by the volume of the element:

η E mes(Ω E ) = η mes(Ω) (2.125)
The elementary size factor can directly be expressed as:

r E = ε 0 mes(Ω E ) η E mes(Ω E ) 1/q
(2.126)

(iv) Minimal element number of the mesh This criterion proposes to minimize the number of elements in the new mesh (see [START_REF] Li | Notes on mesh optimal criteria in adaptive finite element computations[END_REF])):

L(r E i , λ) = Ω E 1 r d E + λ   Ω E r 2q E η 2 E -ε 0   (2.127)
The solution of this minimization can be written as: (2.128) This elementary factor is the same as the one obtain with the hypothesis of uniform elementary error distribution on the new mesh.

r E = ε 1/q 0 η 2 2q+d E Ω E η 2d d+2q E 1/2q

Conclusion

In the chosen adaptive procedure, the main goal is to obtain an optimal mesh in terms of precision and cost. According to [START_REF] Díez | A unified approach to remeshing strategies for finite element h-adaptivity[END_REF], the criterion based on uniform repartition of the error on the current mesh gives meshes where the target error is respected with acceptable computation cost, but the new mesh has the most significant element's number. Besides, the local criterion aiming for the minimal element's number and the criterion based on the uniform repartition of error on the new mesh, both lead to a new mesh with the least number of elements and a measured error on the new mesh very close to the target precision.

The chosen criterion will then be the one minimizing the size of the mesh and resulting in a uniform distribution of the error on the new mesh:

r E = ε 1/q 0 η 2 2q+d E Ω E η 2d d+2q E 1/2q
(2.129) Note: With the chosen local criterion, if the estimated error is too low, then the elementary size factor r E will diverge. In order to avoid this artifact, a maximum element size is imposed during the remeshing process.

Transfer operator

Transfer operators are necessary each time there is a change of discretization. Both nodal data (displacement, nonlocal variables, etc) and integration point data (stress, strain, internal variables) should be transferred onto the new mesh. In order to minimize the perturbation associated with the transfer operation on the computation after remeshing, the transfer field should satisfy the same constrains as the original fields, i.e respecting the local equilibrium, boundary conditions and transfer should preserve the localization band width. In practice, a transfer operator that satisfies all these constrains does not exist and an inevitable error of transfer is always present. Therefore compromises on the choice of the transfer operator have to be made depending on the type of problem.

In this section the most used transfer operators are briefly presented.

Nodal field transfer operator

For nodal fields, the transfer operator is straightforward. First, each node of the new mesh, having X new as coordinates, is localized on the old mesh3 , then the shape functions φ old i of the old mesh are used to compute the value of the field f on the new mesh:

f (X new ) = i φ old i (X new ).f (X i ) (2.130)
where X i is the coordinate of the node i of the corresponding element. Note: For computation including cracks, a special element locator's is needed for nodes on the closed crack lip.

Integration point data transfer operator

(i) Nearest Integration point For each integration point (IP) on the new mesh, the closest IP on the old mesh is searched, and the corresponding value of the field is transferred:

f (X ip new ) = f (X closest ip old
). This transfer operator is valid if the mesh does not change drastically before and after the remeshing procedure. Due to its simplicity, this operator has a negligible computation cost, and preserve high gradient fields after their transfer [START_REF] Simone | From continuous to discontinuous failure in a gradient-enhanced continuum damage model[END_REF].

(ii) Integration point transfer using extrapolation to nodes This transfer operation can be divided into three steps: (i) extrapolation of data from the integration points of the old mesh to the nodes of the old mesh [START_REF] Lee | Error indicators and adaptive remeshing in large deformation finite element analysis[END_REF][START_REF] Perić | Transfer operators for evolving meshes in small strain elasto-placticity[END_REF][START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF], (ii) classical nodal transfer from the old mesh to the new mesh. (iii) interpolation from nodes to integration points of the new mesh using the shape functions.

This method causes a significant diffusion compared to the nearest integration point transfer as reported in [START_REF] Perić | Transfer operators for evolving meshes in small strain elasto-placticity[END_REF]. This diffusion is not desirable especially for damage simulation where the size of the localization band needs to be preserved.

(iii) Nearest integration point with patch smoothing procedure This technique was introduced by [START_REF] Boroomand | Recovery procedures in error estimation and adaptivity. part ii: Adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF]. The first step consists in determining the element E old k of the old mesh containing the considered point IP new . For each node j of the element E old k , a patch P j can be constructed (the elements sharing the considered node). Let us not f , the field to be transferred. On each patch, a smoothed f * is computed using a weighted least square fit with a polynomial function.

The value at the new IP is computed as an weighted average of the nodal smoothed field computed from each patch P:

f (X IP new ) = j ω j f * j j ω j (2.131)
where ω j designates the weight associated with the node j. This type of methods suffer from diffusion of the field due the averaging step. The presented method using smoothing by patch, is considered in literature, the best method in terms of minimizing the diffusion problem due to the use of local patches.

(iv) Moving least square This method was introduced in [START_REF] Villon | Transfert de champs plastiquement admissibles[END_REF] to transfer the integration point fields based on approaching the discrete IP field with a continuous field which is defined locally as polynomial function of order p. In other word, the scalar IP field f (or each scalar component of an IP tensor field) is written in the new mesh at an IP X IP new as a polynomial function P and its parameters a as:

f (X IP new ) = P T a (2.132)
The parameters a of the polynomial function are then obtained by minimizing the following quantity:

J X IP new (a) = 1 2 i∈V (X IP new ) W (X i , X IP new )||P T (X i -X IP new i )a -f (X i )|| 2 (2.133)
where V (X IP new ) represent the vicinity of an IP X IP , and W (X i , X IP new ) correspond to a weight function defined in the vicinity V . This method is widely used [START_REF] Brancherie | Transfert de champs par approximation diffuse avec conservation de l'énergie[END_REF]; Feld-Payet [2010] due to its performance (the computation of the new field is done locally), and With the right choice of the size of the vicinity V for an integration point, the diffusion problem of high gradient field could be limited [START_REF] Patzák | Adaptive resolution of localized damage in quasi-brittle materials[END_REF].

Choice of transfer operator

For the transfer operator, there is no perfect choice valid for any type of problem. The choice of the operator may depend on the constitutive equations. In particular, in case of coupled damage models for which the localization band is an important object of study, the transfer operator should priorities the limitation of diffusion problems. The nearest integration point operator and its variation are valid candidates. However, the relevance of these method is only limited to cases where the mesh sizes do not change significantly.This is not the case with mesh adaption procedures based on error estimation. In this case of procedures, a more elaborate transfer operator as the moving least square per patch is more adapted, the patch size should be tuned in order to limit diffusion problems in case of localizing variables.

In this thesis work, different transfer operators are used depending on the remeshing procedure 4.5 Reequilibrium

Context

The transfer operation needed after mesh adaption creates inevitable transfer errors: neither the global equilibrium nor the constitutive equations are satisfied. Also as the transfer of nodal and integration point variables is independent the transferred strain does correspond to the nodal displacement. These in consistencies degrade the rate of convergence of after remeshing, or can result in a divergence of the resolution algorithm if the transfer creates a significant error.

In order to limit this problem, some measures are needed to facilitate the convergence of the resolution algorithms.

Load increment division

In order to facilitate convergence after remeshing, [START_REF] Boroomand | Recovery procedures in error estimation and adaptivity. part ii: Adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF] proposed to subdivide the load increment until the convergence is obtained. This technique compatible with most FEM resolution algorithms. However, a significant computational time may still be needed, especially with crack increment insertion, in order to find the right load increment that enables convergence. Indeed, for highly non linear problems, this load subdivision may lead to really small increment, thus the convergence rate after remeshing is greatly impacted.

In practice, this method can be used coupled with other more elaborate reequilibrium methods.

Local reequilibrium

In order to limit the transfer error, a minimal number of variable should be transferred [START_REF] Camacho | Adaptive lagrangian modelling of ballistic penetration of metallic targets[END_REF][START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF][START_REF] Boroomand | Recovery procedures in error estimation and adaptivity. part ii: Adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF][START_REF] Javani | Consistent remeshing and transfer for a three dimensional enriched mixed formulation of plasticity and non-local damage[END_REF] to compute the remaining ones on the new mesh. With this strategy, the computed variable satisfy at least local relations such as the constitutive equations. In theory The rate of convergence should be improved because the transfer errors are reduced. However, the authors proposing this method did not provide a thorough analysis of this reequilibrium method and its improvement on the rate of convergence.

Elastic step

Another approach is based on an additional elastic step in order to reestablish the global equilibrium [START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF]. The convergence of this step is somewhat easier than the normal step where plastic and damage variable can evolve. However, from the author experience, although this step may help the rate of convergence after mesh adaption, it does not guarantee the convergence of the next step.

Temporary boundary conditions

Another method consist in applying the force associated with th displacement field and gradually dissipating it over a time increment until the structure returns to equilibrium [START_REF] Boroomand | Recovery procedures in error estimation and adaptivity. part ii: Adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF][START_REF] Javanmardi | Extended finite element method and anisotropic damage plasticity for modelling crack propagation in concrete[END_REF]. This strategy has not been fully studied. From the author experience, in presence of a long crack in the computation mesh and using this reequilibrium strategy, the resolution algorithm may get stuck in the reequilibrium step because of often divergences and load increment subdivisions.

Comments

In order to benefit from the speed up of the proposed remeshing strategy coupled with the error estimator, an efficient reequilibrium strategy is needed. Because, in the absence of an additional correction step(s) and using only the simple load increment division, all the computation time gained from the optimized mesh can be wasted in too small load increments. All the presented method in literature have only been used in simple benchmark simulations.

Summary

Continuous damage models are widely used in the literature to simulate ductile failure with the local approach. These models can be either be qualified as uncoupled damage models, so that the material model is not affected by the evolution of the damage indicator, or qualified as coupled damage models, so that the material behavior is affected by the damage evolution. Coupled damage models are more realistic because the material behavior of a healthy material is different from a damaged material. However, the use of coupled damage models for the simulation of structures using FEM leads to the well-known mesh dependency that occurs during material softening. Thus, for local damage models, the mesh size is considered a parameter of the model.

So-called nonlocal models have been developed to overcome this mesh dependency. These nonlocal models introduce in the material behavior a characteristic length, equivalent to the specific mesh size used in the local simulations. This characteristic length could be interpreted as related with the interactions between the material defaults at the microscopic length scale. It can be introduced by different methods:

-(i) in additional terms to the Helmholtz energy either to enrich the kinematics of the continuous medium (micromorphic approaches) or to constrain the gradient of a local variable (energy-enriched approaches) in order to regularize the mechanical problem.

-(ii) as an additional equation that produces a smoothed form the local variables (integral methods, explicit/implicit gradient).

-(iii) as a regularization of the Griffith's functional energy, as in case of the phase field approach, in order to replace the discrete (crack) transition from the broken area to the intact area by a smooth transition with progressively damaged zone.

In order to use a nonlocal model to simulate the initiation and propagation of ductile cracks over long distances, two main limitations have to be addressed, (i) the continuous damage model cannot describe the kinematics of crack opening and closing, (ii) due to the necessary fine discretization in the localization band, the size of the problem increases with the length of the crack.

On one hand, the first limitation can be addressed by using the so-called continuousdiscontinuous approach, which consists, in the case of a ductile fracture, in representing the crack as a discontinuity in the sufficiently damaged areas. To carry out a continuousdiscontinuous approach, (i) a representation of the crack must be defined, (ii) a method of localization of the discontinuity is necessary to automatically define the position of a continuous crack surface, (iii) and an insertion criterion to select the most appropriate time to transit from the continuous model to the discontinuous model is essential.

On the other hand, the second limitation can be addressed by using efficient mesh adaption. In contrast to fixed-mesh modeling where a fine mesh is assigned in all areas likely to be crossed by a crack, resulting in a significant increase in computational size, only the regions of interest are finely meshed automatically during the simulation. To do this, error estimation can be used, in combination with a global criterion to trigger this adaptive remeshing, and with a local criterion is needed to specify the new mesh size.

In Figure 2.27, the proposed continuous-discontinuous procedure with mesh adaption is summarized. A nonlocal implicit gradient GTN model 4 with two characteristic lengths is Figure 2.27: Flowchart summarizing the proposed continuous-discontinuous procedure with a nonlocal GTN model and discrete crack lips representation using mesh adaption used to simulate ductile failure (see chapter 3) where the damage variable is the total porosity f . At the end of each converged load increment, the incremental error is evaluated. If it exceeds a specific threshold prescribed by the user, the mesh adaption procedure is triggered. Simultaneously, the marching ridges algorithm takes an input f * , the smoothed version of f obtained using the SPR method, and its computed spatial gradient ∆f * and searches for the crack direction to provide an auxiliary mesh describing the crack surface. Then an insertion criterion is tested to determine whether there should be crack initiation and/or propagation, if a crack is found. Mesh adaption is triggered if the insertion criterion is validated the insertion of the found crack increment.

The remeshing procedure, described in chapter 4, consists in generating a new mesh5 respecting the size map calculated using the mesh optimal criterion and possibly including the crack surface. Once the new mesh is generated using meshGems, the mechanical fields necessary to resume the calculation are transferred using the nearest integration point or the moving least square methods. A reequilibrium step is then performed to correct the transferred fields and thus improve the convergence rate after the remeshing procedure.

of the deformed mesh if this latter is adequately meshed 5 Remeshing is carried out on the updated frame so that the deformed structure is meshed with good quality element as opposed to stretched element found in case of the total Lagrangian framework in presence of high strain.

Introduction

Predicting ductile failure is a major challenge when trying to prevent catastrophic failure of structures or to control crack propagation in the case of metal forming. One of the main challenges is the prediction of complex crack paths in metals such as cup cone or slant fracture. To achieve this goal, models able to robustly represent the local degradation phenomena (damage nucleation, growth and coalescence) are needed. Among such models is the micro-mechanics based Gurson-Tvergaard-Needleman (or GTN) model [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]. Based on this seminal framework, numerous extensions have been proposed which are reviewed e.g. in [START_REF] Tvergaard | Material failure by void growth to coalescence[END_REF][START_REF] Besson | Damage of ductile materials deforming under multiple plastic or viscoplastic mechanisms[END_REF][START_REF] Benzerga | Ductile Fracture by Void Growth to Coalescence[END_REF][START_REF] Pineau | Failure of metals I: Brittle and ductile fracture[END_REF]. They are able to account for void shape, void orientation, plastic anisotropy, strain rate sensitivity, coalescence by internal necking. . . All these models must lead to softening up to a point where the material entirely loses its load-carrying capacity so as to be able to model fracture. These softening models were first developed within the usual local framework and thus suffered from the well-known numerical issues such as mesh size and mesh orientation dependence [START_REF] Rousselier | A Methodology for Ductile Fracture AnalysisBased on Damage Mechanics: Anlllustration of a Local Approach of Fracture[END_REF][START_REF] Liu | Mesh-dependence and stress singularity in finite element analysis of creep crack growth by continuum damage mechanics approach[END_REF][START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF]. Due to the loss of ellipticity of the equations in statics, the solution of the mechanical problem is no longer unique. Overcoming these issues is fundamental for the reliability of fracture simulations. Using a finite strain framework is also required as ductile fracture is always accompanied by large deformations. Several solutions have been proposed to overcome mesh dependency. A first solution is based on a local enrichment by embedding a finite thickness band [START_REF] Huespe | A finite thickness band method for ductile fracture analysis[END_REF]. In that case, the band is introduced when loss of ellipticity is detected within one element. The normal to the band is determined following Rice's bifurcation analysis [START_REF] Rice | The localisation of plastic deformation[END_REF]. Another solution is based on implicit gradient methods [START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF][START_REF] Geers | Strain-based transient-gradient damage model for failure analyses[END_REF]] which facilitate the use of integral methods as originally proposed in [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Bazant | Non local continuum damage. localization, instability and convergence[END_REF]]. These methods have been used to model ductile fracture for metals within a finite strain framework [START_REF] Enakoutsa | Numerical implementation and assessment of a phenomenological nonlocal model of ductile rupture[END_REF]Mediavilla et al., 2006b;[START_REF] Linse | Simulation of crack propagation using a gradient-enriched ductile damage model based on dilatational strain[END_REF][START_REF] Hütter | Simulation of ductile crack initiation and propagation by means of a non-local Gurson-model[END_REF][START_REF] Javani | Three-dimensional finite element modeling of ductile crack initiation and propagation[END_REF][START_REF] Seupel | On the identification and uniqueness of constitutive parameters for a non-local GTN-model[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] although they were initially developed for quasi-brittle failure. Using an implicit gradient formulation, the local constitutive equations are preserved at the cost of small adaptations. This makes this formulation particularly attractive in the case of the GTN model, as its initial micro-mechanical foundations are then preserved. Multiple internal lengths were introduced in [Nguyen et al., 2020]. This model uses an implicit gradient formulation written in the initial configuration. An internal length was associated with each damage evolution mechanism: (i) void growth based on the classical GTN model, (ii) internal necking governed by a heuristic extension of the Thomason model [Thomason, 1985a,b] based on the maximum principal stress, (iii) shear-dominated coalescence mechanisms controlled by the maximum shear stress [START_REF] Torki | On void coalescence under combined tension and shear[END_REF]. Although the model proposes different internal length scales, the same characteristic length was assigned to the three nonlocal variables in the presented simulations.

Micromorphic models [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF] can also be used to solve problems related to the loss of ellipticity. They use local state variables and their "micromorphic" counterparts. The free energy of the system depends on the gradient of the micromorphic variables and on coupling terms. Ductile fracture was modeled using this class of models in [START_REF] Brepols | Gradient-extended two-surface damage-plasticity: Micromorphic formulation and numerical aspects[END_REF][START_REF] Diamantopoulou | Micromorphic constitutive equations with damage applied to metal forming[END_REF]. In particular, models based on a microdilatational theory [START_REF] Huetter | A micromechanical gradient extension of gurson's model of ductile damage within the theory of microdilatational media[END_REF] appear to be well suited to represent ductile damage by void growth. More recently, nonlocal gradient enhanced energy (GEE) models, first developed for quasi-brittle fracture [START_REF] Lorentz | A variational formulation for nonlocal damage models[END_REF], were also used to model ductile failure [START_REF] Zhang | Ductile damage modelling with locking-free regularised gtn model[END_REF][START_REF] Chen | Properties of a nonlocal gtn model within the context of small-scale yielding[END_REF]. They only use the gradient of a local state variable. A decomposition-coordination technique is used to treat the nonlocality. The variable of interest is duplicated: a first instance is used at the (global) scale of the structure while a second instance is used at the (local) constitutive law level. As both variables represent the same field, they should be equal. A Lagrange multiplier is introduced to ensure this equality weakly. Finally, the regularization of ill-posed problems can also be achieved by coupling the elastoplastic models with phase-field formulations. This type of formulation was initially introduced for brittle failure [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Hofacker | Continuum phase field modeling of dynamic fracture: variational principles and staggered fe implementation[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Heider | A phase-field modeling approach of hydraulic fracture in saturated porous media[END_REF]. However, many extensions to ductile fracture have been introduced for 2D simulations [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF]Aldakheel et al., 2018], for 3D simulations [Aldakheel et al., 2018;[START_REF] Ambati | A phase-field model for ductile fracture at finite strains and its experimental verification[END_REF][START_REF] Michael | A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF][START_REF] Miehe | Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory[END_REF][START_REF] Hu | A variational phase-field model for ductile fracture with coalescence dissipation[END_REF] and also for 3D simulations with remeshing techniques [Eldahshan et al., 2021b]. More recently, a porous ductile model using a phasefield formulation was introduced in which the critical energy is decomposed into elastic and plastic contributions where the plastic part is described using a GTN model [START_REF] Dittmann | Phase-field modeling of porous-ductile fracture in non-linear thermo-elasto-plastic solids[END_REF].

This work proposes a nonlocal extension of the GTN material model (based on [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]) using an implicit gradient framework incorporating two nonlocal characteristic lengths to capture complex ductile failure patterns, such as cup-cone or slant crack paths. The most common issue with nonlocal formulations is the choice of the characteristic lengths. In this chapter, comparisons with experimental data on steel specimens first enable ontaining an order of magnitude for the model parameters. Then two length scales are distinguished in the model to account for the different spacings that exist between material defaults responsible for damage nucleation and void growth. In the case of ferritic steels used in modern pipelines, two damage mechanisms exist: (i) void growth from MnS and oxides and (ii) damage nucleation at iron carbides (Fe 3 C), which occurs at high strains. Mns and oxides early debond from the metallic matrix so that they can be considered as initial voids. Similar damage behavior is e.g. observed in aluminum alloys where coarse particles are Ferich inclusions and small particles are strengthening dispersoids [START_REF] Bron | Ductile rupture in thin sheets of two grades of 2024 aluminum alloy[END_REF]. In order to guide the choice of these characteristic lengths, an original study is conducted on their effect on the occurrence of cup-cone and slant fracture. This work is structured as follows. First, in section 2, the local version of the model is briefly recalled and its nonlocal modification is introduced. The proposed finite strain formulation is also introduced within an updated Lagrangian framework. In section 3, a parameter fit is carried out for both the local and nonlocal version of the GTN model based on existing literature data. In section 4, the independence of the results to both mesh size and mesh orientation is checked using the same characteristic length for both mechanisms. The relationship between the numerically obtained localization band width and the characteristic lengths is also investigated. This allows the selection of appropriate mesh sizes for given characteristic lengths. Finally, in section 6, an investigation is presented of the effect of different characteristic lengths for damage nucleation and void growth on the formation of cup-cone and slant crack paths and the results are translated in terms of size effects.

Material models and their framework

In the first part of this section, the local version of the considered GTN model used throughout the chapter is briefly recalled. From this basis, the proposed modifications to obtain a nonlocal model with two characteristic lengths are developed. A corotational finite strain formulation is used to deal with large strains (see section 2.2). This allows to use a simple small strain like formulation for the constitutive model based on an additive strain decomposition of the strain tensor.

Material models

The used local model is based on the GTN model presented by [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] but with a different nucleation function using the accumulated plastic strain κ as opposed to the original model which uses the growth porosity f g .

Local GTN model

One assumes the additive decomposition of the strain tensor (ε) into an elastic (ε e ) and a plastic (ε p ) part. The elastic strain tensor and the stress tensor are related using Hooke's law:

σ = E : ε e (3.1)
where E is the fourth-order elasticity tensor. Work hardening is assumed to be isotropic, and the flow stress of the sound material (R(κ)) is expressed as a function of the accumulated plastic strain κ assuming a power law:

R(κ) = K(e 0 + κ) n (3.2)
e 0 , K and n are parameters to be identified (any other function could indeed be used). The yield surface is then expressed following [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] as:

Φ = σ -R(κ)
where σ is an effective stress measure depending on both the stress tensor and the void volume fraction, which is the only damage parameter. In the case of the GTN model, the effective scalar stress is implicitly defined as:

G(σ, f, σ ) = σ 2 eq σ 2 + 2q 1 f cosh q 2 2 σ kk σ -1 -q 2 1 f 2 def = 0 (3.3)
where σ eq is the von Mises stress invariant1 of the Cauchy stress, σ kk designates its trace, q 1 and q 2 are two parameters and f is a function of porosity. Function f is defined as in [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]:

f =    f if f ≤ f c f c + δ(f -f c ) otherwise (3.4)
where f c represents the porosity at the onset of void coalescence and δ ≥ 1 is a parameter which represents the increased deleterious effect of porosity above f c . The damage variable f can be itself decomposed into two parts which reflect the fact that damage is caused by damage nucleation (f n ) and void growth (f g ), so that the total damage is equal to f = f g +f n .

The evolution of these damage variables is given by the following set of equations:

ḟg = (1 -f g )trace( εp ) mass conservation (3.5) ḟn = A n (κ) κ strain controlled nucleation (3.6)
where A n is a coefficient representing the damage nucleation rate. It is expressed as a function of κ so as to represent strain-controlled nucleation [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF][START_REF] Zhang | A complete Gurson model approach for ductile fracture[END_REF]. Note that eq. 3.5 slightly differs from the original one where (1 -f g -f n ) is used and not (1 -f g ). As MnS inclusions early debond from the matrix, it is unnecessary to describe nucleation of voids on these inclusions. In this work, nucleation correspond to damage created on iron carbides which are much smaller than the MnS inclusions. It is assumed that this creates very small voids but high damage. For that reason, the equation for ḟg is modified as growth is only attributed to large voids originated at MnS inclusions.

In practice, this modification has little effect on the overall behavior. The plastic strain rate tensor εp is obtained using the normality rule:

εp = (1 -f ) κ ∂Φ ∂σ = (1 -f ) κ∂σ ∂σ = (1 -f ) κn (3.7)
where n designates the normal to the yield surface. Using this expression, the following equivalence between the macroscopic (left handside) and the microscopic (right handside) plastic dissipations is obtained as:

εp : σ = (1 -f ) κσ (3.8)
as σ is an homogeneous function of order 1 of σ.

It is assumed that the considered material can exhibit a slight strain rate dependence, so that κ is expressed as:

κ = F(Φ) = ṗ0 σ -R σ 0 n (3.9)
where . is the positive part function and n, ṗ0 and σ 0 are material parameters.

Nonlocal GTN model

The proposed extension of this GTN model uses the implicit gradient methodology proposed in [START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF][START_REF] Geers | Strain-based transient-gradient damage model for failure analyses[END_REF][START_REF] Engelen | Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour[END_REF] to regularize two state variables: the volume variation ω = trace(ε p ) and the effective accumulated plastic strain κ.

Their nonlocal counterparts are referred to as ω and κ. Their evolution within the considered material body Ω is governed by the following Helmholtz-type equations:

ω -l 2 ω ∆ω = ω in Ω (3.10) κ -l 2 κ ∆κ = κ in Ω (3.11)
with the following natural boundary conditions:

∇ω. n = 0 on ∂Ω (3.12) ∇κ. n = 0 on ∂Ω (3.13) where ∂Ω and n respectively designate the boundary of the body Ω and its outer normal vector.

Let us note that the gradients are calculated on the current configuration and not on the initial configuration as considered by [START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF]Nguyen et al., 2020]. Besides, this model introduces two characteristic lengths in equations 3.10 and 3.11, l ω and l κ , respectively associated to void growth and nucleation mechanisms.

The nonlocal variables ω and κ are then used to formulate the evolution of the damage variables f g and f n as:

ḟg = (1 -f g ) ω (3.14) ḟn = A n (κ) κ (3.15)
These two equations now replace equations 3.5 and 3.6. All other equations remain unchanged. In particular, the hardening law R(κ) remains identical, which would not be the case with micromorphic models [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF] or gradient enhanced energy models [START_REF] Chen | Properties of a nonlocal gtn model within the context of small-scale yielding[END_REF] as it is expressed as a function of κ in those cases. In the present case, the flow stress then remains unaffected by the gradient of κ.

Implicit local resolution

The state variables describing the material behavior are then the elastic strain tensor (ε e ), the accumulated plastic strain (κ), the porosity due to void growth (f g ), the porosity due to damage nucleation (f n ) and the volume variation ω. Their evolution laws can be expressed as functions of the rates of the input variables ε, ω, κ as:

εe = ε -(1 -f g ) κn (3.16) κ = F(Φ) (3.17) ḟg = (1 -f g ) ω (3.18) ḟn = A n (κ) κ (3.19) ω = trace ((1 -f g ) κn) = (1 -f g ) κtrace(n) (3.20)
The implicit resolution of this system of equations for each integration point, detailed in A, requires its integration over a finite time step ∆t. This leads to a system of equations relating the residuals of these equations to the increments of the state variables ∆V S = (∆ε e , ∆κ, ∆f g , ∆f n , ∆ω), for a given set of the increments of input variables ∆V IN = (∆ε, ∆ω, ∆κ). The resolution of this system directly gives the increments of output variables [∆σ, ∆ω, ∆κ]. It then becomes possible to numerically compute the consistent tangent matrix by considering that any given small variation of the input variables δV IN leads to a modification of the state variables so that the residual equations remain null (see details in B). The consistent tangent matrix can be expressed as a block-matrix as:

∆V OUT =
K mat = ∂∆V OUT ∂∆V IN =           ∂∆σ ∂∆ε ∂∆σ ∂∆ω ∂∆σ ∂∆κ ∂∆ω ∂∆ε ∂∆ω ∂∆ω ∂∆ω ∂∆κ ∂∆κ ∂∆ε ∂∆κ ∂∆ω ∂∆κ ∂∆κ           (3.21)
The terms of this matrix are required to express the elementary stiffness matrix of the associated finite elements.

Finite strain framework

Both local and nonlocal GTN models must be formulated within a finite strain framework. This is done using a rotating frame concept.

(i) Corotational framework To do so, a corotational formulation [START_REF] Sidoroff | Some issues about anisotropic elastic-plastic models at finite strain[END_REF] is used. The corotational frame is defined using the rotation tensor Q such that:

Q = W .Q with Q(t = 0) = 1 (3.22)
where the W tensor designates the skew-symmetric part of the velocity gradient tensor L.

The symetric part of the velocity gradient tensor is denoted D so that: L = D + W . This velocity gradient tensor is defined from the transformation tensor F = ∂ x/∂ X (where X is the position of a material point in the initial configuration and x its position in the current configuration) and its derivative as: L = Ḟ .F -1 . Knowing the rate of deformation D, one gets the material deformation rate through the following expression:

ε = Q T .D.Q (3.23)
Knowing the corotational Cauchy stress (i.e., the Cauchy stress expressed in the coordinate system that rotates with the material), the Cauchy stress in the unrotated frame is given by:

Σ = Q.σ.Q T (3.24)
Let us note that the rotated rate of σ is the Jaumann rate of Σ.

Proof:

Σ = Q.σ.Q T + Q. σ.Q T + Q.σ. QT = W .Q.σ.Q T + Q. σ.Q T + Q.σ.Q T .W T = W .Σ -Σ.W + Q. σ.Q T using W T = -W .
(ii) Finite element formulation Considering cases where the load is applied sufficiently slowly so that the inertial forces can be neglected, the weak form of the equilibrium equations then reads:

∀ v Ω Σ : L dΩ = Ω Σ : D dΩ = W ext , with L = ∂ v ∂ x (3.25)
where v is a virtual velocity field, W ext represents the virtual power of external forces and D is the symmetric part of L .

In the case of a nonlocal formulation, based on equations 3.10, 3.11 and 3.12, the weak form corresponding to the nonlocal variables is expressed as:

∀ ẇ Ω (w -w) ẇ + l 2 w ∇w. ∇ ẇ dΩ = 0 ∀ ṗ Ω (p -p) ṗ + l 2 p ∇p. ∇ ṗ dΩ = 0 (3.26)
where ẇ and ṗ designates a virtual rate of the nonlocal field. Gradients are evaluated in the current configuration as ∇• = ∂•/∂ x.

These equations must then be discretized both in time and in space, as detailed in C. For the spatial discretization, quadratic shape functions are used for displacement degrees of freedom (DOFs), and linear shape functions are used for the nonlocal variables (The notation P2P1P1 is used in this chapter to refer to this choice of interpolation order).

The discretization of the set of equations 3.25 and 3.26 leads to a global system built thanks to the assembly of elementary reactions. This system relates the unknown variables to the external forces through the assembly of elementary stiffness matrices. The calculation of the different terms of these elementary matrices requires the evaluation of the consistent tangent matrix K mat , as detailed in D. The resolution of this system is performed implicitly using a Newton-Raphson algorithm and leads to the simultaneous determination of the DOFs associated with the displacements and the nonlocal variables (as opposed to a staggered resolution scheme).

3 Fitting of the model parameters

Model parameters

The constitutive model of [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF] being relatively close to the local one presented in this chapter, it is possible to use the same elasto-visco-plastic parameters and most of the GTN parameters (see table 3.1). The low initial porosity (f 0 ) corresponds to that of a modern line pipe steel. It corresponds to the MnS inclusion volume fraction as these particles easily debond from the matrix. Table 3.1: Set of material parameters from [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF] Young modulus E 210 GPa Poisson ratio ν 0.3 Isotropic hardening (power law) K 795 MPa e 0 0.002 n 0.13 Gurson criterion q 1 1.5 q 2

1. f 0 1.5 10 -4 Viscosity ṗ0 1s -1 σ 0 55 MPa n 5

The main difference between the constitutive model of [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]] and the local model in this chapter is the damage nucleation rate function A n (κ) in equation 3.6. In the model described in [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF], A n (κ) is a function of the growth porosity f g , whereas in this chapter, this function is chosen to be a constant (A n ) when the accumulated plastic strain reaches a threshold κ c :

A n (κ) =    A n if κ > κ c 0 otherwise (3.27)
The model was changed as the initial version did not allow forming cup-cone fracture using the nonlocal formulation. This leaves thus two material parameters, A n and κ c , to be fitted.

Besides, the nonlocal model requires the identification of the two nonlocal characteristic lengths l ω and l κ .

Tests used for the fit of the nucleation law

The experimental database consists in an axisymmetric tensile test presented in [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] and a plane strain test described in [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]. Both tests were carried out on the same material, i.e., X70 HSLA (high strength low alloyed) ferritic-pearlitic steel, and are experimental results from [START_REF] Rivalin | Développement d'aciers pour gazoducs à haute limite d'élasticité et ténacité élevée : Mécanique et mécanismes de la rupture ductile à grande vitesse[END_REF]. They have been selected because they result in complex crack paths, with cup-cone fracture for the axisymmetric specimen and slant fracture for the plane strain specimen. In addition, the tests have been used in the literature [START_REF] Scheider | Simulation of cup-cone fracture using the cohesive model[END_REF][START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] as a reference to identify and study ductile fracture models. Finite element simulations are carried out for both tests in order to fit the parameters (trial and error procedure). Meshes, dimensions as well as boundary conditions are shown in fig. 3.1a and fig. 3.1b. Mesh design is similar in both cases. In the case of the plane strain specimen, the entire specimen is meshed to enable the formation of a slanted crack path. In both cases, elements have quadratic shape functions with eight nodes and reduced integration (4 integration points). In the following S 0 represents the initial minimal cross-section of each specimen; S 0 = π /4φ 2 0 for the axisymmetric specimen (with φ 0 the initial diameter at the center of the specimen); S 0 = e 0 w 0 for the plane strain specimen where e 0 represents the specimen thickness and w 0 the width of the specimen. The normalized force (F/S 0 ) is then plotted as a function of the minimum diameter variation (∆Φ/Φ 0 ) or thickness variation (∆e/e 0 ).

Parameter fit for the local model

In the case of local models, it is important to keep the same mesh size in areas where cracks propagate [START_REF] Liu | Mesh-dependence and stress singularity in finite element analysis of creep crack growth by continuum damage mechanics approach[END_REF][START_REF] Skallerud | Effects of finite element mesh on the numerical prediction of ductile tearing[END_REF]. Because the diameter of the tensile bar (10 mm) is twice the thickness of the plane stress specimen (5 mm) and because symmetry is not accounted for in the latter case, the same number of elements (N h ) is used to discretize the minimum cross-sections in order to keep the same mesh size. The initial aspect ratio of these elements (r h ) is set to 6:1 so that it leads to approximately square elements at the onset of fracture. Using N h = 60, the initial mesh size is consequently: 83µm × 14µm.

The fitted material parameters are:

A n = 0.2 and κ c = 1.2 (3.28) Indeed, using these values, the global responses for both test cases show a good agreement with both reference results, as can be seen in fig. 3.2. In the case of the axisymmetric specimen (figure 3.2a), it can be noted that the fitted local model response is identical to the reference one before crack initiation, which corresponds to the sharp load drop. In this regime, the damage has minimal effect on the overall behavior. Crack initiation occurs at a slightly higher diameter reduction (-∆φ/φ 0 = 0.453) compared to the reference (-∆φ/φ 0 = 0.44). Besides, as load drop is faster for the new fitted parameters, the total [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] and [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]) indicates a correct parameter fit. Total porosity is displayed at total failure for both simulations: a cup-cone failure (left) and a slant failure (right) can be respectively observed.

diameter reduction (-∆φ/φ 0 ) c = 0.467) is ultimately slightly smaller than in the case of the reference simulation (-∆φ/φ 0 ) c = 0.48. The difference is less than 3 % and is acceptable. As can be seen on the image of total porosity at the end of the simulation (fig. 3.2a), a cupcone crack path is obtained.

In the case of the plane strain specimen (fig. 3.2b), the computed forces before crack initiation are slightly higher than the reference ones. This can be explained by the fact that in the reference simulation [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF], a Hill anisotropic model was used to describe plasticity. Such behavior is not accounted for in the present material model. In this case too, the thickness reduction (-∆e/e 0 ) c = 0.59 predicted using the new model parameters is slightly different from the reference (-∆e/e 0 ) c = 0.566. As in the case of the tensile test, the difference is acceptable. As can be seen on the image of total porosity at the end of the simulation (fig. 3.2b), a V-shaped crack path is obtained. Note that both V-shaped and S-shaped crack paths can be experimentally obtained [START_REF] Besson | Modeling flat to slant fracture transition using the computational cell methodology[END_REF].

Parameter fit for the nonlocal model

Using the local model, the band width is about the element height in the deformed configuration. For the above-considered cases, this size is l local b ≈ 70µm. The characteristic length for the nonlocal model was then selected so that the resulting band width (l nl b ) is approximatively equal to l local b . It is therefore assumed that the mesh size for the local model allows representing the material internal length following the early work by [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF]. Having l nl b ≈ l local b is a simple way to establish a correspondence between both models. Is is shown below (section 5.3) that the band width is related to the nonlocal characteristic length by the following relation:

l nl b ≈ 1.5l c .
Using this relation, one gets l ω = l κ ≈ 40µm. Using these values, 90 elements are required to discretize the specimen radius or thickness with at least three elements in the band width. Using fewer elements results in simulations being mesh dependent. This results in an initial element height equal to 10µm. As for the damage nucleation rate parameter A n , since the nonlocal model will always lead to higher ductilities than the local model with h m in different that the resulting localization bandwidth for a given set of material parameters, it was necessary to impose a larger value: A n = 0.4. Nucleation parameters for the local and nonlocal models are compared in table 3.2. Note that the number of elements used to mesh the minimum crosssection differs using the above hypotheses. Using N h = 60 for the nonlocal model would result in slightly mesh-dependent simulations.

This fit leads to a relatively good agreement of the global responses for both test cases with the reference results as well as results from the literature [START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF], as can be seen in figures 3.3 and 3.4. For the simulation of the smooth axisymmetric bar (see figure 3.3), the engineering stress-diameter reduction curve is particularly close to the reference from [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF], and the total failure diameter reduction is well predicted. As can be seen on the images of total porosity taken at different stages of the simulation on the same figure, a cup-cone crack path is also obtained. et al., 2001;[START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] and the local model indicates a correct parameter fit. Total porosity is displayed at various times of the simulation to illustrate the formation of the expected cup-cone crack path.

In contrast, the global response for the plane strain simulation (presented in figure 3.4) shows a slight overestimation of the strain at crack initiation compared to the reference and a higher total failure strain. This difference can be explained by the fact that the crack path remains flat (as can be seen on the images of the total porosity maps taken at different stages of propagation). This result is not in agreement with experimental data as well as the simulation using the local framework, which both exhibit slant fracture. The effect of the characteristic lengths on the crack path is discussed below in section 6 where it is shown that slant fracture may indeed be obtained using the nonlocal framework and a proper choice for the characteristic lengths.

Study of the convergence properties of the nonlocal model

In this section, the convergence of the solution with respect to mesh size and mesh orientation is first checked. This allows to determine the mesh size needed to obtain convergence as a function of the internal lengths, which are assumed to be equal. Using the appropriate mesh size, it is then possible to establish a relation between the internal lengths and the numerically obtained band widths. [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF][START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] and the local model indicates a correct parameter fit. Total porosity is displayed at various steps of the simulation to illustrate the propagation of a flat crack instead of a slanted one. This illustrates the need to properly calibrate the nonlocal characteristic lengths.

Effect of mesh size

Several meshes with the same element initial aspect ratio r h = 6 and the same orientation but different element sizes are used to study the effect of mesh size.

(i) Local model In the case of the local model, the considered set of number of elements in the width of the central section N h is: 45, 60, 75 or 90. The engineering stress (F/S 0 ) is plotted as a function of the variation of the minimum diameter -∆φ/φ 0 for different mesh sizes in figure 3.5. On the same graph, the distribution of total porosity f at total failure is shown for the different meshes. Let us note that all the simulations predict crack initiation for the same diameter reduction and that the diameters at full failure (F/S 0 = 0) only slightly differ. In this particular case, the well-known mesh size dependency (see e.g. [START_REF] Liu | Mesh-dependence and stress singularity in finite element analysis of creep crack growth by continuum damage mechanics approach[END_REF][START_REF] Rousselier | A Methodology for Ductile Fracture AnalysisBased on Damage Mechanics: Anlllustration of a Local Approach of Fracture[END_REF]) has no significant effect on the global response as the load drop is very sharp. Indeed, whatever the mesh size, the highly damaged area is localized within one row of elements.

(ii) Nonlocal model In the case of the nonlocal model, the considered set of number of elements in the width of the central section N h is: 20, 45, 60, 75, 90,, or 180. In figure 3.6, the engineering stress-diameter reduction curves and the total porosity f at total failure are shown for the different meshes. The global response appears to be converged, and the simulation leads to a cup cone crack path as soon as N h ≥ 45. Let us note that, for the coarsest mesh (N h = 20), the global response is very different than for the other meshes close to full failure. In that case, there is no crack bifurcation so that the minimum diameter always plastically deforms. In the case of cup-cone failure, the material undergoes elastic unloading at the notch root due to crack bifurcation at some stage of loading so that the minimum diameter does no longer change. Obviously, for N h = 20, the mesh is too coarse to properly capture the width of the highly damaged zone, and the simulation is, in that case, mesh size dependent. In all cases where cup-cone fracture is obtained, there exist two possible symmetric crack paths when flat propagation ends. One path is selected, but the other path can still be visualized. In the case of the smallest mesh size (N h = 180) a zigzagging crack path is obtained. In that case, the crack first runs downwards (-45 • ) and then upwards (+45 • ). Such a behavior is often experimentally observed [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF].

Effect of mesh orientation

In this section, three meshes having the same number of elements in the central part and same initial ratio r h = 6 : 1 are used. However the mesh is slightly tilted with respect to the horizontal axis. The tilt angles are: 0 • , 5 • and 10 • . Note that because of large strain, the final tilt angle is much larger than the initial one.

(i) Local model In the case of the local model, the number of elements in the central part is N h = 60. In figure 3.7, the global responses and the total porosity f at total failure are shown for the different meshes. Although the global responses are very similar, it is interesting to note that each mesh leads to a different crack path. Indeed, when there is not tilt, a classical cup-cone fracture path is obtained as shown above. For a tilt angle of 5 • , the crack first runs upwards from the center of the specimen and then runs downwards after reaching the mid-radius. For a tilt angle of 10 • , the crack runs upwards following one row of elements. In both cases, cup-cone fracture is not observed and the predicted crack path is only formed by cones. The local model thus clearly suffers from a strong dependence on mesh orientation. This implies that the crack path can hardly be predicted using local models. They can be useful if the predicted crack path corresponds to the expected one. In that case, the local stress and strain history can be investigated to analyze/understand crack formation. This can be done by prescribing the crack path using the computational cell methodology [START_REF] Xia | Ductile crack growth -I. A numerical study using computational cells with microstructurally-based length scales[END_REF][START_REF] Besson | Modeling flat to slant fracture transition using the computational cell methodology[END_REF]. The result also indicates that automatic meshing can hardly be used to generate the mesh because its design cannot be fully controlled. Simulations with automatic remeshing during the calculation are also not possible. (ii) Nonlocal model In the case of the nonlocal model, the number of elements in the central part is N h = 90. In figure 3.8, the engineering stress-diameter reduction curves and the total porosity f at total failure are shown for the different meshes. All the simulations give the same global response and lead to a similar cup-cone crack path. In all cases, the band width of the damage localization band is the same. Therefore, it can be concluded that results are converged and insensitive to mesh orientation. This paves the way for the use of automatic meshing techniques of complex parts and remeshing.

Analysis of the damage localization bands

In this section, the relation between the characteristic length and the width of the localization band is investigated. To do so, it is first necessary to define how to measure the width of a localization band. The finite element discretization needed to obtain a converged band width is then investigated. Based on these results, the relation between the characteristic length and the band width is studied for both axisymmetric and plane strain specimens as this relationship may depend on the stress state.

Definition of the band width

In the literature, there is no agreement, to the author knowledge, on a definition of a measure of the width of the localization bands. [START_REF] Mazière | Numerical aspects in the finite element simulation of the Portevin-Le Chatelier effect[END_REF] study localized bands in a nickelbased superalloy which exhibits a Portevin-Le Chatelier effect. Their proposed method is based on the analysis of the local plastic strain rate. It is applied to FE simulations but is designed to mimic actual experimental techniques (see references in the chapter). In [START_REF] Labergere | Strain rate distribution and localization band width evolution during tensile test[END_REF], the authors study a copper sheet and also use the strain rate profile The global responses are similar for all the orientations. However, the total porosity for each mesh size at total failure shows a different crack path for each mesh orientation, which indicates mesh dependence.

to determine the band width. The profile is fitted by a pseudo-Voigt function for which the band width is a fitting parameter. Gorodetskyi et al. [2017] propose a method able to detect discontinuities such as displacement jumps which can be applied to both FE simulations and experimental DIC results. In [START_REF] Mehenni | Spatiotemporal correlations in the Portevin-Le Chatelier band dynamics during the type B -type C transition[END_REF], plastic strain localization is experimentally studied as a function of the applied strain rate in the case of an Al-Mg alloy exhibiting a PLC effect. The bandwidth is defined "as the width at middle height of the localized strain band". In this study, it is proposed to define the band width (referred to as l v b for variable v = ω or v = κ) as the width over which the considered variable reaches half its maximum value. In that sense, this definition corresponds to that of [START_REF] Mehenni | Spatiotemporal correlations in the Portevin-Le Chatelier band dynamics during the type B -type C transition[END_REF]. In the following, the case of the axisymmetric specimen will only be described as similar results are obtained for the plane strain specimen. Profiles in the deformed configuration along the symmetry axis for both ω and ω are exemplified in fig. 3.9 where the band width for ω is shown by the blue arrows and the band width for ω by the red arrows. The profiles for ω can easily be plotted as this variable is defined at nodes and can be interpolated. Plotting the profiles for the local variable ω first requires extrapolating values at Gauss points to nodes (red curve in figure 3.9). In all cases, the band width for the local variable is about the size of one single element in the deformed configuration. In the following, the width of the localization bands is determined when the maximum of ω reaches 0.5 along the profile. In that case, the material points are close to full failure. Applying this definition to measure the band width for the local and nonlocal volume variations along the vertical symmetry axis for the axisymmetric specimen leads to consistent relative measures. Indeed, as shown in fig. 3.9, the band width measured on ω is larger than the band width measured on its local counterpart ω, even if the maximum of the nonlocal variable is lower than the maximum of the local variable, as expected using an implicit gradient nonlocal formulation. Besides, the total porosity for each mesh size at total failure shows an identical cup-cone crack path for each mesh orientation. These observations indicate mesh independence.

Conditions to obtain a converged band width

In this part, the evolution of the band width with element size is studied. As the present model makes use of gradients of ω and κ computed in the current configuration, it is important to evaluate the localization band width in the deformed configuration as proposed above. The element height after deformation and when fracture starts is referred to as H r1 in the following (see fig. 3.10). The element height along the element row just above is referred to as H r2 .

Due to different loading histories, the element elongation at fracture can differ greatly from one simulation to another. It is therefore important that the initial element height (h min ) is small enough to assure that the element height at fracture (H r1 ) is still able to represent the localization band. The convergence of the band width should therefore be checked with respect to H r1 . Different values for the number of element along the thickness of the specimens are used: N h = 540, 270, 180, 90, 60, 45, 30, 20. This study is performed for two values of the characteristic length, taken here equal for both nonlocal variables (l κ = l ω = 80µm and l κ = l ω = 40µm). Both axisymmetric and plane strain specimens are considered. The considered variable is the nonlocal volume variation but similar results can be obtained if the nonlocal effective accumulated plastic strain is used to measure the localization band width.

Figure 3.11 shows the evolution of the band width size as a function of the element height in the current configuration (H r1 ), as defined in fig. 3.10. Each set of points can be approximated by a line, which means that it can be assumed that the band width is an affine function of the mesh size. Using this approximation, the band width for H r1 → 0 can be extrapolated to a non-zero value which is an estimated value of the band width l ∞ b free from discretization error. In practical applications, this error is always present.

Relation between the band width and the internal lengths l κ and l ω

Let us note that in fig. 3.11, for a given geometry, the line corresponding to the largest characteristic length (i.e. l κ = l ω = 80µm) has approximately the same slope as the line corresponding to a smallest characteristic length (i.e. l κ = l ω = 40µm) but leads to a higher band width estimate l ∞ b . In order to validate this observation for other characteristic lengths and to establish a simple relation between the characteristic length and the band width l ∞ b , the same series of simulations was performed for additional values for of l ω = l κ : i.e. 20µm, 60µm, and 100µm. Plotting the evolution of the measured band width l ∞ b versus the imposed characteristic length (see fig. 3.12), it can be observed that for a given geometry the measured points approximately lie on a line. Using this regression, it is observed, as expected, that a null band width is obtained when the characteristic length is also null. Let us note that the slope in this linear relationship is different for each geometry. However, the difference between the relation obtained for plane strain and axisymmetric elements remains small: l ∞ b ≈ 1.54l κ,ω in the first case and l ∞ b ≈ 1.69l κ,ω in the second case. This simple relation presents a major advantage as it enables to estimate, for a given geometry, and based on the choice of the characteristic length, the width of the localization band. It is then possible to choose the number and size of the elements in the localization band (respectively denoted a and H r1 ) in order to properly capture its gradients with: a × H r1 < l ∞ b . In all the nonlocal simulations of this chapter (except for fig. 3.6), the factor a is at least equal to 3, i.e. there are at least three P2P1P1 elements in the band width at the onset of failure. .12: Evolution of the estimated band width free from discretization error for both the axisymmetric and plane strain tensile tests as a function of the characteristic length (l κ = l ω ). There appears to be a proportional relation between the band width and the characteristic length, with slightly different coefficients for each geometry. This means that it is possible to estimate the size the band width for a given characteristic length and a given geometry.

Using two characteristic lengths to model the cup-cone and slant fracture

In this section, two different characteristic lengths (l κ ≤ l ω ) are used to model the occurrence of cup-cone and slant fracture.

Procedure

Several simulations were performed for different characteristic lengths varying between 20µm and 100µm and different ratios l ω /l κ (i.e. 1, 2, 3, 4 and 5). Let us note that larger ratios were not considered in order to avoid too large computational costs due to mesh refinement associated with the requirement to have 3 elements in the band width at the onset of fracture

(3 × H r1 < l ∞ b )
. The relative values of the two characteristic lengths are chosen to be consistent with the physical degradation of metallic materials, and especially steels. Indeed, in steels, two damage mechanisms exist: (i) void growth from sulfides/oxides and (ii) damage nucleation at iron carbides Fe 3 C which occurs at high strains [START_REF] Tanguy | Plastic and damage behavior of a high strength X100 pipeline steel: experiments and modelling[END_REF]. As the spacing between sulfides/oxides is larger than between carbides, it makes sense to only consider l ω ≥ l κ . Figure 3.13: Contours of total porosity at total fracture for different characteristic lengths l κ and l ω for the axisymmetric specimen. Cup-cone crack paths are observed only for l κ < 80µm. The characteristic length l κ should then be relatively small compared to the specimen size to observe a cup-cone fracture. The larger the characteristic lengths, the wider the localization band (which is linked to the displacement at failure) and the further from the symmetry axis the bifurcation..

Cup-cone fracture

For the axisymmetric specimen, the contours of total porosity at full fracture for different characteristic lengths are displayed in fig. 3.13. It appears that large values for l κ prevent the occurrence of crack bifurcation for the considered specimen. As seen in fig. 3.13, simulations with l κ ≥ 80µm predict a flat crack while all simulations with l κ < 80µm predict the cupcone crack path. For a given value of l κ , it can be noticed that the larger l ω , the wider the localization band. Similarly, for a given value of l ω , the band width increases with increasing l κ . This is indeed expected since a large characteristic length leads to a larger band. In addition larger bands tend to increase the length of the flat part of the crack path. The flat central crack then bifurcates to form the slanted crack path at +45 • or -45 • . In one case (l ω = 60µm and l κ = 20µm) a zigzagging crack path is obtained, but it was impossible to determine why such a crack path was formed, although this phenomenon is often experimentally observed.

Slant fracture

For the plane strain specimen, the contours of total porosity at total fracture for different characteristic lengths are displayed in fig. 3.14. For l κ ≥ 40µm flat fracture is always obtained so that only one crack path is shown (l ω = l κ = 40µm). For l κ = 20µm slant fracture is always obtained (V-shape). In that case, the effect of l ω is small. The band width only slightly increases with increasing l ω . Figure 3.14: Contours of total porosity at total fracture for different characteristic lengths l κ and l ω for the plane strain case. Slant fracture is observed only for l κ = 20µm. The effect of l ω appears to be limited in this case.

Size effect

It is interesting to revisit the above observations from a more physical point of view and to consider, for a given set of material lengths, how specimen size impacts the formation of the cup-cone and slant fracture. Indeed increasing the material lengths for given specimen size is equivalent to decreasing the size of the specimen for given material lengths. Similar crack paths are obtained as long as the ratios l ω /L and l κ /L are constant where L is a characteristic dimension of the simulated structure (diameter or thickness in the present case).

From the results corresponding to the axisymmetric specimen (fig. 3.13), it can be concluded that it is easier to form a cup-cone crack path for large specimens (case corresponding to l κ < 80µm for a given geometry) than for smaller specimen sizes. From the results corresponding to the plane strain specimen (fig. 3.14), it can be concluded that only sufficiently large plane strain specimens break following the slant crack path. This opens the way to the determination of material lengths by using homothetic specimens so as to obtain different crack paths. Lengths could be fitted to represent the different fracture pattern. This will obviously require the use of very small specimens in the case of metals where lengths are expected to be of the order of the mean spacing between damage initiation sites.

Conclusions

In this work, a nonlocal GTN model for ductile fracture was proposed which uses two material lengths. The model is based on an implicit gradient formulation applied to the plastic volume variation (material length l ω ) and the accumulated plastic strain (material length l κ ). The model allows to regularize void growth and strain-controlled nucleation. The model parameters are fitted to reproduce the global response for tests carried out on a pipe line steel using a tensile bar and a plane strain specimen. These experiments have been already used in the literature to test several models and numerical techniques. Each test is representative of a characteristic crack path: cup-cone fracture (tensile test) and slant fracture (plane strain test).

The material parameters were first fitted assuming that both lengths are equal (l ω = l κ ) in order to reproduce the macroscopic behavior of both specimens. It was checked that mesh independence is obtained if a sufficiently fine mesh is used. In particular, it was shown that a local model is strongly dependent on mesh design (mesh orientation in the present study) whereas the nonlocal model is not. This implies that remeshing cannot be applied if a local model is used but with a nonlocal model, it can be a solution to reduce computational cost in cases where long crack propagation needs to be simulated.

A novel methodology was proposed to estimate, for a given geometry, the band width associated with a given characteristic length free from discretization error. By applying this methodology for different characteristic lengths, it was possible to establish a linear relation between the estimated band width and the characteristic length. This piece of information was then used to select the proper discretization of the localization bands for given values of the characteristic lengths in order to obtain converged solutions. In practice, the element size at failure is recommended to be three times smaller than the characteristic lengths.

Finally, two distinct characteristic lengths were used to simulate both test cases. It was first concluded that the characteristic length controlling strain controlled damage nucleation plays an important role on the occurrence of crack bifurcation: it must be small enough relative to the specimen size in order to observe a cup-cone crack path or a slanted crack path. It was also noted that increasing the characteristic lengths leads to an increase of both the band width and the radius at which the cup-cone crack path starts to bifurcate and propagate either at +45 • or -45 • . Larger values favor flat crack advance. This study also enabled to conclude that crack path should be affected by the size of the specimens. Pure flat fracture should be favored when using (very) small specimens. It could theoretically be possible to identify them separately. This paves the way to the more direct determination of material lengths by using homothetic specimens to obtain different crack paths. The fitted lengths should then be able to represent both the macroscopic behavior and the crack path for all sample sizes. Résumé en franc ¸ais L'approche continue-discontinue est présentée avec une stratégie pragmatique de remaillage pour optimiser le maillage et donc le temps de calcul. Dans un premier temps, de nouvelles composantes de la simulation de la propagation de fissures ductiles à l'aide de maillages 3D sont développées. En particulier, un nouvel algorithme d'initiation de fissures est proposé afin d'initier des fissures de forme correcte. En outre, un nouveau critère pour la transition continu-discontinu est développé afin de déterminer la taille de l'incrément de fissure à insérer dans le maillage en utilisant le remaillage. Le problème de la détérioration du taux de convergence après le remaillage est abordé, et une méthode de rééquilibrage est ensuite proposée pour surmonter ce problème. L'efficacité de l'approche continue-discontinue développée est testée par rapport à la discrétisation spatiale et à la discrétisation temporelle.

Enfin, des applications de la méthode développée utilisant différents maillages 2D et 3D de spécimens utilisés dans des tests expérimentaux sont présentées.

Taken from: A. El Ouazani Tuhami, S. Feld-Payet, S. Quilici, N. Osipov, J. Besson. New contributions for robust and cost-efficient continuous-discontinuous transition: application for complex crack path simulations., in preparation 1 Introduction

To predict crack initiation or when linear fracture mechanics cannot predict complex crack propagation, it can be interesting to use continuum damage mechanics models describing local material degradation phenomena. However, the continuum damage mechanics models have their limits. On the one hand, they cannot provide an accurate estimation of the opening and growth kinetics of the crack [START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF]. On the other hand, In the case of industrial problems with a large computational size, the propagation length is limited by the computational size of the continuous model due the fine discretization [START_REF] Chen | Modeling of ductile fracture using local approach : reliable simulation of crack extension[END_REF]. Indeed, the use of regularized methods [START_REF] Bazant | Non local continuum damage. localization, instability and convergence[END_REF]; [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF]; [START_REF] Lorentz | Ill-posed boundary conditions encountered in 3d and plate finite element simulations[END_REF]; [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]; [START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF] is primordial in case of ductile rupture, in order to avoid pathological mesh dependency Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF]; DE [START_REF] Borst | Fundamental issues in finite element analyses of localization of deformation[END_REF], and these types of regularized methods necessitate a fine refinement of the mesh inside localization bands. To overcome these two problems, different propositions exist in the literature for continuous-discontinuous models to fracture. They consist in inserting a strong discontinuity in the FE model based on the material degradation field. In the work of [START_REF] El | An anisotropic mesh adaptation strategy for damage and failure in ductile materials[END_REF], an anisotropic mesh adaption strategy was used to describe ductile fracture using bubble P1+/P1 tetrahedral elements. An implicit gradient formulation was used to overcome the mesh dependency. To drive the remeshing procedure, an error estimator based on the damage and damage rate was used. Moreover, to represent the crack path and avoid stretched elements in the vicinity of the crack tip, simple element erosion was used. A cup-cone simulation using a 3D mesh was presented to showcase their proposed remeshing strategy. However, the obtained meshes were refined over the whole crack path, which does not solve the issue of large computational size. In addition, the crack path resulting from a simple element erosion was indented, and the bandwidth of the removed element was not controlled. In the work of [START_REF] Javani | Three-dimensional finite element modeling of ductile crack initiation and propagation[END_REF], an implicit gradient method is also used to regularize the ill-posed mechanical problem. A global remeshing procedure was used to avoid element distortion occurring at high strains. An extension to 3D simulation of the crack path tracking algorithms and the discrete crack representation introduced in [START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF] was detailed. Also, a method capable of initiating 3D crack in the volume or near the body surface was also developed. Then, simulations in 3D showcasing simple crack paths were presented. However, their proposed method using the smoothed local maximum search over different radii to track crack paths is not capable of detecting crack branching. In the work of [START_REF] Yang | A 3d h-adaptive local remeshing technique for simulating the initiation and propagation of cracks in ductile materials[END_REF], a local h-adaptive remeshing procedure was developed. An improved transfer operator based on point selection was developed in order to optimize remeshing computational cost. However, the presented simulations refined the whole crack path; thus, the meshing procedure is not efficient for industrial simulations. Also, with the use of simple element deletion, the crack path is imprecise due to the indented crack path and to uncontrolled bandwidth of deleted elements. In the work of [START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF], an implicit gradient method was applied with multiple internal lengths to a multi-surface model with internal necking and damage mechanisms under shear. For damage to crack transition, cohesive bands were inserted between the element interfaces. The position and insertion time of the cohesive band was determined using a bifurcation analysis. Then, 2D simulations of cup-cone and slant fracture were presented. However, as the cohesive band was introduced between element interfaces, the resulting crack path was mesh-dependent. In addition, fine meshing is needed along the whole crack path; thus, the extension to real industrial simulation is limited due to the high computational cost. Another type of approach based on the XFEM method, and that does not necessitate fine meshing, was presented in [START_REF] Nikolakopoulos | Progressive failure of ductile metals: Description via a three-dimensional coupled czm-xfem based approach[END_REF]. A simulation of the cup-cone crack path was showcased using a coarse mesh. Thus this approach is not limited by the high computational cost. However, the framework was restricted to the small strain hypothesis due to the complexity of the proposed formalism to simulate ductile failure. In a more recent work, Eldahshan et al.

[2021a] proposed a phase-field approach coupled with remeshing and discrete crack insertion. The main contribution of this work was the development of a crack path tracking algorithm and crack insertion method based on remeshing adapted for parallel computations. Various numerical applications were presented to illustrates the capacities of the developed algorithm (CIPFAR: Crack Insertion and Propagation using the Phase Field and Adaptive Remeshing) to simulate ductile failure. However, the proposed continuous-discontinuous model uses a "bubble finite element" formulation (P1+/P1) based on the small strain hypothesis, which is not adapted for ductile failure even with the use of an updated Lagrangian formalism and remeshing. Also, the proposed remeshing strategy (detailed in Eldahshan et al. [2021b]) consists in refining the mesh in the process zone and does not offer the efficiency brought by mesh coarsening in zones already crossed by the discrete crack.

In this chapter, different novel contributions are proposed to arrive at a continuousdiscontinuous strategy enabling efficient simulations of complex crack path advance at a reasonable cost. The complete strategy is described in section 2, with more details on the new contributions. Then, the gain in robustness and cost efficiency offered by the proposed contributions is evaluated in section 3. Finally, the ability of the proposed strategy to capture complex crack paths, such as cup-cone and slant failure, or deal with challenging simulations, is illustrated in section 4 on several 2D and 3D cases. All simulations are performed with the Z-set software [START_REF] Foerch | Polymorphic constitutive equations in finite element codes[END_REF]; [START_REF] Besson | Large scale object-oriented finite element code design[END_REF]. The implicit gradient GTN model with two characteristic lengths of El Ouazani Tuhami et al. [accepted], presented in the previous chapter section 2 is used within a finite strain framework. A fully implicit resolution scheme is used.

continuous-discontinuous damage transition strategy 2.1 Regarding the constitutive behavior

Although this chapter proposes numerical contributions for the continuous-discontinuous transition, it is necessary, in order to illustrate the principles and evaluate the robustness of the proposed numerical tools, to consider real application cases. As our main applications cases deal with failure of ductile media, the implicit gradient GTN model with two characteristic lengths of El Ouazani Tuhami et al. [accepted] is used. This model is able to describe all three stages of ductile failure within the finite strain framework and is especially suited to model complex crack path such as the one associated with cup cone failure. It is therefore perfectly suited to serve for this chapter. However, the proposed contributions could be used, with slight adjustments, in combination with different models. That is why the considered constitutive behavior is only briefly recalled in the previous chapter in section 2. The only requirements for the model are:

• There should be a continuous scalar field, f, D, ..., representing material degradation;

• It is possible to estimate a critical value f crit for f corresponding to complete material degradation (or in other words, a material point with f crit = f r can be considered as broken);

• Material degradation is more or less localized in a band of width l b .

In the present case, the best scalar field candidate for the considered nonlocal GTN model is the total porosity variable f t .

Summary of the strategy

The chosen continuous-discontinuous damage transition strategy can be summarized by the following choices:

1. Mesh adaption to produce optimal meshing: (i) before crack initiation, the map size of the mesh is set by the user and actualized after a chosen time increment ∆t remeshing .

(ii) after crack initiation, an automatic mesh refinement of the active process zone is carried out, which guaranties mesh coarsening in already cracked zones;

2. Discontinuity representation method using mesh adaption by representing the crack as a strong discontinuity in the mesh;

3. Strategy to define where to insert a crack increment: marching ridges algorithm for propagation in 2D and 3D and 2D crack initiation and new strategy for 3D crack initiation;

4. Criterion to determine when to insert a crack increment: new criterion based on the crack increment's surface;

5. Procedure to retrieve equilibrium after the model has been changed due to crack increment insertion: transfer with finite element interpolation for nodal data, and nearest integration point association for integration points data and new reequilibrium procedure consisting in a relaxation stage using the material behavior for a few increments as opposed to the elastic step reequilibrium where only the elastic behavior is considered.

Each component of the proposed procedure is detailed in the following. The sequence of the corresponding operations is represented in figure 4.1.

Figure 4.1: Flowchart describing the proposed continuous-discontinuous approach with the different involved components and their inputs/outputs. After resolution of the mechanical problem, the considered scalar field related to material degradation f is smoothed thanks to the superconvergent patch recovery method. Then a crack path tracking algorithm is used to determine where the next crack increment should be inserted and its geometry. The crack increment is not inserted until the insertion criterion is verified. Only then is a new mesh generated thanks to the mesh generator meshGems based, on the one hand, on the intersection of the structure's mesh and the crack increment mesh, and on the other hand on the tracking of the active process zones which are meant to be refined. The variables are then transferred onto the new mesh. To facilitate the return to equilibrium and ensure a satisfactory convergence rate, a reequilibrium step is proposed. Computation of the next load increment can then be performed.

2.3

Where to insert a crack increment?

Crack path tracking algorithms

To define the location of a continuous and sufficiently regular crack surface, crack path tracking algorithms considering an underlying spatial scalar field representing material degradation are preferred here. Indeed, these algorithms use only local information near the crack front, thus the computational cost is negligible compared to the FEM computation.

This low computational cost is essential for a procedure like the one proposed as it is called after each converged load increment. Some rarely addressed challenges for these algorithms are crack initiation completely inside the volume of a structure and crack branching. The ability of the algorithms to overcome these challenges can serve as differentiation points. In this regard, the marching ridges algorithm introduced in Feld-Payet et al. [2015] enables to deal with both situations and is thus chosen in this work for the crack propagation phase. However, although this algorithm provides a unified way to deal with crack propagation and crack initiation, a more cost effective solution which leads to a more realistic first crack increment is proposed for crack initiation in 3D. In both cases, the algorithm provides an explicit surface of the crack increment that can be used directly by the cutting algorithm [START_REF] Chiaruttini | Advanced remeshing techniques for complex 3d crack propagation[END_REF] to introduce the crack into the computational mesh. Let us start with the initiation of the crack in the computational mesh which leads to the creation of a first crack front. The crack propagation, modeled by the advancement of the crack front is then tackled.

Crack initiation

(i) In 3D The geometrical approach proposed in this chapter consist in approximating the first crack increment for 3D cases by an elliptical surface. This surface is defined as the ellipse shape containing the orthogonal projection of the positions of the broken integration points (i.e. points having reached the critical value f crit ), see figure 4.2.

Figure 4.2: Ellipse shape representing the first crack increment build from a cluster of points having reached the critical value f crit . The ellipse (in red) is defined, from its center G corresponding to the barycenter of the broken points, so that it contains all the projections on its surface of the considered points (in blue).

In practice, the center of the ellipse is classically chosen to be the barycenter of the broken points. The plane of the elliptical shape is obtained thanks to the Singular Value Decomposition (SVD) theory [START_REF] Stewart | On the early history of the singular value decomposition[END_REF]. According to this approach, if the coordinates of the considered points are stored in a matrix M coord ∈ R 3×N broken , where N broken is the number of broken points, then this matrix can be decomposed as follows:

M coord = V ΣW t (4.1)
where V ∈ R 3×3 and W ∈ R N broken ×N broken are orthogonal matrices and Σ is the matrix with diagonal coefficients corresponding to the eigenvalue in ascending order (so that it is uniquely defined). Then, the columns of the matrix V , respectively V 1 , V 2 and V 3 , represent the preponderant vectors describing the cluster of broken points. In particular, the plane of the ellipse is defined by

(G, V 1 , V 2 ), as illustrated in figure 4.2.
Finally, the radii of the ellipse, denoted a and b, are choses so that the ellipse contains all the projections of the broken points on the ellipse plane. To do so, it is useful to consider the matrix M * coord = ΣW t which contains the coordinates of each broken point in the new coordinate system

(G, V 1 , V 2 , V 3 )).
The radii of the ellipse can then simply be defined by the following:

• for the largest radius:

a = max i (|M * coord [0, i]|) for i ∈ [0, N broken [;
• and for the smallest radius: ) In 2D For 2D cases, the marching ridges algorithm provides an efficient and sufficiently cost efficient approach to determine the first crack increment. Its principle is recalled in the next section.

b = max i (|M * coord [1, i]|) for i ∈ [0, N broken [. (ii

Crack propagation using the marching ridges algorithm

(i) Principle in 2D: In 2D, the idea of the marching ridges algorithm is to start from a point of the crack and approximate the next crack increment by a line segment. As, the crack is assumed to appear where material degradation is maximum, the segment should be placed on the ridge of the considered scalar field representing material degradation f . As with any crack path tracking algorithm considering a scalar field, the first step is a smoothing operation performed on the considered scalar field f . Following Feld-Payet et al.

[2015], this is done using the super-convergent patch recovery (SPR) technique [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF]. The smoothed field is then defined using the same interpolation functions as for the displacements. Consequently, the evaluation of the smoothed degradation variable f SP R or its spatial gradient at any point is straightforward. Then, it is necessary to define a starting point x 0 : in 2D, a simple choice is to consider the global maximum of the scalar field for crack initiation and the current crack front for crack propagation. This point can then be used as the center of a polar grid used to evaluate the projection of the spatial gradient of the considered field f SP R onto the tangential vector e θ at an evaluation radius R and for a finite number of angles θ separated by the angular accuracy ∆θ. The next points on the ridge then correspond to the points (R, θ) for which the derivative of f SP R with respect to θ is positive before (i.e for θ -∆θ/2) and negative after (i.e for θ + ∆θ/2):

e θ (θ -∆θ/2).∇f SP R (x 0 + R e r (θ -∆θ/2)) > 0 e θ (θ + ∆θ/2).∇f SP R (x 0 + R e r (θ + ∆θ/2)) < 0 (4.2)
(ii) Adaptation to 3D for crack propagation In Feld-Payet et al. [2015], a simple method to adapt this planar search for 3D crack propagation was proposed. It consists in performing the 2D search in planes that are orthogonal to the current crack front. In practice, the current crack front is a line that can be approximated by a finite number of segments. Let us note that this discretization can be completely independent of the FE discretization: e.g. in this work, the crack front line is discretized into segments of length equal to the evaluation radius R. For each segment, it is straightforward to define an orthogonal plane for the 2D search. For increased robustness, several parallel planes can be defined for one segment so as to be able to average the resulting direction(s), especially for large segments. The starting point x 0 is then the intersection of the considered plane with the current crack front line. Finally, an average direction is computed for each point connecting two segments using their respective directions (see figure 4.3). This 3D adaption method is actually quite generic and can be used with any crack path tracking algorithm developed for 2D cases.

A major difference in this work compared to Feld-Payet et al. [2015], is that the crack increments do not have the same length for all the segments discretizing the current crack front in 3D, which broadens the scope of the methodology. Indeed, for each point between two segments, the propagation length is given by the distance to the iso-value f = f crit . Each point of the current front is then associated to a new end point that serves as control point to build a Bezier curve of degree 2 . Finally, the next crack front is defined by projecting each node of the current crack front onto the Bezier curve, which completely defines the crack increment.

As material degradation progresses, the iso-line f = f crit usually gets further from the crack front, thus increasing the size of the crack increment until the insertion criterion allows its insertion.

Insertion criterion

An insertion criterion determines whether a crack increment is ready to be inserted into the mesh based on the local fields. This insertion should be performed when material degradation is sufficiently important.

In 2D

In case of 2D simulations, the insertion criterion can be easily defined since the crack increment is a line segment (for both initiation and propagation), especially when inserting crack increments of constant length. It then suffices to evaluate the considered material degradation related field f and to allow crack insertion when it exceeds the critical value f crit along the entire crack increment, as in [START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF]. This is also the choice made here for 2D problems. The main challenge is then to define an appropriate crack increment size to avoid too frequent remeshing if the crack increment is too small, or delayed insertion if the crack increment is too large, resulting in stretched elements near the crack front. Following [START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF], a crack increment length ∆a build as a multiple of the localization bandwidth l b : ∆a = m c × l b , with m c is a multiplication factor that can be tuned (a study on the effect of the crack increment length ∆a on the crack path is presented in section 3.1.2 leading to m c ≈ 4).

A more elaborate criterion for 3D structures is needed since the crack can be represented by a complex surface.

For crack initiation in 3D

For crack initiation, a minimal size of the elliptical crack surface should be determined in order to avoid inserting too small increments into the mesh. In order to be as consistent as possible with the criterion proposed for 2D cases, it is proposed here to proceed to insertion only if both radii of the ellipse (build to contain all the projections of the points with f ≥ f crit ) are greater than ∆a.

For crack propagation in 3D

When the crack increment does not have a simple geometrical shape, as it may be the case during crack propagation, a new criterion has to be defined. We propose here to consider the area of the crack increment, A real , and to compare it with the area considered when using constant crack length increments, which had provided good results in [START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF]. A progression with constant crack increments corresponds to an isotropic propagation of the current crack front. The area of the uniform crack increment is thus denoted A unif . The real crack increment is validated for insertion if its area A real is greater than the area of the uniform crack propagation A unif (see figure 4.4).

Mesh adaption 2.5.1 Discontinuity representation

In order to obtain a realistic representation of cracks' opening/closing, a discrete representation of crack lips in the computation mesh is preferred. In the present case, the discrete representation is obtained thanks to mesh adaption since this operation is already used even before crack initiation in order to avoid element distorsion.

Insertion of the discrete crack into the volumetric mesh here relies on a cutting operation of the external boundary of the mesh by the crack surface. A robust cutting algorithm introduced in [START_REF] Chiaruttini | Advanced remeshing techniques for complex 3d crack propagation[END_REF] is used to perform this operation. Note that the same cutting algorithm is used to perform discrete crack propagation in fatigue problems [START_REF] Proudhon | 3d simulation of short fatigue crack propagation by finite element crystal plasticity and remeshing[END_REF]; [START_REF] Fessler | Relation between crack growth behaviour and crack front morphology under hold-time conditions in da inconel 718[END_REF]. This cutting algorithm produces a surface mesh incorporating both the outer surface of the structure mesh and the crack surface. Then, an additional refinement operation, using meshGems [START_REF] Frey | Mesh generation: Application to finite elements: Second edition[END_REF], is performed in order to fill in the volume of the created surface mesh and generate a suitable mesh with acceptable element quality and adapted discretization. This refinement operation is detailed in the next section.

Local refinement of active localization bands

Through refining specific areas and coarsening the mesh in others, the objective is to achieve a compromise between accuracy and cost. This strategy requires to distinguish Active Damage Zones (APZ) from areas where damage does not significantly evolve to only refine the APZ. Inside the APZ, the minimal element size h min necessary for the convergence of the mesh is imposed(see El Ouazani Tuhami et al. [accepted]). Outside the APZ, the sizes of the elements are linearly interpolated between the minimal and maximal sizes, respectively h min and h max , according to the distance of the considered element from the APZ; this enables to minimize the computational cost.

In practice, the APZ is defined differently according the simulation stage:

• Before crack initiation, the APZ is defined by a circle in 2D, or a sphere in 3D, centered on the potential locations of crack initiation. The radius of the circle L AP Z is defined as the maximum distance between two points reaching f crit and f crit /2 respectively. Depending on the size of the APZ, this radius can be adapted by considering a larger area through a more general definition of the radius involving f crit /m AP Z with m AP Z ≥ 2. This choice is made so that an adequate fine mesh is attributed to zones of crack initiation at early stages of the damaging process (see figure 4.5-a). Note that this procedure is based on prior knowledge of the potential locations of crack initiation. However, adaptive mesh procedures based on error estimation (see chapter 5, Feld-Payet [2010]) may automate the mesh adaption in this stage of the simulation, but fine tuning of the error estimator parameters and the map size computation are needed to achieve the same goal;

• After crack initiation: the same definition is used, but the center of the APZ changes. In 2D, the center is considered as the farthest point from the crack front, taken along the last crack increment direction, where the considered scalar field f reaches f crit /m AP Z .

Then the distance between this point and the crack front is the radius L AP Z (see figure 4.5-b). In 3D, the center line of the APZ is constructed as the 1D line intersecting the points reaching the f crit /m AP Z in front of the crack front (also a line in 3D). Thus the APZ is the cylinder revolving around the center line and its radii is define by the distance between the center line and the crack front. Note that in the propagation phase, as the damage already evolved and reached the critical value f crit , the APZ refinement is done automatically as opposed to the initiation stage. Using this procedure, only the active process zone ahead of the crack front is finely meshed. The mesh used for areas already crossed by a crack are coarser depending on their distance to the crack front.

Mesh adaption within a finite strain framework

Since the FE simulation is performed within a finite strain framework, it is important to ensure that the aspect ratios of the elements stay correct throughout the simulation. To do so, mesh adaption is performed on the current configuration and at a given frequency ∆t remeshing set by the user. This first criterion is considered only before crack initiation.

Transfer operator

In order to resume the simulation after each discretization change, it is necessary to transfer, from the old mesh to the new mesh, the data stored at the nodes and integration points.

Nodal data transfer operator

The procedure for transferring the nodal variables is standard: each node of the new mesh is located in an element of the old mesh, and then the new nodal value is calculated using the element's shape functions. In the present case, since the remeshing is performed on the current configuration, the displacement vector u is set to 0 on the new mesh, so no transfer is needed for the displacement field. However, since the considered model is based on an implicit gradient formulation, the nonlocal nodal variable are transferred using the described method.

Integration points data transfer operator

For the state variables stored at the integration points, the nearest integration point transfer technique is used [START_REF] Simone | From continuous to discontinuous failure in a gradient-enhanced continuum damage model[END_REF]. This technique has been chosen to preserve the high gradients of the localized variables, since no smoothing operation is used. Let us note that a prerequisite to use this technique is that the meshes do not change significantly before and after remeshing. This is the case with the proposed mesh adaption procedure thanks to a regular refinement (due to the mesh adaption frequency ∆t remeshing or to the insertion of sufficiently small crack increments) of the APZ which moves progressively as the crack advances, with a gradually evolving mesh size outside the APZ.

Reequilbrium using a relaxation step

Independent transfer operations on the nodal and integration point variables create unavoidable errors: in particular, neither the global equilibrium equations nor the local constitutive equations are satisfied and the transferred strains are not consistent with the nodal displacements. These inconsistencies in the transferred fields degrade the convergence rate after remeshing and it may lead to a large number of increments if a load increment division procedure [START_REF] Boroomand | Recovery procedures in error estimation and adaptivity. part ii: Adaptivity in nonlinear problems of elasto-plasticity behaviour[END_REF] is used. Too significant errors may even lead to divergence of the resolution algorithm.

In this chapter, the author propose a novel reequilibrium method based on the elastic balancing step introduced in [START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF]. The idea of this latter method is to retrieve equilibrium thanks to an additional step with frozen boundary conditions during which there should be no evolution of state variables like plasticity and damage. This balancing step then consists in solving the mechanical problem using the transferred nodal and integration data, while considering an elastic behavior. However in presence of damage and crack insertion, even with this balancing step, from the author experience, degradation of convergence rate and risk of divergence are still a problem. Part of the problem may come from the fact that plasticity and damage variables are kept constant.

For the proposed reequilibrium procedure, only adapted for viscous material, the boundary conditions are also frozen but during a non null time increment ∆t reeq during which the mechanical problem is solved, with several sub-increments, considering the material behavior as opposed to the elastic behavior used in [START_REF] Varas | Continuous and discontinuous modelling of ductile fracture[END_REF]. This procedure can be interpreted as a viscous relaxation of the material. Plasticity and damage are then free to evolve, but this evolution should actually be rather negligible due to the fixed boundary conditions and the choice of a small time increment for the reequilibrium procedure. Indeed, the reequilibrium time increment ∆t reeq is taken as a fraction (1/10, 1/100 or 1/1000) of the time increment required for convergence during the ductile failure process. Note that, the closer ∆ reeq to the time increment required for convergence, the better the rate convergence, but in the simulations presented in this chapter, a smaller value of this parameter is chosen to limit the state fields variation during this reequilibrium step. This time increment is divided onto sub-increments which are taken very small at first and they rapidly increases, until the imposed time ∆t reeq has elapsed. In average, we have observed, for the considered model, a number of increments during the reequilibrium stage following a crack increment insertion around 8. After this reequilibrium increment, the usual boundary conditions are considered again and the usual convergence rate should be retrieved, or only slightly diminished.

Summary of the parameters of the continuous-discontinuous procedure

In table 4.1, the parameters of the the continuous-discontinuous procedure used for all the applications in this chapter are listed.

Evaluation of robustness and cost efficiency

In this section, first, the robustness of the considered continuous-discontinuous approach, and in particular of the proposed new ingredients to build the crack surface is illustrated on 2D nonlocal simulations of an axisymmetric tensile test for different mesh sizes and different crack increment sizes ∆a. Then the gain in cost offered by the use of the mesh adaption based on the tracking of the APZ and the reequilibrium procedure is illustrated on the same application. In order to illustrate the robustness of the continuous-discontinuous approach, and in particular of the proposed new ingredients to build the crack surface, let us compare the results obtained with different mesh sizes for the same structure. Let us consider the axisymmetric round tensile specimen (see figure 3.1a) corresponding to one of the experimental tests in [START_REF] Rivalin | Développement d'aciers pour gazoducs à haute limite d'élasticité et ténacité élevée : Mécanique et mécanismes de la rupture ductile à grande vitesse[END_REF], the same experimental test was used in literature as benchmark simulation [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF]; [START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF]; [START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF]; El Ouazani Tuhami et al. [accepted].

In order to obtain converged results, the localization band should be discretized by several elements. It is thus necessary to first determine the width of the localization band l b . In the present case, it is possible to refer to El Ouazani Tuhami et al. [accepted] Consequently, meshes with with minimum size h min ∈ {10, 20, 40}µm are considered. For comparison purposes, a simulation with a very fine fixed mesh from with N h = 180 elements in the minimal cross section and an initial element aspect ratio equal to 6 : 1, without crack insertion, is also considered. In Figure 4.6, the evolution of the engineering stress is displayed as a function of the diameter reduction for all the considered meshes and the total porosity f t for each mesh size at total failure is also shown. For meshes with minimal element size h min ≤ 40µm ≤ l b , the global response is very similar to the simulation using a very fine mesh (N h = 180) without crack insertion, which is assumed to provide converged results. A cup-cone crack path is obtained for all the considered mesh. It is clear that the crack lip opening is more realistic using the proposed procedure compared to the simulation without crack insertion. In the last case, the crack is represented in a smeared manner, thus the position of the crack lips or crack front is simply not defined.

Crack increment size effect

The crack increment length ∆a is an important parameter of the proposed procedure because it affects both the discretization of the crack path and the number of mesh adaptions (through the insertion criterion). On the one hand, its choice is guided by considerations similar to the ones guiding the choice of the elements' size for any finite element computation, i.e. it should be :

1. small enough to enable a proper discretization, in this case, of the crack path; With minimal mesh sizes h min ∈ {10, 20, 40}µm, the proposed continuous-discontinuous method has been used whereas the simulation with a very fine mesh is without any mesh adaption or crack insertion. For all the considered meshes, the global responses are similar and the total porosity at complete failure shows a cup-cone failure. This similarity indicates that the proposed procedure is rather robust to the choice of the mesh size.

2. large enough in order to minimize computational cost, in particular to limit the cost associated to a change of discretization by limiting the number of crack insertions.

Let us note that in the present case, choosing large increment also means less transfer related errors, whereas for standard finite element computation, larger elements usually lead to more discretization error.

On the other hand, this choice should also be guided by an additional constraint when considering a scalar variable presenting a plateau instead of a sharp ridge: the crack increment length should not be smaller than plateau width, otherwise, the crack path may zigzag. This is the case here since the porosity variable ends up saturating locally due to the choice of regularization through an implicit gradient formulation. That is why only crack increment length larger than the localization bandwidth are considered in this study. Let us then compare 3 simulations with the same minimum element size h min = 20µm but different crack increment lengths: ∆a ∈ {80, 160, 240}µm. As shown in figure 4.7, the evolution of the engineering stress as a function of the diameter reduction is the same for the 3 crack increment lengths. It can thus be concluded that the global response is rather robust to the choice of the crack increment length.

This also means that it is necessary to use different criteria to choose an appropriate crack increment length. By looking at the different contours of total porosity at complete failure in figure 4.7, one can see that for the largest increment, i.e. ∆a = 240µm, the crack path is not completely smooth, which can be explained by excessive stretching of the elements near the crack front due to delayed insertion. Smaller sizes should thus be considered. The choices ∆a = 160µm and ∆a = 80µm lead to rather similar crack paths. However, in the first case, only 20 remeshes are necessary to obtain the complete propagation in comparison to the 38 remeshes for the smallest crack increment length. In order to minimize computational costs, the value ∆a = 160µm should be preferred. This crack increment length corresponds to four times the evaluated localization bandwidth (i.e. ∆a = 4l b ), which is the order of magnitude recommanded by the author. -∆φ 0 /φ for different crack increment lengths ∆a, using the proposed continuousdiscontinuous method, in the case of the axisymmetric specimen. It appears that ∆a has no effect on the global response. However, the largest increment length, ∆a = 180µm, leads to a cup-cone crack path which is not completely smooth compared to the other crack increment lengths.

Cost efficiency

Gain from refining only the active process zones

In this section, the gain obtained thanks to the refinement of only the active process zones is evaluated. To do so, the three following approaches are compared in terms of the number of nodes:

1. the proposed continuous-discontinuous procedure with refinement of only the active process zones;

2. an approach with mesh refinement of the all the process zones without coarsening after a process zone is no longer active;

3. a simulation with a fixed mesh (i.e. without any mesh adaption) of rectangular quadratic elements with 239 123 nodes. This last case is meant to be representative of a situation where the exact crack path is not known a priori, so minimum size elements to achieve convergence are used in the entire zone of possible crack propagation.

For both cases with mesh adaption, triangular elements are used with a minimum element size h min = 20µm. Let us specify that all three approaches lead to similar global responses and identical crack paths. In Figure 4.8, the evolution of the total number of nodes as the meshes are adapted is given for all the three approaches. Let us note that the smallest number of nodes is obtained with the proposed procedure: about 3 times less than the mesh adaption procedure without proposed approach refinement without coarsening fixed mesh (rectangular element) Figure 4.8: Evolution of the number of nodes as a function of the number of mesh adaption procedures for the proposed procedure with refinement of only the active process zones (in blue), mesh adaption with only refinement without coarsening of the process zones (in orange) and a fixed mesh (in green). The firt approache consistently leads to a smaller and rather constant number of elements. coarsening and 13 times less than the fixed mesh. Besides, with the proposed approach, the total number of nodes is approximately the same for all remeshing instances. In other words, the total number of nodes does not depend on the length of the crack path, in contrast with the second procedure without mesh coarsening for which the number of nodes increases as the crack propagates. The proposed approach consisting in refining only the active process zones is thus very cost effective.

Gain from the reequilibrium procedure

In this section, in order to evaluate the cost effectiveness of the reequilibrium procedure, let us compare the computational cost for three simulations:

1. a simulation with the full proposed approach and in particular, the reequilibrium procedure;

2. a simulation with mesh adaption only in the active process zone, but without the reequilibrium procedure;

3. a fixed mesh similar to the one presented in the previous section 3.2.1.

Since each iteration involves the update of the local internal variables and the inversion of the global matrix, which are quite expensive operations, especially with a non linear problem, the accumulated number of iterations is an appropriate measure of the computational cost.

In figure 4.9, the evolution of the accumulated number of increments is given as a function of the imposed vertical displacement U d 2 , from U d 2 = 3mm. Let us specify that the global response obtained in all the three considered cases is identical. The lowest number of increments is obtained with the fixed mesh (i.e. without any field transfer): only 690time increments are needed to obtain total failure. But even with this fixed mesh, there is an increase of the number of increments after crack initiation (first integration points reaching f crit ). This trend is also observed with the simulations with mesh adaption. It can thus be assumed that, whatever the considered approach is, the severity of the nonlinearities associated with the damage evolution during the fracture process may be responsible for a part of the difficulties to converge. However, additional difficulties appear when adapting the mesh to insert a crack increment, as the number of increments seems to increase as soon as the crack initiates (compared to the relatively small number of increments before). There is indeed a cost that has to be paid in order to introduce, in the continuous model, a discrete and well defined crack. However, this cost can be limited thanks to the proposed reequilibrium procedure. Indeed, the number of time increments is two times smaller than without reequilibrium. .9: Evolution of the cumulated number of time increments as a function of imposed vertical displacement for the simulation of the axisymmetric tensile test on a fixed mesh (in black), with mesh adaption but without reequilibrium (in red) and with mesh adaption and the proposed reequilibrium procedure (in dark red). The higher number of time increments is obtained for the simulation with mesh adaption without reequilibrium. With the reequilibrium procedure, the total number of increment is approximately divided by two. The lowest number of time increments is obtained for the fixed mesh for which there is no discrete crack represented, so no transfer needed. In all cases, the crack propagation phase is accompanied by an increase in the number of increments.

4 Numerical applications

Introduction

The purpose of this section is to demonstrate the effectiveness of the proposed continuousdiscontinuous approach, and in particular of the new ingredients, on several cases from the literature. The considered cases are selected among the most used specimens for material identification and either lead to a non-trivial crack paths (like tensile tests on notched round or round bar specimens) or present some difficulty for its simulation (like a tensile test on a CT specimen resulting in relatively long crack propagation with very large strains). Numerical applications are presented in various frameworks: 2D axisymmetric, 2D plane strain and 3D simulations. The material characteristic presented in table 3.1 are used unless specified otherwise.

Cup-cone failure prediction for different 2D axisymmetric tensile tests

In this section, more cases presenting a cup-cone failure are tackled: round specimens with different notch radii ({2, 4, 10}mm) and without any notch (referred to as round bar and assimilable to a notched specimen with an infinite radius R n ) are simulated with 2D axisymmetric meshes: see figure 4.10-a for a description of the geometry and boundary conditions. The interest of considering several radii is to study the effect on the crack path of different stress triaxiality states at the center of the specimen, since triaxility level increases as the notch radius decreases. The results are first compared in terms of evolution of the engineering stress as a function of the diameter reduction -∆φ 0 /φ in figure 4.11. It appears .11: Evolution of the engineering stress as a function of diameter reduction -∆φ 0 /φ for all the considered axisymmetric specimens, using the proposed continuous-discontinuous approach. The final diameter reduction decreases with the notch radius. For each global response, three particular instants of the simulation are highlighted with a numbered arrow in order to track damage evolution and crack path in figure 4.12.

the crack is initiated (stage 1) at the center where stress triaxiality is maximum. Then for all the specimens a cup-cone crack path is predicted (stage 3). It is interesting to note that the distance from the symetry axis at which the cup-cone crack path bifurcates increases as the notch radius decreases: this indicates that the cup-cone bifurcation diameter is a function of the triaxiality at the center of the specimen.

Figure 4.12: Total porosity f t contours at the three instants of the simulation highlighted in figure 4.11, for all the considered axysimmetric specimens. For all the specimens, there is a cup-cone crack path at stage 3. The diameter corresponding to the bifurcation increases as the notch radius decreases. Let us note that, for the round bar specimen, striction causes the boundary to appear as curved at these advanced stages of material degradation.

Slant fracture prediction for 2D plane strain tensile tests

In this section, the proposed continuous-discontinuous procedure is applied to predict slant fracture for a tensile test selected from [START_REF] Rivalin | Développement d'aciers pour gazoducs à haute limite d'élasticité et ténacité élevée : Mécanique et mécanismes de la rupture ductile à grande vitesse[END_REF], the geometry and boundary conditions are presented in figure 4.10-b. It is interesting in this case to compare the results with a fixed mesh without crack insertion since, contrary to the cup-cone fracture simulations in figure 4.6, slight differences are found in this case. Let us note that both simulations are performed assuming plane strain conditions. Besides, both discretizations are chosen sufficiently fine (i.e. with a minimum element size of h min = 20µm for the triangular meshes generated with the proposed procedure and with N h = 180 rectangular elements in the section for the fixed mesh) in order to have converged results.

In figure 4.13, the evolution of the engineering stress is plotted as a function of the thickness reduction -∆e 0 /e 0 for both the simulation using the proposed procedure and the simulation with a fixed mesh. Similar global responses for both simulations are obtained up to the latest stages of the simulation: at this point, the simulation with crack insertion ends up with a slightly larger final thickness reduction. In order to better understand this difference, it is interesting to consider the total porosity contours in the same figure. In both cases, a slant crack path is predicted. However, with the fixed mesh, the slant crack path is symmetric with regards to the central vertical axis, whereas with the proposed continuousdiscontinuous procedure, an antisymmetric slant crack path first appears before a bifurcation towards another shear band is observed. The occurence of the crack bifurcation seems to coincide with the deviation of the global curves: the difference between the crack paths may then explain the difference in the global responses at the latest stages of the simulation.

Let us note that this bifurcation of the slant fracture towards another shear band was also reported in the numerical work of [START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] using a different ductile model and a different continuous-discontinuous approach. However, in the work of Leclerc et al., only one half of the specimen was modeled, which means that the symmetry conditions impose a symmetric crack path with regards to the central vertical axis. From these results, it can be assumed that, rather than the symmetric or antisymmetric characteristic of the crack path, it might be the insertion of a true discontinuity in the model that favors the occurence of crack bifurcation. However, more simulation with fixed meshes with different characteristic lengths would be necessary to support this assumption. .13: Evolution of the engineering stress as a function of the thickness reduction -∆e 0 /e 0 of the plane strain specimen for both a simulation with the proposed continuousdiscontinuous procedure and a simulation with a fixed mesh. The global curves are similar up to the latest stages of the simulation where the simulation with the proposed continuousdiscontinuous procedure ends up with a larger final thickness reduction. Total porosity is displayed at various times of the simulation using the proposed method, and at total failure for the fixed mesh simulation. The proposed method predicts an axisymmetric slant crack path with a bifurcation whereas the fixed mesh simulation predicts a symmetric slant crack path. The bifurcation coincide with the deviation of the global responses.

Blunting and long crack propagation for a 2D plane strain CT test

In this section, a 2D plane strain simulation of a tensile test on a CT specimen is used to showcase the capability of the proposed continuous-discontinuous procedure, and in particular of the mesh adaption procedure of only refinement of APZ, to predict crack tip blunting before crack initiation and to simulate long crack path propagation accompanied by very large strains. The geometry and boundary conditions can be found in figure 4.10c For this geometry, since no prior convergence study was performed to determine the maximal acceptable size for converged results according to the methodology proposed in El Ouazani Tuhami et al. [accepted], two meshes with different minimum element sizes h min ∈ {20, 40}µm are used to verify the mesh convergence. In figure 4.14, the evolution of the engineering stress is plotted as a function of Crack MOuth Displacement (CMOD) for the two considered meshes. Both curves are identical which tends to indicate that mesh convergence is achieved. Besides, based on the contours of total porosity in the same figure, it can be seen that mesh adaption enables to maintain a good quality for the elements near the crack tip, even in the presence of large strains. Indeed, before crack initiation, Before crack initiation (marked with a green cross), crack tip blunting is well described thanks to mesh adaption which enables to preserve good quality elements near the crack tip. Then the crack propagates over a relatively long distance with still well shaped fine elements near the crack tip and a coarsened mesh along the crack lips.

crack tip blunting is gradually developing and well captured by the simulation. This good representation of the crack tip blunting would be more difficult to obtain with a fixed mesh due to stretched element near the crack tip. During crack propagation, good quality elements are preserved in the vicinity of the crack tip while the non-active process zones behind the crack tip are coarsened. This enables to keep a the number of nodes throughout the crack propagation process approximately the same and thus to limit computational cost. Furthermore, thanks to crack insertion into the computation mesh, a stable and relatively long crack propagation until total failure is obtained with, all along, a well represented crack opening.

3D simulation of cup-cone failure

In this section, a 3D simulation of a cup-cone failure on a simple round bar (without any notch) is presented to illustrate the efficiency of the proposed ingredients that are specific to 3D, i.e. the proposed geometrical approach based on a scalar damage field to initiate a crack with an elliptical crack surface and the proposed insertion criterion based on the area of the crack increment. The considered initial mesh is presented in figure 4.16-a: only a quarter of the specimen is simulated. In spite of the gain offered by the new ingredients (see 3.2) for the post-processing operations related to crack increment evaluation and insertion, 3D simulations are still quite expensive during the increments' calculation due to integration of the constitutive behavior for a very large number of elements and the resolution of the linearized mechanical problem. In order to reduce this last part of computational costs which are not the subject of this chapter, larger elements are considered here, with h min = 80µm (i.e. four times larger than for the previously presented 2D simulations of cup-cone in 4.2). Consequently, larger minimum characteristic lengths are also used to have a sufficient number of elements to discretize the localization band (here l κ = l ω = 40µm, see previous chapter section 5.2). The results are confronted with the ones obtained performing, on the one hand, a 2D axisymmetric simulation using the proposed continuous-discontinuous procedure and a finer mesh (h min = 20µm) and on the other hand, a 2D axisymmetric simulation with a fixed mesh with N h = 180 rectangular elements discretizing the section. In figure 4.15, the evolution of the engineering stress is plotted as a function of the diameter reduction -∆φ 0 /φ for all the three considered simulations. The global responses from the three simulations are very similar, which tends to validate the results of the 3D simulation. Contours of the total porosity f t are also displayed at different times for the 3D simulation to illustrate the cup-cone crack path. The final cup-cone fracture surface obtained with the 3D simulation is shown in figure 4.16-b. It is clear that, thanks to the proposed continuous-discontinuous procedure, the crack lips are particularly well described and smooth. This is an advantage of discrete crack determination, which requires the use of a crack path tracking algorithm. Let us note that contrary to continuous-discontinuous procedures that use element erosion prior to mesh adaption, with the proposed procedure there is abolutely no material loss. Let us note that the cup-cone fracture surface obtained with the 3D simulation is in good qualitative agreement with experimental results of [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF] (reproduced in [START_REF] El | An anisotropic mesh adaptation strategy for damage and failure in ductile materials[END_REF]).

Conclusion

In this chapter, some new contributions were proposed to enable a robust and cost-efficient transition between a continuous model evaluating local material degradation and a model representing propagation of a discrete crack. These contributions have been evaluated on complex crack path simulations (i.e. either non-trivial crack paths or challenging simulations , a new procedure was proposed in order to facilitate the return to equilibrium and to improve the convergence rate after change of discretization. As a consequence, several simulations of specimens used for material identification with both plane strain and axisymmetric elements and different numerical parameters could be performed and presented in this chapter. It has been shown that this reequilibrium procedure enables to limit the number of increments during crack propagation, which contributes to the cost-efficiency of the method. Another new element, that also enables to limit computational cost, is the refinement of only the active process zones. This feature has a double advantage: not only are the degradation phenomena leading to complete failure is well described, but also, the number of elements (and thus the computational cost) remains approximately constant. Finally, two more contributions, on rarely tackled subjects, are proposed for simulations in 3D. The first one deals with crack initiation completely inside the structure: a simple geometrical crack initiation methodology based on SVD analysis of the coordinates of broken points has been proposed to generate an elliptical shaped first increment. This solution offers a gain in computational cost compared to the method of Feld-Payet [2010]; Feld-Payet et al. [2015]. Finally, a new insertion criterion was proposed to select both the most appropriate time to insert a given crack increment and the crack increment's surface. These new features enabled to successfully simulate the challenging case of a cup-cone failure in 3D.

The obtained results with the proposed continuous-discontinuous procedure have been validated in several ways. First, whenever possible, comparisons with simulations on a fixed very fine mesh were performed. Different minimum element sizes where tested in order to show the insensitivity of the procedure to the mesh size. Several crack increment lengths were considered in order to illustrate the insensitivity of the global response. This study also enabled to propose guidelines, based on the crack path representation, for the choice of the optimal crack increment length. Finally, in all cases, it was verified that the simulated crack surfaces were, in comparison with numerical simulations or experimental observations from the literature, in qualitative good agreement.

The proposed methodology thus enables to consider continuous-discontinuous transition for a vast number of cases. Let us note that among the Four proposed ingredients, three can be used with a different discontinuity representation: indeed, the simple geometrical approach to initiate elliptical crack shapes in 3D meshes, the new 3D insertion criterion and the reequilibrium procedure could all be used in combination with X-FEM or even to perform controlled element erosion. Thanks to all the proposed contributions, the cost associated to crack insertion can be significantly limited. However, the resolution of the mechanical problem still remains expensive in 3D: it would then be interesting to consider parallel computing. Let us note that, if several cracks were to appear in different areas, parallel computing could also be considered to perform the post-processing operations related to crack detection, insertion and propagation in each area independently. However, if a crack is to propagate at the interface of two adjacent domain, problem of consistency of the crack front may arise. Therefore, the cutting and the crack path tracking algorithm, must be adapted to the parallel framework in order to take full advantage of parallel computing.

Résumé en franc ¸ais

L'objectif principal est d'évaluer la capacité prédictive de la méthode développée. Pour ce faire, une base de données expérimentale issue d'une précédente thèse [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF] est utilisée pour effectuer des comparaisons entre simulations et tests expérimentaux. Dans un premier temps, le materiau étudié ainsi que la base de données expérimentales sont brièvement décrits. Ensuite, un modèle de comportement est développé afin de reproduire les phénomènes non linéaires observés dans les essais expérimentaux. Le modèle nonlocal à deux longueurs caractéristiques du chapitre 3 est utilisé dans la simulation présentée. Enfin, les simulations des essais expérimentaux utilisant l'approche continue-discontinue développée sont comparées aux résultats expérimentaux. Un bon accord est trouvé.

1 Inputs from to the previous work of [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF][START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] 1.1 Summary of the work

The aim the Ph.D. work of [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF] work was to develop and implement a robust and efficient numerical method enabling the numerical prediction, by the FEM, of ductile failure of car parts made of thin metal sheets during crash simulations. This numerical method was verified using a representative ductile material often used in the automotive industry: the DP450 "dual-phase" steel sheets. The first part of the work consisted in the characterization and the constitutive modeling of the DP450 steel. Characterization tests (approximately 120) has been carried out over a wide range of loading rates (quasi-static and dynamic), stress triaxialities, Lode parameter, and at different temperatures. Based on the exploitation of the results, a numerical local constitutive model have been established taking into account the different observed phenomena that have an influence on crack initiation and propagation: plasticity, strain-rate dependence (viscosity and self-heating) and damage [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF]. As the local coupled damage model suffers from mesh sensitivity, the second part of the work consisted in implementing two regularization methods: the first one based on the micromorphic [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF] approach and the second one based on the implicit gradient method [START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF]. Both approaches were implemented using an explicit resolution scheme adapted to dynamic simulations. After comparison, the implicit gradient method appeared to be more practical, robust and easy to identify and was retained as the best solution. Simulations were then carried out on several geometries, with dynamic loading, for which mesh independence (size or orientation) was demonstrated. In addition, simulations were generally in good agreement with experimental results [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF].

In this chapter, the experimental database developed in the work of [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF] is used, in order to evaluate the predictive capability of the continuous-discontinuous approach presented in the previous chapter. This chapter is structured as follows: First, in the section 1.2, the studied material is presented, as well as the experimental tests used in the characterization of the mechanical behavior. Then, the main experimental results are discussed in section 1.4, and the constitutive equations used to model the observed phenomena are presented in section 2. It should be noted that a model, adapted to the quasi-static part of the experimental database, and very similar material parameters are used compared to those used in the work of [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF]. Finally, simulations of different specimen geometries under quasi-static loading using the proposed continuousdiscontinuous method are performed. In these simulations, meshes with minimal element size obtained from the mesh study presented in [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF] are used. In addition, comparison of global responses and crack paths between experiments and simulations are presented in section 4.

Material: DP450 "dual-phase" Steel

The material of this study is a dual-phase DP450 steel which is a material often used in the automotive industry. This material is suitable for the metal forming process due to its large maximum elongation and its low yield stress. This material also provides good shock absorption which is essential in car crash conditions. This DP steel is composed of a small amount of martensite within a ferrite matrix (see fig. 5.1). Alloying elements are given in tab. 5.1 The material have been supplied as a 1.18mm thick sheet obtained by rolling process. In the following the rolling direction will be referred to as L (0 • ), the transverse direction as T (90 • ) and the diagonal direction as D (45 • ). with virtual extensometers. Measurement points for some geometries are marked with • on 5.2. All tests were repeated two or three times to check reproducibility, which was found to be very good. The summary table of the performed tests for all specimens, loading rates and temperatures is given in table 5.2.

Main experimental results

From experimental observations, the phenomena that have an influence on crack initiation and propagation:

• Plasticity : isotropic with hardening,

• Strain rate sensitivity : increase in strength with the increasing strain rate,

• Self-heating : heat generation at high strain rate that leads to a loss of macroscopic ductility,

• Damage : nucleation, growth and coalescence of cavities responsible for fracture.

The details of the full experimental analysis are presented in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF].

Constitutive equations

In this section, the constitutive equations used in the model are detailed. A simpler model than the one developed in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF][START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF] is chosen because only the quasi-static part of the experimental data is considered. [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF]) [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] 0.002mm.s -1 0.003mm.s -1 0.004mm.s -1 0.009mm.s

-1 NT1 2 L 2 T - - - NT2 2 L 2 T - - - V45 2 L 2 T • - - - M- shape - 2 L 2 T - - CCP - - 2 L 2 T • FN - - - 2L
2

.1 Hardening model

To model the hardening observed on the experimental results, Davaze proposed a modified Voce law [START_REF] Voce | The relationship between stress and strain for homogeneous deformations[END_REF] with a linear hardening term:

R(κ) = R e + Hκ + Q(1 -e -bκ ) (5.1)
where R e is the initial tensile yield stress, H the linear hardening coefficient, κ is the cumulated plastic strain, Q and b are the two Voce parameters. Their value is set according to the work of [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF]] see table 5.3

Nonlocal damage model

In the experimental study from [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF], it has been observed that failure of DP450 steel is due to the nucleation and growth of cavities. To represent these phenomena, the Gurson-Tvergaard-Needleman (GTN) coupled damage model [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] was used, with the effective stress of the GTN model written as:

   F (σ * , R) = σ * -R G(σ, σ * , f * ) = σ 2 eq σ 2 * + 2q 1 f * cosh( 3q 2 2 σm σ * ) -1 -(q 1 f * ) 2 def. σ * = 0 (5.2)
where σ eq is the von Mises equivalent stress and f * is a function of the porosity and defined as:

f * = f if f < f c f c + δ(f -f c ) if f ≥ f c (5.3)
In the work of [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF], a local GTN model was identified for a fixed mesh size to fit the experimental tests. Then in the subsequent work of [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF], an implicit gradient approach with regularization of the cumulative plastic strain κ was proposed. For this nonlocal model, a characteristic length (equal to 200µm) has been identified.Let us note that all material parameters of the local model had been preserved for the nonlocal version, except the controlled strain nucleation parameters: A n and κ c , in order to fit the experimental data base. This model has been implemented in an explicit resolution scheme in order to be suitable for dynamic simulations. The high-speed loading tests from the experimental database were simulated, with the implicit gradient model, and good agreement with the experiments was found.

In this section, the nonlocal model two characteristic lengths introduced in chapter 3 is adapted to simulate the low-velocity tests from the experimental database. In addition to nucleation and growth porosity, shear damage, as introduced in [START_REF] Nahshon | Modification of the Gurson model for shear failure[END_REF], is also used to properly describe the development of damage at low triaxiality, which is the case for the M-Shape specimen. In summary, the total porosity is expressed as follows:

f = f n + f g + f sh
(5.4) In this equation, the quantity e * ,incr scalar is defined as:

||e * ,incr scalar || = Ω (f * -f h ).(∆f * -∆f h )dΩ 1 2 =   Ω E ||e * ,incr scalar E || 2   1 2
(5.11)

The smoothed fields f * and ∆f * are calculated using the SPR technique [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF]. The accuracy of the error obtained from the SPR technique inside the domain is significantly better than that obtained near the external boundary [START_REF] Ródenas | Improvement of the superconvergent patch recovery technique by the use of constraint equations: the spr-c technique[END_REF] (see figure 5.5).

Figure 5.5: Cylinder subjected to internal pressure. Error estimator using the recovered stress σ * 22 (from [START_REF] Ródenas | Improvement of the superconvergent patch recovery technique by the use of constraint equations: the spr-c technique[END_REF]).

In the case of a thin metal sheet with 3D meshes, the majority of the nodes belong to the external boundary. This lead to a poor accuracy of the error estimator overall. To overcome this problem, [START_REF] Ródenas | Improvement of the superconvergent patch recovery technique by the use of constraint equations: the spr-c technique[END_REF] proposed a modified superconvergence patch recovery technique called SPR-C. This modified technique consists in applying an additional constraint in the classical 2.108 construction of the smoothed field. This additional constraint ensures that the smoothed field respects the boundary conditions of the mechanical problem. In [START_REF] Ródenas | Improvement of the superconvergent patch recovery technique by the use of constraint equations: the spr-c technique[END_REF], this constraint has been added using a Langrange multiplier.

Here, another approach, only adapted to 3D thin sheet meshes is proposed. This solution is based on the assumption that for a thin sheet simulation using a 3D mesh, the variation of the considered field f1 is small in the sheet thickness.

So at first, the scalar field f is projected into the 2D mesh used to extrude the thin sheet 5.6. The error is then evaluated on this 2D mesh, with less apex nodes connected to the

r E = ε 1/q 0 η 2 2q+d E Ω E η 2d d+2q E 1/2q
(5.12) Finally, the 2D map size is extruded uniformly along the thickness of 3D sheet mesh so to obtain the 3D map size needed for the mesh generator meshGems [START_REF] Frey | Mesh generation: Application to finite elements: Second edition[END_REF].

Continuous-discontinuous parameters

Using the definition of the bandwidth introduced in chapter 3, the obtained bandwidth for all specimen 5.2 is approximately: l b = 400µm. Most of the continuous-discontinuous parameter depend on the localization bandwidth and are summarized in 5.4: Note that with the chosen local criterion eq.5.12, if the estimated error is too low, then the elementary size factor r E will be unreasonably large. In order to avoid this artifact, a maximum element size h min is imposed during the remeshing process. Considering the identified parameters (see table 5.3), the element failure is obtained for f = f r = 0.2. However, the critical value f crit for crack insertion is chosen smaller than f r (i.e. 0.16), so to obtain a better convergence rate after crack insertion: because for bigger values of f crit , the crack insertion is delayed and the mesh becomes extremely stretched near the crack tip. 5.4.

In the presented simulations the crack increment length ∆a and f crit are chosen only considering numerical reasons: (i) to avoid frequent remeshing; (ii) to obtain an accurate crack direction and (iii) preserve an overall good convergence rate. However, if a correlation image analysis were performed on experimental test coupled with a crack path tracking algorithms [START_REF] Feld-Payet | Crack path tracking from full field measurements: a novel empirical methodology[END_REF], the crack tip and propagation velocity could be precisely estimated. In this case, the increment crack length ∆a and f crit could be experimentally calibrated.

Quasi-static simulations of the experimental data base

In this section, simulations of the quasi-static tests are carried out using the proposed continuous-discontinuous approach (see flowchart 2.27). Experiments and simulations are systematically compared in order to evaluate the prediction capability of the proposed continuous-discontinuous approach.

Inputs

In order to simulate all the sample geometries presented in figure 5.2, the constitutive model is implemented following a fully implicit scheme in the object-oriented FE software Z-set [START_REF] Foerch | Polymorphic constitutive equations in finite element codes[END_REF][START_REF] Besson | Large scale object-oriented finite element code design[END_REF]]. The material model and the implicit gradient finite element are implemented according to the formalism presented in 3. The simulations are performed using unstructured tetrahedral meshes generated with meshGems [START_REF] Frey | Mesh generation: Application to finite elements: Second edition[END_REF]. The meshes consist of 10-nodes elements with 4 Gauss points.

The usual symmetry conditions are used so that half of the samples are meshed (half thickness). The other symmetries are not taken into account in order to be able to represent the experimentally observed non-symmetrical crack paths for the M and FN shaped specimens.

Results

Quasi-static simulations at room temperature are presented and results are shown in figure 5.7 (standard flat), figure 5.8 (NT1), figure 5.9 (NT2), figure 5.10 (V45), figure 5.11 and figure 5.12. An overall good prediction is observed for all specimens. For the load-displacement curves, experimental scatter is represented by the filled area; this scatter corresponds to both material anisotropy and actual scatter (which is very limited). Experimental macroscopic transverse strains curves obtained by shadow tracking (see figure 5.4) are also compared in figures 5.7b, 5.8b, 5.9b, 5.10b. In cases where cracks are initiated at the notch root, transverse strains can only be evaluated before the onset of cracking as the shadow tracking technique is not accurate enough to detect the crack tip. This situation prevails in all cases except for the standard flat specimens in which cracks are initiated at the center of the specimen. For all simulations, the porosity field f after complete failure is plotted.

In all cases a very good agreement is obtained between simulated and experimental loaddisplacement/notch opening curves up to full specimen failure (i.e. load = 0). This first shows that the hardening model is appropriately identified. Crack initiation location is also well represented for all specimen: in particular initiation at the center of the specimen for standard flat specimens is well captured by the model (5.7c). However, in figure 5.7b, the engineering stress is plotted as a function of the transverse strain. The same results are obtained with the same parameters (taken from Davaze's work) and same loading on the standard flat specimen without remeshing. This to say, that the early fracture predicted from in the simulation is not due the continuous-discontinuous procedure but to the parameters fit.

Therefore, the damage parameters needed to be re-adjusted taken into account the whole experimental database so that to obtain a better fit.

M-shape simulated failure given in 5.11b, occurs almost simultaneously at both shear ligaments.

Due to its geometry, four crack initiation spots exist in the FN specimen. This leads to crack path dissymmetry and to several possible paths. The crack initiation locus then depends on the boundary conditions (due to e.g. a small misalignment in the griping system or specimen positioning). Note that, these differences have no influence on the macroscopic load-displacement curves. To be as representative as possible of the actual experimental conditions, displacements measured using DIC were used as boundary conditions on horizontal lines located at ∓17mm mm from the center line of the specimens. The simulation is then able to represent the observed crack path for the FN specimen (se figures 5.12b and 5.12c).

Let us note that for all the presented simulations, the meshes are only refined in the area where the damage is developing thanks to the use of the incremental error estimator on the total porosity variable f . In addition, mesh adaption allows the opening of the crack lips to be well represented in all simulations, and to obtain an accurate crack path thanks to the remeshing procedure with discrete crack insertion.

Note that the CCP simulation was note presented due to convergence problems encountered during the propagation process. Even with the use of the reequilibrium step which successfully preserved a good convergence rate for all other tests, the CCP simulation diverges during the reequilibrium step after the fifth crack increment insertion. The cause of the divergence is still not solved. This is all the more surprising that, in principle the CCP specimen is relatively close to the V45 specimen which was successfully simulated. .9: For the NT2 specimen, the load-notch opening is obtained from simulation and experiment are compared. Also, the map of total porosity f t is displayed at complete failure. .12: For the FN specimen, the load-displacement obtained from simulation and experiment are compared. Also, the map of total porosity f t is displayed at complete failure.

Conclusion

In order to evaluate the prediction capability of the developed continuous-discontinuous procedure, the experimental data base of [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] is used for experimentalsimulation comparisons. Only the quasi-static tests were considered in the study. At first, a simplified version of the constitutive equation, considering only non linear phenomena that occurs at low speed loading tests, was proposed. The implicit gradient nonlocal approach developed in 3 was used in order to obtain mesh insensitive results during the remeshing process. The identified parameters and characteristic length used in [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF] were, as much as possible, reused in the presented simulations. Note that, the identified parameters in [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF] considered only the dynamic part of the experimental data base, and in this chapter their pertinence is evaluated for the static part of the same experimental data base.

In addition, a remeshing strategy, coupled with an incremental error estimator, was adapted for thin sheet simulations with 3D meshes and then used in the proposed continuousdiscontinuous procedure. This remeshing strategy is more adapted than the active process zone remeshing procedure which was not adapted for the FN simulation for example, because of the size of the active process zone.

Finally, the result of the simulations were compared to experimental responses and good agreement was found for all specimen except for the standard flat specimen; where the crack initiation was not well predicted. This is due to the fact that, the chosen parameters from Davaze work did not take into account the quasi-static experiments. Therefore, a better fit is necessary for the damage parameters in order to obtain good correlation for all the experimental data base.

The opening of the crack lips was well represented due to the crack insertion into the simulation mesh. Also optimal meshes were produced thanks to the use of the incremental error estimator based on the total porosity variable. However, convergence problems were encountered for the CCP specimen even though this test is very similar to the V45 which was successfully simulated with a good convergence rate Furthermore, because the crack lips opening and the crack tip position were well presented in the simulations using the proposed continuous-discontinuous procedure, experimental calibration of the parameters (f crit , ∆a) can become possible if the crack tip position could be accurately measured, e.g. following the procedure described Feld-Payet et al. [2019]. In any case, these results already constitute an encouraging step toward a better understanding and prediction of the ductile failure for laboratory tests then for industrial structures.

Chapter 6

Conclusion and perspectives

This thesis was dedicated to developing a general framework for ductile crack initiation and propagation over long distances. In order to achieve this goal, the developed framework is based on the continuous-discontinuous modeling of ductile failure. First, the ductile failure process is described using a physically-based model, the GTN model, in order to model the damage of the structure as a continuous scalar field. The continuous aspect of the procedure was the main subject of chapter 3, with the development of a nonlocal GTN model with two characteristic lengths. Then, once damage reaches a critical value, a transition from the failure modeling as damage to a discrete crack is carried out using mesh adaption. This transition can be used to describe both crack initiation and propagation. Indeed, as the damage continues to grow ahead of the initial crack, it drives the crack propagation process. The transition from the continuous model to a discreet crack using mesh adaption was presented in chapter 4. Finally, the developed approach was applied in chapter 5 to an existing experimental database in order to evaluate its predictive capability.

1 Implicit nonlocal GTN model with two characteristic lengths adapted for mesh adaption

In this work, a nonlocal GTN model based on an implicit gradient formulation (applied to both the accumulated plastic strain that drives void nucleation and the plastic volume variation driving void growth). The implementation of this model was detailed for the updated Lagrangian framework. Indeed, throughout the failure process, large plastic strains are present in the localization bands where cracks would initiate. Thus, the use of finite strain formulation becomes a necessity. The implementation of the model was validated with a thorough mesh convergence study proving mesh insensitivity was achieved. The model parameters were fitted to reproduce the global response for tests on a pipeline steel using a tensile bar and a plane strain specimen. Experimental results were taken from the literature. Each test is representative of a characteristic crack path: cup-cone fracture (tensile test) and slant fracture (plane strain test). Some guidelines regarding the choice of the characteristic length were proposed. For that purpose, a novel methodology was proposed to measure the localization bandwidth free from discretization error. Then by applying this methodology, a linear relationship between the characteristic length and the localization bandwidth was identified. These results were then used to select the proper discretization for given values of the characteristic lengths to obtain converged solutions. Finally, the effect of the characteristic length was demonstrated on the formation and shape of the cup-cone crack path. This study concluded that the cup-cone crack path should be affected by the size of the specimens. Pure flat fracture should be favored when using (very) small specimens. It could theoretically be possible to identify the internal length directly. This result paves the way for a more direct determination of material lengths by using homothetic specimens so as to obtain different crack paths. The fitted lengths should then be able to represent both the macroscopic behavior and the crack path for all sample sizes.

2 Damage crack transition using mesh adaption

The proposed continuous-discontinuous procedure is based on the method presented in Feld-Payet [2010] using mesh adaption in order to introduce a discrete crack into zones that have reached a critical value for damage. Some new contributions were proposed and implemented in the Z-set software to enable a more robust and cost-efficient transition between a continuous model evaluating local material degradation and a model representing the propagation of a discrete crack. (i) A new procedure was proposed in order to facilitate the return to equilibrium and to improve the convergence rate after a change of discretization (remeshing). It has been shown that this reequilibrium procedure limits the number of increments during crack propagation, which contributes to the cost-efficiency of the method.

(ii) Another new element is the refinement of only the active process zones. This feature has a double advantage: not only are the degradation phenomena leading to complete failure are well described, but also the number of elements (and thus the computational cost) remains approximately constant. (iii) Besides, two contributions related to 3D simulation were proposed. The first one deals with crack initiation entirely inside the structure: (a simple geometrical crack initiation methodology based on SVD analysis of the coordinates of broken points has been proposed to generate an elliptically shaped first crack increment.

(iv) Finally, a new insertion criterion was proposed to select both the most appropriate time to insert a given crack increment and the crack increment's surface. The proposed methodology thus enabled to consider continuous-discontinuous transition for a vast number of cases and enabled to successfully simulate the challenging case of a cup-cone failure in 3D.

3 Toward a predictive continuous-discontinuous model for ductile failure

To evaluate the predictive capability of the proposed continuous-discontinuous approach to failure, the experimental database constructed in the work of [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF] is used to conduct experiments-simulations comparisons.

Only the quasi-static experimental tests were considered. For the continuous part of the approach, the nonlocal model detailed in chapter 3 was adapted and reused with the same parameters calibrated in the work of [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF]. Note that, the calibrated parameters in [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF] only considered the dynamic part of the experimental database and its relevance for the quasi-static experiments is evaluated in this work.

Concerning the continuous/discontinuous transition, the strategy detailed in 4 is used to perform the presented simulations. A more suitable remeshing strategy driven by an incremental error estimator was considered. Indeed, because of the significant size of the APZ, the remeshing strategy using the incremental error estimator produced more suitable meshes, thus a better cost-precision trade-off.

Finally, the experimental results were simulated using the proposed approach, and a good agreement, with respect to the global responses and crack paths, was found for all specimens except for the prediction of the crack initiation for the standard flat specimen. This result indicates that the parameter fit, carried out in the work of [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF], on the "dynamic" experiments of the database can be improved to reproduce with increased accuracy crack initiation and propagation for the entire experimental database.

Future work

The presented work and obtained results demonstrate the accuracy and robustness of the proposed continuous-discontinuous approach. Nevertheless, several improvements could be proposed to enhance the efficiency of the proposed method. Thanks to the results in chapter 3, a proposition was made to possibly identify the characteristic length(s) using homothetic specimens and studying the change from a straight crack path to a cup-cone crack path when increasing the specimen size. Indeed, the choice of the characteristic length may be restricted due to numerical reasons as using very small lengths would require significant computational efforts.

In chapter 4, even though the fine remeshing is only applied to the APZ, the 3D simulation of the axisymmetric specimen was conducted using a larger minimal mesh size than the one used for 2D simulation due to the size of the computation. To overcome this issue, several improvements could be the subject of future work:

1. Adaptation of the continuous-discontinuous procedure to parallel computing. Parallel computing and parallel remeshing are problems that are well addressed in the literature; however, the building brick needed for the continuous-discontinuous procedure needs to be adapted for the parallel computing framework. In particular, the question of robustness of the crack path tracking algorithm should be addressed when the crack crosses several decomposed domains.

2. The sequential mumps solver was used for the resolution of the mechanical problem. With these types of solvers, significant computational resources are needed.

Another less expensive alternative is the use of an iterative solver. However, a good preconditioner adapted for remeshing is then essential to ensure an efficient convergence rate.

3. Using a lower approximation for the finite element by adapting the "bubble elements" to the nonlocal framework. This could be well adapted for less critical industrial applications. Indeed, P1+/P1/P1 nonlocal element would provide a less accurate result, but larger structural simulations could be performed.

Another challenge was encountered in the simulations of the experimental database. One particular simulation of the CCP specimen was not possible due to divergence during the reequilibrium step of the resolution algorithm after the fifth crack insertion. A possible explanation is a significant transfer error due to the volumetric locking problem. Therefore, a locking-free nonlocal element formulation would be needed to eliminate this source of error. This also means additional nodal field(s) of unknowns in the finite element formulation. So in order to preserve the efficiency of the proposed approach, a bubble formulation is essential, especially for industrial applications.

Thanks to the proposed continuous-discontinuous procedure, the position of the crack front is explicitly known throughout the simulation. This information could be used in the fitting procedures of the models if the crack position and shape are experimentally determined.

• Derivatives of R n :

∂R n ∂∆f n = 1 (A.17)

• Derivatives of R ω :

∂R ω ∂∆ε e = -(1 -f )∆κI : N : E (A.18)

∂R ω ∂∆κ = -(1 -f )trace(n) (A.19) ∂R ω ∂∆f g = -(1 -f )∆κ trace(n f ) (A.20) ∂R ω ∂∆f n = -(1 -f )∆κ trace(n f ) (A.21) ∂R ω ∂∆ω = 1 (A.22)
The partial derivatives of σ are obtained considering the stationarity of the G function (eq. 3.3) : G = 0 and δG = 0. One then has: .23) and .24) and

χ = 1 h = ∂G ∂σ , n = ∂σ ∂σ = -h ∂G ∂σ = -hν, ∂σ ∂f = -h ∂G ∂f , ( A
N = ∂n ∂σ = -h ∂ 2 G ∂σ 2 -h 3 ∂ 2 G ∂σ 2 ν ⊗ ν + h 2 ∂ 2 G ∂σ∂σ ⊗ ν + ν ⊗ ∂ 2 G ∂σ∂σ (A
n f = ∂n ∂f = -h ∂ 2 G ∂σ∂f + h 2 ∂G ∂f ∂ 2 G ∂σ∂σ + h 2 ∂ 2 G ∂σ ∂f -h ∂ 2 G ∂σ 2 ∂G ∂f ν (A.25)
Appendix B

Consistent tangent matrix

Once the solution of the system of equations A.1 to A.5 is found, it becomes possible to numerically compute the consistent tangent matrix. For any given small variation of the input variables δV IN , the state variables will be modified so that R remains null. Indeed, in the case of κ and ω, the output and state variables are similar. All other terms are null.

The calculation of K uu follows usual derivations for finite strain formulation and leads to: The term K κu is the most complex of the three terms to calculate. Indeed it requires to calculate the derivative of Ḟ κ with respect to the displacements, which involves three terms: -J κN , Jl 2 l ( ĠT .G + G T . Ġ).κ e and J (κ -κ)N + l 2 κ G T .G.κ e def = JT k . Considering the first term, one has:

K uu = Ω e
-J κu N = -J ∂∆κ ∂∆ε .B D . ue
To compute the second term, one first notices that G = F -T .G 0 (note that in the case the Voigt notation is not used).

3. 9

 9 Band width measure for the nonlocal volume variation ω when ω max = 0.5. Evolution of the local and nonlocal volume variations along the symmetry axis. The double arrows show the measure of the corresponding band widths l ω b and l ω b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 (a) Definition of H r1 and H r2 . (b) Contours of the nonlocal volume variation ω when ω max = 0.5. Arrows here define current element heights H r1 and H r2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.3 Use of the marching ridges algorithm for 3D structures with 2D polar grid search in planes orthogonal to the current crack front segments, figure a), as proposed in Feld-Payet et al. [2015]). . . . . . . . . . . . . . . . . . . . . . 93 4.4 The proposed insertion criterion for crack propagation in 3D structures compares the area of the crack increment with constant length for each segment, A unif , and the area obtained with different lengths, A real : insertion is allowed only when A real ≥ A unif . . . . . . . . . . . . . . . . . . . . . . . 94 4.5 Automatic detection of the active process zone (APZ) for the remeshing procedure: the center and radius definitions change a) before and b) after crack initiation. Using this procedure, only the active process zone ahead of the crack front is finely meshed. The mesh used for areas already crossed by a crack are coarser depending on their distance to the crack front. . . . . . . . 95

FDΣE 2 N T 4 4 P

 244 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Transpose operator of a second order tensor as: a = a T ij = a ji A = a ⊗ b By definition : A ijkl = a ij b kl A = a⊗b By definition : A ijkl = a ik b jl A = a⊗b By definition : A ijkl = a il b jk ȧ Rate operator of a variable with respect to time : ∂a ∂t trace(a) Trace operator of a seconde tensor variable: a ii δa Differential operator ∂a ∂b Partial derivative of a variable a with respect to b ∇a First spacial gradient of a scalar variable: ∇a = ∂a ∂x i ∇ n a N th spacial gradient of a scalar variable: ∇ n a = ∂ Lode operator of a second order tensor θ L (.) Lode angle of a second order tensor Bandwidth of a localized variable l κ Characterestic length with respect to the variable κ l ω Characterestic length with respect to the variable ω G c Elastic part obtained fromm an decomposition of the strain tensor ε p Plasric part obtained fromm an decomposition of the strain tensor Finite strain symbols x Position of a material point on the current configuration X Position of a material point on the initial configuration Transformation Symetric part of the velocity gradient tensor: D = 1 2 (L + L T ) W Unsymetric part of the velocity gradient tensor: W = 1 2 (L -L T ) Q rotation tensor of the corotaional frame: Q = W .Q with Q(t = 0) = 1 Cauchy stress in the unrotated frame: Σ = Q.σ.Q T Resolution scheme notations a n Varibale at the time t n : a n = a(t n ) ∆a Increment of a varible such as: a n+1 = a n + ∆a R a Residu related to the evolution equation of the variable a V IN T Variables to be integrated J Jacobian necessary for a Newton-Raphson scheme: R(V IN T , t n+1 ) = R(V IN T , t n ) + J.∆V IN T V IN Input variable to the material integration step V OU T Output variable to the material integration step K mat Consistent tangent matrix Error estimation e Discretisation error measured on the displacement field e σ Discretisation error measured on the stress field η ex Norm of the exact error η Norm of an the global error estimator η E k Norm of an the error estimator restricted to the finite element E k v The exacte value of a variable v h The approximation of the variable using finite element method γ Efficienty index of an error estimation e rdc Error estimator using the constitutive equation e * Error estimator using gradient recovery method e * ,incr Incremental error estimator using gradient recovery method η * ,incr Relative incremental error estimator using gradient recovery method f * Approximation of the variable f using super convergent patch recovery method Mesh adaption h E Size of the element E on the current mesh h * Size the element E on the new mesh r E Elementaty size factire: r E = Shape function i on an element at the current mesh φ new i Shape function i on an element at the new mesh X ip old Coordinates of an element integration points at the current mesh X ip new Coordinates of an element integration points at the new mesh ∆t remeshing User paramater to set remeshing frequency Damage to crack transition symboles f crit Critical value of a variable used for the crack insertion criterion M coord Cordinates of points that have reached the critical value R Radius of the search algorithm (Marching Ridges) ∆θ Precision of the search algorithm (Marching Ridges) m c Multipication factor used in the definition of the increment crack length: ∆a = m c * l b A unif Area of the uniform crack increment A real Area of the real crack increment m AP Z Factore to define the process zone from a localized varibale ∆t reeq Time increment for the reequilibrium step Thickness reduction of a plane strain specimen N h Number of the finite element in the diameter/thickness of a meshed specimen r h Ratio of the hight/width of a finite element at the refined zone of the mesh H r1 Height of the finite element at the onset of fracture of an element at the refined zone h min Size of finite elements at the refined zone of the mesh h max Size of finite elements at the coase zone of the mesh Notch Tensile round bar with a radius 10 N T 2 Totch Tensile round bar with a radius Notch Tensile round bar with a radius D. work is realized in the framework of a collaboration between the Centre des Matériaux (CdM) of the École Nationale Supérieure des Mines de Paris, Onera -the french aerospace lab, and Transvalor's Z-set development team.

4 . 3

 43 Mesh adaption procedure . . . . . . . . . . . . . . . . . . . . . . . 4.4 Transfer operator . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Reequilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 . 1 :

 21 Figure 2.1: Cross section of a damaged material

Figure 2 . 2 :

 22 Figure 2.2: Gurson elastic domain for different porosities.

Figure 2 . 3 :

 23 Figure 2.3: A representative general ellipsoidal void with semi-axes a 1 , a 2 , a 3 and the corresponding orientation vectors n (1) , n(2) , n(3) ; (b) a prolate void; and (c) anoblate void ( from[START_REF] Cao | A model for ductile damage prediction at low stress triaxialities incorporating void shape change and void rotation[END_REF]).

Figure 2 . 4 :

 24 Figure 2.4: Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for mesh orientations in the case of the axisymmetric specimen.The global responses are similar for all the orientations. However, the total porosity for each mesh size at total failure shows different a crack path for each mesh orientation, which indicates mesh dependence.

Figure 2 . 6 :

 26 Figure 2.6: Typical nonlocal weighting functions: Gaussian (a) and bell-shaped polynomial functions (b) (from[START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF]).

Figure 2 . 9 :

 29 Figure 2.9: Crack surface and damaged elements: (a) crack surface; (b) damaged elements from [Yun et al., 2019]. Controlled kill element was used to represent the crack in a fracture simulation of a quasi-brittle material.

Figure 2 .

 2 Figure 2.10: Scheme displaying a) intra-element cracking and b) inter-element cracking.

Figure 2 .

 2 Figure 2.12: Damage distribution reported on the deformed configuration from [Feld-Payet, 2010]. Intra-element cracking with adaptive remeshing for crack insertion.

Figure 2 .

 2 Figure 2.16: polar grid search around the crack front in 2D from [Feld-Payet et al., 2015]

Figure 2 .

 2 Figure 2.17: Schematic representation from [Feld-Payet et al., 2015], of the damage profile plotted for two different evaluation radii (left) and the corresponding circles centered at the point of maximum damage, thus already on the ridge (right)

Figure 2 .

 2 Figure 2.21: Scheme of the median-axis-search from [Tamayo-Mas and Rodríguez-Ferran, 2015]

Figure 2 .

 2 Figure 2.24: (a) Visualization of the optimal sampling points of a function (b) and its gradient for 1D (quadratic ) (from[START_REF] Zienkiewicz | The finite element method[END_REF]).

Figure 2 .

 2 Figure 2.25: Order of super-converging sampling points for different C 0 elements, (from [Zienkiewicz and Taylor, 2000])

Figure 2 .

 2 Figure 2.26: super-convergent point on patches; sampling point for the smoothed field (from [Zienkiewicz and Taylor, 2000])

  Figure 3.1: Dimensions and boundary conditions for the tensile test simulations of the axisymmetric specimen 3.1a and the plane strain specimen 3.1b.

Figure 3 . 2 :

 32 Figure 3.2: Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for the axisymmetric specimen (left) and the thickness reduction -∆e/e 0 for the plane strain specimen (right). The comparison with global curves from the literature (respectively from[START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] and[START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]) indicates a correct parameter fit. Total porosity is displayed at total failure for both simulations: a cup-cone failure (left) and a slant failure (right) can be respectively observed.

Figure 3 . 3 :

 33 Figure 3.3: Evolution, for the nonlocal model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for the axisymmetric specimen with a central spatial discretization of N h = 90. The confrontation with global curves from the litterature[START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] and the local model indicates a correct parameter fit. Total porosity is displayed at various times of the simulation to illustrate the formation of the expected cup-cone crack path.

Figure 3 . 4 :

 34 Figure 3.4: Evolution, for the nonlocal model, of the engineering stress as a function of the thickness reduction -∆e/e 0 for the plane strain specimen with a central spatial discretization of N h = 90. The confrontation with global curves from the literature[Besson et al., 2003; Huespe et al., 2012; Leclerc et al., 2020] and the local model indicates a correct parameter fit. Total porosity is displayed at various steps of the simulation to illustrate the propagation of a flat crack instead of a slanted one. This illustrates the need to properly calibrate the nonlocal characteristic lengths.

Figure 3 . 5 :

 35 Figure 3.5: Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for different spatial discretization in the case of the axisymmetric specimen. The similarity of the global responses does not enable to observe a dependence on the mesh size. The total porosity for each mesh size at total failure shows a cup-cone failure with a slight deviation near the axis of symmetry.

Figure 3 . 6 :

 36 Figure 3.6: Evolution, for the nonlocal model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for different spatial discretization in the case of the axisymmetric specimen. The global responses are similar for the mesh sizes N h > 20. The total porosity for each mesh size at total failure shows a cup-cone failure with a straight first horizontal part except for the coarsest mesh size N h = 20, which exhibits a flat crack path. The finest mesh size N h = 180 enables the crack path to bifurcate twice, which can also be experimentally observed. These results indicate mesh independence.

Figure 3 . 7 :

 37 Figure 3.7: Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for various mesh orientations in the case of the axisymmetric specimen. The global responses are similar for all the orientations. However, the total porosity for each mesh size at total failure shows a different crack path for each mesh orientation, which indicates mesh dependence.

  Figure 3.8: Evolution, for the nonlocal model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for various mesh orientations in the case of the axisymmetric specimen. The global responses are similar for all orientations. Besides, the total porosity for each mesh size at total failure shows an identical cup-cone crack path for each mesh orientation. These observations indicate mesh independence.

Figure 3 Figure 3 .Figure 3 .

 333 Figure 3.9: Band width measure for the nonlocal volume variation ω when ω max = 0.5. Evolution of the local and nonlocal volume variations along the symmetry axis. The double arrows show the measure of the corresponding band widths l ω b and l ω b .

  Figure 3.12: Evolution of the estimated band width free from discretization error for both the axisymmetric and plane strain tensile tests as a function of the characteristic length (l κ = l ω ). There appears to be a proportional relation between the band width and the characteristic length, with slightly different coefficients for each geometry. This means that it is possible to estimate the size the band width for a given characteristic length and a given geometry.

  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

Figure 4 . 3 :

 43 Figure 4.3: Use of the marching ridges algorithm for 3D structures with 2D polar grid search in planes orthogonal to the current crack front segments, figure a), as proposed in Feld-Payet et al. [2015]).

Figure 4 . 4 :

 44 Figure 4.4: The proposed insertion criterion for crack propagation in 3D structures compares the area of the crack increment with constant length for each segment, A unif , and the area obtained with different lengths, A real : insertion is allowed only when A real ≥ A unif .

Figure 4 . 5 :

 45 Figure 4.5: Automatic detection of the active process zone (APZ) for the remeshing procedure: the center and radius definitions change a) before and b) after crack initiation.Using this procedure, only the active process zone ahead of the crack front is finely meshed. The mesh used for areas already crossed by a crack are coarser depending on their distance to the crack front.

  and the linear relation found between the characteristic length l c and the bandwidth size l axi b model, there are two characteristic lengths. Only the smallest one (i.e. l κ = 20µm, see the previous chapter, section2) should be considered to evaluate the bandwidth. The smallest bandwidth should thus be, according to equation 4.3: l b = 45.6µm.

  0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0h min = 40 m h min = 20 m h min = 10 m meshh min = 10 m meshh min = 20 m meshh min = 40 m meshfixedmesh, N h = 180

Figure 4 . 6 :

 46 Figure 4.6: Evolution of the engineering stress as a function of the diameter reduction -∆φ 0 /φ for different minimal mesh sizes obtained when simulating a tensile test on an axisymmetric specimen. With minimal mesh sizes h min ∈ {10, 20, 40}µm, the proposed continuous-discontinuous method has been used whereas the simulation with a very fine mesh is without any mesh adaption or crack insertion. For all the considered meshes, the global responses are similar and the total porosity at complete failure shows a cup-cone failure. This similarity indicates that the proposed procedure is rather robust to the choice of the mesh size.

Figure 4 . 7 :

 47 Figure 4.7: Evolution of the engineering stress as a function of the diameter reduction -∆φ 0 /φ for different crack increment lengths ∆a, using the proposed continuousdiscontinuous method, in the case of the axisymmetric specimen. It appears that ∆a has no effect on the global response. However, the largest increment length, ∆a = 180µm, leads to a cup-cone crack path which is not completely smooth compared to the other crack increment lengths.

Figure 4

 4 Figure 4.9: Evolution of the cumulated number of time increments as a function of imposed vertical displacement for the simulation of the axisymmetric tensile test on a fixed mesh (in black), with mesh adaption but without reequilibrium (in red) and with mesh adaption and the proposed reequilibrium procedure (in dark red). The higher number of time increments is obtained for the simulation with mesh adaption without reequilibrium. With the reequilibrium procedure, the total number of increment is approximately divided by two. The lowest number of time increments is obtained for the fixed mesh for which there is no discrete crack represented, so no transfer needed. In all cases, the crack propagation phase is accompanied by an increase in the number of increments.

Figure 4 .Figure 4

 44 Figure 4.10: Specimens' geometry specifications and boundary conditions.

  fixed mesh N h = 180 proposed procedure h min = 20 m

Figure 4

 4 Figure 4.13: Evolution of the engineering stress as a function of the thickness reduction -∆e 0 /e 0 of the plane strain specimen for both a simulation with the proposed continuousdiscontinuous procedure and a simulation with a fixed mesh. The global curves are similar up to the latest stages of the simulation where the simulation with the proposed continuousdiscontinuous procedure ends up with a larger final thickness reduction. Total porosity is displayed at various times of the simulation using the proposed method, and at total failure for the fixed mesh simulation. The proposed method predicts an axisymmetric slant crack path with a bifurcation whereas the fixed mesh simulation predicts a symmetric slant crack path. The bifurcation coincide with the deviation of the global responses.

Figure 4 .

 4 Figure 4.14: Evolution of the engineering stress as a function of Crack MOuth Displacement (CMOD) for two meshes of the CT specimen with different minimum element sizes: h min ∈ {20, 40}µm. The global response is identical for both meshes which tends to indicate mesh convergence. Contours of the total porosity f t are displayed at different times of the simulation.Before crack initiation (marked with a green cross), crack tip blunting is well described thanks to mesh adaption which enables to preserve good quality elements near the crack tip. Then the crack propagates over a relatively long distance with still well shaped fine elements near the crack tip and a coarsened mesh along the crack lips.

  Figure 4.15: Evolution of the engineering stress as a function of diameter reduction -∆φ 0 /φ for the simulation of a bar with round section with the proposed approach in 3D (in red), the proposed approach with 2D axisymmetric elements and a smaller minimum element size (in green) and a fixed mesh with 2D axisymmetric elements (in blue). All the global responses are similar, which thend to indicate the validity of the used 3D mesh. Contours of the total porosity are displayed at different times for the 3D simulation to illustrate the cup-cone crack propagation.

Figure 4 .

 4 Figure 4.16: a: initial mesh of quarter of the specimen for the 3D simulation. b: fracture surfaces obtained from the 3D simulation. The predicted cup-cone crack surfaces are smooth, well defined and qualitatively in agreement with what can be expected from the experimental results obtained by Besson et al. Besson et al. [2003].

  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5 . 1 :Figure 5 . 2 :Figure 5 . 3 :Figure 5 . 4 :

 51525354 Figure 5.1: Microstructure of DP450 steel after nital etching from [Davaze et al., 2020].

Figure 5 . 6 :

 56 Figure 5.6: Scalar incremental strategy adapted for thin 3D meshes.

  Load-Engineering strain curve at ( ε= 0.001 s -1 ) Load-Engineering trans strain curve at ( ε= 0.001 s -1 ) (c) f t map at total failure

Figure 5 . 7 :

 57 Figure 5.7: For the standard flat specimen, the load-Engineering transverse strain obtained from simulation and experiments are compared. Also, the map of total porosity f t is displayed at complete failure.

  Figure 5.8: For the NT1 specimen, the load-notch opening is obtained from simulation and experiment are compared. Also, the map of total porosity f t is displayed at complete failure.

  Figure 5.10: For the V45 specimen, the load-notch opening is obtained from simulation and experiment are compared. Also, the map of total porosity f t is displayed at complete failure.

2 κ

 2 {C τ } .B D dΩ + Ω e B T L .M R (Σ).B L dΩ = K M uu + K G uu (D.10)where K M uu is the part related to the material non linearity and K G uu the geometrical non linearity.Derivation of K uu :In order to compute K uu , one first considers variations with respect to u e . One can show thatB L = M R (F -1 ).B F and B T L = B T F .M R (F -T ) Proof: {L} = B L . ue = Ḟ .F -1 = M R (F -1 ).{ Ḟ } = M R (F -1 ).B F . ueOne then has:B T L = B T F .M T R (F -1 ) = B T F .M R (F -T ) The product ḂT L .{Σ} can then be written as:ḂT L .{Σ} = B T F .M R ( Ḟ -T ).{Σ} = B T F . Σ. Ḟ -T (D.11)One has:Ḟ -T = -(F -T ⊗F -1 ) : Ḟ T = -F -T ik F -T = -Σ.L T .F -T = -M R (F -T ). Σ.L T Finally ḂT L .{Σ} = -B T F .M R (F -T ). Σ.L T = -B T L . Σ.L Tand Ḟ uu can be expressed as (note that J = Jtrace(L)): Σ.L T + Σu + trace(L)Σ dΩ 0 (D.13)The Truesdell rate of Σ is:Σ τ = Σu -L.Σ -Σ.L T + trace({C τ } .{L} + M R (Σ).{L}) dΩ 0 {C τ } + M R (Σ)) .B L . ue dΩ 0 (D.18) = Ω e B T L . ({C τ } + M R (Σ)) .B L dΩ . ue (D{C τ } .B D dΩ + Ω e B T L .M R (Σ).B L dΩ = K M uu + K G uu (D.20)where B L was replaced by B D considering the symmetries of C τ .Then, to compute the K uκ term, let us focus on the part of the derivative of Σ that depends on κa sub block matrix of the consistent tangent matrix K mat (eq. B.4) as :∂Σ ∂κ = ∂σ ∂κ (D.24)The bloc K uω is computed in a similar way by replacing κ by ω from equation D.21 to equation D.24.Calculation of Kκu , K κκ , K κω Calculation of the second bloc line can be obtained by writing the variation of F κ : ĠT .G + G T . Ġ .κ e + Jl 2 κ G T .G. κe + J (κ -κ)N + l 2 κ G T .G.κ e dΩ 0 (D.25) and considering the derivative κ as the sum of the partial derivative of κ regarding the variables u e , κ e and ω e : κ = κu + κκ + κω (D.26)

  j κ e k = F -T ij G 0 jk κ e k so that G = F -T .G 0 . Therefore Ġ = Ḟ -T .G 0 remembering that Ḟ -T = -L T .F -T Ġ = -L T .F -T .G 0 = -L T .G and ĠT = -G T .LLet V be the operator linking the Voigt representation of a tensor to the matrix representation such that a = V.{a} ora ij = V ijk a k and V * such that a T = V * .{a}. Indeed V * ijk = V jik . Using this notation, G T . Ġ.κ e = -G T .L T .G.κ e = -G T . (V * .(B L . ue )) .G.κ e = -G T . (V * .(B L . ue )) . ∇κ G T ij Ġjn κ e n = -G T ij V * jkl B L lm ue m ∇ k κ = -G T ij V kjl B L lm ue m ∇ k κ = -G T ij ∇ k κV kjl B L lm ue m G T . Ġ.κ e = -G T .( ∇κ.V.B L ).u e Similarly ĠT .G.κ e = -G T .( ∇κ.V * .B L ). ue m and ĠT .G.κ e + G T . Ġ.κ e = -G T .( ∇κ.(V + V * ).B L ). ue = -2G T .( ∇κ.V.B D ). ueFinally, one now considers the last term JT k .JT k = Jtrace(L)T k = J({1}.B D . ue )T κ = JT κ ⊗ ({1}.B D ). ue (D.27)In the end, one finally gets: .( ∇κ.V)+ T κ ⊗ {1} .B D dΩ (D.28)where ∂∆κ ∂∆ε is computed using a sub block matrix of the consistent tangent matrix K mat (eq. B.4) as:Calculation of K ωu , K ωκ , K ωω Calculation of the last bloc line can be obtained exactly in the same way as the second bloc line by writing the variation of F ω . It leads to the same terms only replacing κ by ω.

  2.1 Cross section of a damaged material . . . . . . . . . . . . . . . . . . . . . . 2.2 Gurson elastic domain for different porosities. . . . . . . . . . . . . . . . . . 2.3 A representative general ellipsoidal void with semi-axes a 1 , a 2 , a 3 and the corresponding orientation vectors n

  Scheme of the median-axis-search from [Tamayo-Mas and Rodríguez-Ferran, 2015] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.22 3D crack initiation algorithm from

	List of Figures
	2.21

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

  ) . . . . . . . . . . . . . . . . . . . . . . . . 57 2.26 super-convergent point on patches; sampling point for the smoothed field (from [Zienkiewicz and Taylor, 2000]) . . . . . . . . . . . . . . . . . . . . . 58 2.27 Flowchart summarizing the proposed continuous-discontinuous procedure with a nonlocal GTN model and discrete crack lips representation using mesh adaption . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

  Simulation of ductile failure est expliqué puis illustré par une simulation sur le modèle local GTN. A la fin de cette première partie, le modèle GTN est choisi comme modèle de dommage à utiliser pour ce travail. Ensuite, les approches de régularisation des modèles de dommages sont exposées avec leurs forces et faiblesses. Enfin, la méthode du gradient implicite est choisie comme méthode de régularisation à utiliser pour le modèle GTN. L'approche de transition entre l'endommagement et la fissuration est présentée et chaque composant de cette méthode est détaillé. A la fin de cette section, une approche continue-discontinue basée sur des techniques de remaillage est choisie pour simuler la rupture ductile. À la fin du premier chapitre, un cadre pour la simulation de la rupture ductile est développé, les composants existants du cadre sont présentés, et les composants qui doivent être développés sont renvoyés au chapitre qui les détaille.

	1 Ductile failure
	1.1 Ductility and ductile fracture

Une revue de la littérature sur les modèles d'endommagement continus et les modèles micromécaniques pour la rupture ductile est présentée. Le problème de localisation qui est l'inconvénient majeur de ces modèles d'endommagement couplés dans leur forme locale Chapter 2. Bibliography: Materials can be classified according to their ductile or brittle behavior: generally ceramics, glasses, concretes and composites are grouped in the category of brittle materials category while metals and polymers belong to the category of ductile materials. This classification is not absolute: at high temperatures a material can tolerate higher deformations and is thus more ductile. At very low temperatures metals (such as ferritic steels) or polymers can become brittle. Other factors can affect the ductility of a material such as corrosion and irradiation.

  from: A. El Ouazani Tuhami, S. Feld-Payet, S. Quilici, N. Osipovc, J. Besson. A two characteristic length non-local GTN model: application to cup-cone and slant fracture, in preparation.

	Chapter 3
	Implicit gradient nonlocal model with
	two characteristic length

Contents 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Material models and their framework . . . . . . . . . . . . . . . . . . . 2.1 Material models . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Finite strain framework . . . . . . . . . . . . . . . . . . . . . . . . Taken

Table 3 .

 3 2: Material parameters for the local and nonlocal models.

		local model nonlocal model
	κ c	1.2	1.2
	A n	0.2	0.4
	l κ = l ω	-	40µm

Table 4 .

 4 1: The parameters of the the continuous-discontinuous procedure. Notation l b designates the localization bandwidth that can be measured following the methodology described inEl Ouazani Tuhami et al. [accepted].

	Remeshing parameters
	Minimal element size	h min	l b /2
	Size of APZ	m AP Z	2
	Remeshing frequency	∆t remeshing	adjustable
	Marching ridges algorithm
	Evaluation radius	R	4l b
	angular precision	∆θ	5 •
	Insertion criterion
	Critical value for insertion	f crit	0.95f c = 0.6333
	Crack increment length	∆a	4l b
	Reequilibrium time increment	∆t reeq	10 -4 s
	3.1 Robustness
	3.1.1 Mesh size independence

Table 5 .

 5 1: Nominal chemical composition (weight %) from[START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF]

	C	Mn	Si	P	S	Cu	Al	Fe
	0.08	1.6	0.4	<	<	< 0.2 < 1.	bal.
				0.05	0.01

Table 5 .

 5 2: Table summarizing the number of tests per configuration realized in the experimental campaign at room temperature of 20 • C. Tests also performed at 100 • C are market with a red dot •. Those used for the plasticity constitutive model identification are marked with , and those for the damage constitutive model are market with . (from

Table 5 .

 5 3: Set of identified parameter from[START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF]

	3.1 Error estimator on thin metallic sheets
	1
	2	(5.10)

Let us recall from chapter 2, the scalar incremental error estimator based on the Super convergent Patch Recovery (SPR) adapted for non linear problems: η * , incr scalar = ||e * ,incr scalar || E * where E * = Ω f * .∆f * dΩ

Table 5 .

 5 4: Continuous-discontinuous parameters used in the simulations presented in 4

	Error estimator and remeshing
	parameters
	Minimal element size	h min	l b /2 = 200µm
	Threshold error	ε 0	0.01
	Remeshing frequency	∆t remeshing	adjustable
	Marching ridges algorithm
	Evaluation radius	R	2l b = 800µm
	angular precision	∆θ	5 •
	Insertion criterion
	Critical value for insertion	f crit	0.8f r = 0.16
	Crack increment length	∆a	2l b = 800µm
	Reequilibrium time increment	∆t reeq	1.e -4s

  The outputs variable are expressed as function of the state variables only so that:The matrices ∂R/∂∆V IN and ∂∆V OUT /∆∂V S are computed as block-matrices. The calculation is straightforward with:

												One therefore
	has:	δR = δ0 =	∂R ∂∆V IN	.δV IN +	∂R ∂∆V S	.δV S =	∂R ∂∆V IN	.δV IN + J .δV S	(B.1)
	so that					δV S = -J -1 .	∂R ∂∆V IN	.δV IN	(B.2)
		δV OUT =	∂∆V OUT ∂∆V S	.δV S = -	∂∆V OUT ∂∆V S	.J -1 .	∂R ∂∆V IN	.δV IN	(B.3)
	so that the consistent tangent matrix is equal to:
													∂∆σ	∂∆σ	∂∆σ	
		K mat = -	∂∆V OUT ∂∆V S	.J -1 .	∂R ∂∆V IN	=	        	∂∆ε ∂∆ω ∂∆ε ∂∆κ	∂∆ω ∂∆ω ∂∆ω ∂∆κ	∂∆κ ∂∆ω ∂∆κ ∂∆κ	        	(B.4)
												∂∆ε	∂∆ω	∂∆κ
		∂R e ∂∆ε	= -I,		∂R g ∂∆ω	= -(1 -f ),	∂R n ∂∆κ	= -A n	(B.5)
	and				∂∆σ ∂∆ε e	= E,	∂∆κ ∂∆κ	= 1,	∂∆ω ∂∆ω	= 1	(B.6)

(1 -κ) in McClintock analysis.

In case of heterogeneous media, σ * can be constructed for each homogeneous part of the structure as will be described in the next section, while adding a constraint: imposing that the normal stress across the interfaces between different materials remain continuous[El Ouazani Tuhami et al.,

2018].2 p designates the degree of the mesh shape functions and h the fineness of the mesh

i.e. one must determine the element of the old mesh containing the considered node

The model is written using the updated Langrangian framework, in order to obtain an accurate description

Summary

Note that the model can be easily extended to plastically anisotropic materials by replacing the von Mises by any stress measure accounting for anisotropy[Benzerga and Besson,

[START_REF] Hambli | Finite element model fracture prediction during sheet-metal blanking processes[END_REF][START_REF] Tanguy | Plastic and damage behavior of a high strength X100 pipeline steel: experiments and modelling[END_REF][START_REF] Shinohara | Anisotropic ductile failure of a high-strength line pipe steel[END_REF].

Continuous-discontinuous approachIn this section, the continuous-discontinuous approach described in chapter 2.27 is used to simulate the quasi-static part of the experimental database, in order to evaluate the predictive capability of this approach. This approach is very similar to the one used in the chapter

4, but with a different remeshing strategy. In the present case, remeshing is based on the data provided by an incremental error estimator. This choice was made for cost reasons. Indeed for the considered experimental database, the active process zone is relatively large and finely meshing this whole region with the minimal mesh size would be too expensive, especially for 3D simulations. On the contrary, error estimation coupled with an appropriate remeshing criterion enables to obtain good accuracy-cost compromise.

In the present case, the total porosity is the field chosen for the error estimation

Chapter 5

Application to an existing experimental database Contents 1

Inputs from to the previous work of [START_REF] Davaze | Numerical modelling of crack initiation and propagation in ductile metallic sheets for crash simulations[END_REF][START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] 1 with the evolution of each damage mechanism given as: ḟn = A n (κ) κ (5.5) ḟg = (1 -f ) ω (5.6) ḟsh = k w f (1 -L(σ) 2 ) σ dev : εp σ eq (5.7) In these equations, k w is a parameter to identify and L is the Lode parameter as defined in 2.28.

Strain rate dependence

In the work of [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF], the strain rate effects are modeled using the following flow expression with the Johnson-Cook model for the viscous part and the Zhao for the thermal softening part:

Thermal softening

Strain-rate effects (5.8) Because only the quasi-static applications are considered in this chapter, the thermal softening effect can be neglected. Therefore, another approach is chosen to model the strain rate dependence:.

where ṗ0 and σ 0 are two material parameters. And σ is the GTN effective stress as defined in 2.26.

Model parameters

The identified parameters from [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF] are summarized in table 5.3. The viscosity parameters σ 0 , n and ṗ0 are selected to correctly fit the quasi-static part of the data base. Note that the identified characteristic length from [START_REF] Davaze | A non-local damage approach compatible with dynamic explicit simulations and parallel computing[END_REF] is attributed to both l κ and l ω .

Appendix A

Implementation of the nonlocal GTN model

Integrating the set of equations 3.16 to 3.20 over a finite time step ∆t using a fully implicit scheme is equivalent to solving the following set of non linear equations with respect to the increments of the state variables ∆V S = (∆ε e , ∆κ, ∆f g , ∆f n , ∆ω) for a given increment of the input variables ∆V IN = (∆ε, ∆ω, ∆κ) :

Using a fully implicit scheme all variables in the previous system are evaluated at the end of the time increment. Solving the system requires the evaluation of its Jabobian matrix which is formally expressed as: .6) where R = (R e , R κ , R g , R n , R ω ). The Jacobian matrix can be computed block-wise as follows (the zero terms are omitted).

• Derivatives of R e :

)

• Derivatives of R g :

Appendix C

Global problem discretization 1 Spatial discretization

Elements have nodal DOFs corresponding to the displacements and nonlocal variables κ and ω. They are represented as vector u e , κ e and ω e . Linear shape functions are used to interpolate the nonlocal nodal variables in the elements (i.e. κ e and ω e ), so that:

where N is a matrix formed with the linear shape functions. The gradients of κ and ω are computed as:

where the matrix G is formed with the derivatives of the linear shape functions with respect to the final configuration. The gradient with respect to the initial configuration is computed as ∇ 0 κ = G 0 .κ e where G 0 is formed with the derivatives of the linear shape functions with respect to the initial configuration.

Standard quadratic shape functions are used to interpolate the nodal variables associated with the displacements in the elements (i.e. u e ). The transformation tensor F and the velocity gradient tensor L are computed as:

where B F is a matrix formed with the derivatives of the quadratic shape functions with respect to the initial configuration whereas B L uses the derivatives of the same shape functions with respect to the final configuration. The deformation rate D is expressed as B D . ue .

Temporal discretization

The finite strain formulation is obtained using a mid-point integration scheme over the time step [t, t + ∆t]. The transformation tensor at t + 1 2 ∆t is first computed as

2 ∆F using the standard polar decomposition. The tranformation increment is then computed as:

where ∆D and ∆W are respectively the symmetric and skew-symmetric parts of ∆L. The rotation Q is updated as (based on the mid-point integration of eq. 3.22)

Proof: Eq. 3.22 is: Q = W .Q. The mid-point integration corresponds to:

Solving for ∆Q leads to eq. C.5

The strain ∆ε is updated as:

The constitutive equations are then integrated (see A) using this strain increment. The resulting stress tensor is rotated back as

. The 2 nd (e.g. ∂∆σ/∂∆ω) and 4 th (∂∆σ/∂∆ε) order tensors in the consistent matrix (eq. B.4) are also rotated back using Q 1 2 . The rotated quantities are indicated with a | . In particular ∂∆σ/∂∆ε| = C J is associated with the Jaumann rate of the stress tensor Σ. The tangent operator corresponding to the Truesdell rate is then given by:

Details of the terms involved in the global system 1 Elementary reactions

From the discretized form of the weak formulation (3.25 and 3.26), elementary reactions associated unknowns are given as:

where notation {Σ} indicates that the tensor Σ is expressed using Voigt notations. Let us underline that the integrals are taken over the current configuration. Integration is performed using standard Gauss integration. Let us introduce here the function M R such that the product c of two second order tensors a and b, usually written c = a.b, can be expressed using Voigt notations as:

The matrix M R (b) thus depends linearly on b.

Elementary stiffness matrix

It is then necessary to evaluate the elementary stiffness matrix which is computed as a block matrix.

Calculation of the first bloc line of the elementary stiffness matrix can be obtained through the derivative of F u with respect to u e , κ e and ω e :

Derivation of Ḟ u leads to:

In the central term, the derivative tensor Σ can be separated into three parts corresponding to variations relative to the different unknowns:

R ÉSUM É

La pr édiction de l'endommagement sous des chargements accidentels constitue un enjeu majeur pour garantir la s écurit é des installations industrielles : centrales nucl éaires, pipelines, etc. Il est alors important de pr édire correctement l'initiation et la propagation des fissures. L'objectif de ce travail est de d évelopper une strat égie fiable et efficace pour pr édire l'initiation et la propagation des fissures sur de longues distances en utilisant la m éthode des él éments finis (FEM) dans le cadre d'une approche locale. La premi ère partie de la th èse consiste à d évelopper un mod èle d'endommagement ind épendant du maillage et adapt é aux m étaux tels que les aciers de construction. Pour ce faire, un mod èle GTN non local à gradient implicite avec deux longueur internes a ét é d évelopp é et impl ément é dans le logiciel Z-set. Ce mod èle est adapt é pour repr ésenter les deux m écanismes d'endommagement actifs dans les aciers : (i) la nucleation, (ii) la croissance des cavit és. Cependant, l'utilisation de mod èles non locaux est limit ée par le co ût de calcul, en cas des discr étisations fines d'un long chemin de fissure. Ce probl ème est abord é en utilisant des techniques de remaillage. Les m êmes courbes de contrainted éformation et les m êmes trajets de fissure ont ét é obtenus en utilisant le remaillage par rapport à la simulation à maillage fixe avec 20 fois moins de DDLs. En outre, une approche continue-discontinue est adopt ée pour repr ésenter avec pr écision la cin ématique de la propagation de la fissure. Cette nouvelle approche est appliqu ée à une base de donn ées exp érimentale existante, et un bon accord entre les simulations et les exp ériences est trouv é.

MOTS CL ÉS

Fissuration ductile, mod èle GTN non local, gradient implicite, techniques de remaillage automatique, insertion de fissures discr ètes ABSTRACT The prediction of damage development under hypothetical accidental loading constitutes a major stake to guarantee the safety of industrial installations: nuclear power-plants, pipelines, etc. It is then important to correctly predict crack initiation and propagation. The aim of this work is to develop and implement a reliable and efficient strategy to predict crack initiation and propagation over long distances using the Finite element Method (FEM) in the scope of local approach ro failure. The first part of the thesis consists in developing a mesh objective GTN damage model adapted for metals such as construction steels. To do so an implicit gradient nonlocal model with two internal length scales was developed and implemented in the Zset software. The two internal length scales model is tailored to represent both damage mechanisms active in steels: (i) void nucleation on strengthening particles, (ii) void growth. However, the use of nonlocal models is limited by the computational cost, especially for long distance propagation computations, due to the fine discretization needed in the localization band. This problem is addressed by using mesh adaption. Same stress-strain curves and crack-paths were obtained using mesh adaption in comparison to fixed mesh simulation with 20 times fewer DOFs. In addition, a continuous-discontinuous approach is adopted to accurately represent the kinematics of the crack propagation (i.e. displacement jumps). This new approach is applied to an existing experimental database, and good agreement between simulations and experiments is found.

KEYWORDS

Ductile failure, nonlocal GTN model, implicit gradient, automatic mesh adaption, discrete crack insertion