Balaji, Lukasz, Eduard, Mijo, Benoît, David Laurine Lydie

Céline

Anaïs

Pierre Anthony

Liang, Jacky, Sitou, Mouldi... Un Autre Pierre

Matt, Fanny Ben Thibault

Noar, Mia Un Claire

Mica Autre Alex

Emilie, Dimi, Solveig, Guillaume, Sylvain, Sophie, Fabien, Lol, Naza... Un Autre Thibault

Keywords: viscosity procedure, Computational results for different viscosity parameters

Après ces années de dur labeur, l'heure des remerciements est enfin venue. La liste des personnes qui m'ont soutenues durant cette période riche en émotions, en échanges scientifiques, en tuyaus techniques variés, est bien longue pour en faire une synthèse exhaustive.

Je commence par remercier les membres du jury qui ont bien voulu participer à l'examen et au rapport de cette thèse, Messieurs les rapporteurs Carl Labergère et Luc Davenne, ainsi que les examinateurs, les Professeurs Adnan Ibrahimbegović et Jean-Michel Bergheau avec qui les échanges d'après soutenance ont été très riches scientifiquement.

Je continue avec les personnes qui ont suivi de très près ce travail, mon directeur de thèse Jean-Marc Roelandt avec qui les échanges scientifiques ont été très enrichissants et se sont toujours fait dans une convivialité très plaisante. J'adresse un remerciement tout particulier à ma directrice de thèse Delphine Brancherie dont les qualités humaines, les connaissances scientifiques ainsi que les compétences techniques ne sont pas étrangères à l'aboutissement de ce travail. Son soutien moral pendant les longues heures/journées/semaines/mois de debugging a été irréprochable et indispensable pour me faire perséverer dans ces moments difficiles. Merci Delphine. En espérant pouvoir travailler à l'avenir avec toi dans d'aussi bonnes conditions que durant ma thèse.

Merci également à Dominique Deloison ainsi qu'à Fabrice Congourdeau d'Airbus Group pour les échanges que l'on a pu avoir et qui ont apporté un regard industriel à ce travail, ce qui a notamment permis d'aborder certaines problématiques avec une approche différente de l'approche académique, apportant donc une dimension supplémentaire à ces travaux. Je remercie également Pierre Feissel pour son expertise technique linuxienne de premier choix, qui m'a notamment permis d'éviter de perdre la moitié de mes données au milieu de ma thèse.

Merci aux différents collègues de labo pour les nombreux moments conviviaux passés: Ludo,

Introduction xv

Chapter 1 Ductile fracture modeling: a review [START_REF] De Souza Neto | A fast, one-equation integration algorithm for the Lemaitre ductile damage model[END_REF]

Introduction

Understanding ductile fracture phenomena has been essential to the mechanical engineering community for the past forty years, especially in the context of applications such as the design of aeronautical and space structures, or metal forming processes.

Existing numerical methods, such as standard Finite Element methods, tend to encounter difficulties to objectively predict strain localizations and crack propagations involved in the ductile fracture process, since the material models usually entailed when modeling ductile fracture include a softening behavior. This specific behavior is of crucial importance to account for the loss of carrying capacity inevitably involved when modeling structure responses until fracture.

Indeed, it is now a well-known fact that, in a Finite Element modeling context, the "post-peak" phase of a structure response subjected to severe loadings until fracture are crucially depending on the chosen Finite Element discretization of the structure. The element size as well as its type and orientation are critical parameters. This is essentially due to the fact that standard Finite Element methods are initially designed to model volumic phenomena and not surfacic phenomena such as fracture. Modeling ductile fracture actually requires the introduction of an additional parameter which is a characteristic length that corresponds to the domain in which fracture occurs. In a standard Finite Element context, this characteristic length automatically becomes the element characteristic size, unavoidably inducing a crucial dependence of the response on this element size. These considerations will be dealt with in details in the following chapter in section 1.3.1.

For the past forty years, the fracture mechanics community has been adressing this problem by proposing novel methods. These methods all have advantages and drawbacks as we will see in section 1. [START_REF] Armero | New finite elements with embedded strong discontinuities in the finite deformation Range[END_REF]. A criterion that we use to choose and enhance one of these methods is the amount of programming that this enhancement involves at the different stages of the standard Finite Element method.

Another reason for this choice is the fact that, in the short term, one of the objective of this work is the possible implementation in an industrial context. Since the vast majority of industries that rely on numerical simulation are using industrial Finite Element codes, and since for most of these codes the programmer does not have access to its whole architecture, it becomes critical to limit the enhancements involved in the method to the levels of the architecture that the programmer has access to. Mainly, these levels are the element level and the material point level. This means that we will limit our enhancements to these levels of the FE architecture.

Having considered these requirements, one of the frameworks that fulfills these conditions is the Strong Discontinuity Approach (SDA) framework whose foundings and basics will be discussed in section 1.3.4. As specified above, this method has the advantage to only involve local enhancements. Indeed, the corresponding numerical enhancements remain limited to the elements involved in the localization zone, each enhancement being independent from an element to another. This latter consideration is of great importance since we will see that it allows the non-linear Finite Element resolution procedure to remain unchanged from standard non-linear computations. In particular we will see that this is due to the fact that no additional degrees of freedom are considered at the global level of the code architecture, contrary to X-FEM based methods for instance. However this approach has drawbacks, for instance its inability to ensure the continuity of the crack opening and sliding from one element to another. This framework has mostly been the subject of brittle and quasi-brittle fracture modeling, and few authors have intended to adapt this method to the ductile fracture context. One of the main contributions of this work is to propose a model that accounts for coupled plasticity and damage in the context of the Strong Discontinuity Analysis, both of these dissipative mechanisms being of great importance when modeling ductile fracture, as we will see in the following chapter in sections 1.1 and 1.2.

The structure of the manuscript is as follows. In a first chapter, the physical mechanisms involved in ductile fracture are discussed, as well as the different techniques that have been developed to model these phenomena, without any ambition of being exhaustive, but rather with the goal of emphasizing those that have most influenced the ductile fracture scientific community. In this chapter, the strain localization phenomenon is also highlighted from a numerical point of view, along with the different techniques that have been developed to circumvent this phenomenon, again without any claim of proposing an exhaustive review, since the amount of methods is too important to be mentionned in this manuscript.

In a second chapter, the SDA theoretical formulation and its numerical implementation are described in the context of Lemaitre's phenomenological coupled damage-plasticity model (see [START_REF] Lemaitre | A continuous damage mechanics model for ductile fracture[END_REF]). A four-node quadrangular element enhanced with a linear displacement jump is described, and numerical examples showing the capabilities of the method are presented. A critical damage value is used as a criterion to introduce the discontinuity surface. All numerical developments presented have been implemented in FEAP, an open source Finite Element software developed by R.L. Taylor at the beginning of the 90's (see [START_REF] Taylor | FEAP website[END_REF] for FEAP's website), relying on the books written by O.C. Zienkiewicz and R.L. Taylor [START_REF] Zienkiewicz | The Finite Element Method: Its Basis and Fundamentals[END_REF][START_REF] Zienkiewicz | The Finite Element Method For Solid and Structural Mechanics[END_REF].

xvi In a third and final chapter, we aim at extending the formulation proposed in the previous chapter to the large strain hypothesis. This is one of the main objectives of this work since ductile fracture often involves large strain scenarios, thus the limitation to the small strain hypothesis becoming irrelevant. In this chapter and in particular, we will see that the corresponding numerical enhancements do not involve considerable changes in the structure of the code from the small strain hypothesis case, even if the hypothesis itself is significantly different.

Finally we will conclude and suggest a few perspectives with respect to the present work.

Chapter 1

Ductile fracture modeling: a review Ductile fracture modeling has been the subject of many different approaches, mainly phenomenological and micromechanical approaches. Without the ambition of being exhaustive, these different techniques are reviewed in this chapter. The strain localization phenomenon is also highlighted from a numerical point of view, along with the different techniques that have been developed to circumvent this phenomenon

Ductile fracture physical mechanisms

Ductile fracture is composed of three distinct stages. The first stage deals with void nucleation in which voids are initiated at material defects, that is voids or inclusions (in most cases) related to the material microstructure. In a second phase, voids grow due to plastic deformation of the matrix surrounding the void, with a growing dependent on the stress triaxiality ratio. At high stress triaxiality values, voids grow until they coalesce to form a micro-crack, which means that several voids merge.

This situation is prevalent when only one population of inclusions is present in the microstructure as it was shown by several authors [START_REF] Bron | Déchirure ductile des tôles minces en alliage d'aluminium 2024 pour application aéronautique[END_REF][START_REF] Tanguy | Modélisation de l'essai charpy par l'approche locale de la rupture : application au cas de l'acier 16MND5 dans le domaine de transition[END_REF]. These authors have found that with these alloys that contain only one population of inclusions, and at high stress triaxiality a macro-crack appears due to the coalescence of primary voids that are nucleated at inclusions. This phenomenon is reffered to as internal necking.

However, when the microstructure contains several populations of inclusions at different length scales, another fracture mechanism may prevail at low stress triaxialities. This physical process was reffered to as void sheeting by Garrison et al. [START_REF] Garrison | Ductile fracture[END_REF]. In this case, primary voids remain small due to low stress triaxiality and secondary voids are nucleated at narrow strain localization bands to form a macro-crack.

In most industrial cases and considering highly sophisticated metallic alloys, the actual fracture process is a combination of these two mechanisms, which depends on several material and loading characteristics such as the types of voids, inclusions, their spacing, their sizes, the stress triaxiality ratio, the plastic strain ...etc... In a more phenomenological description, ductile fracture is composed of three phases: an elastic phase, a pre-peak phase in which plasticity and damage is increasing until a global peak is reached, and a post-peak phase in which strain localization occurs, as well as the formation of a macro-crack leading to its propagation until complete fracture.

Models for ductile fracture

For the past thirty years, authors have intended to describe and model these mechanisms by different approaches, in order to be able to implement these models in industrial numerical tools. Very few of them have intended or succeeded in modeling the whole mechanism, that is the three different stages of the ductile fracture process. They have often focused on modeling void growth. This large amount of studying and modeling these processes has led to the rise of new challenges in numerical modeling which we will intend to describe in the following section, after having presented the main different categories of ductile fracture modeling.

Many approaches have been used to model the ductile fracture process. As hard as it is to classify these approaches, we may divide them into two main categories which rely on the type of description that is used, whether it is a micromechanical description (with the three stages involving void nucleation, growth and coalescence) or a phenomenological description relying on the effects of the mechanisms at the macro-scale (elasticity, then plasticity and damage, and finally strain localization and fracture).

The first category deals with micromechanical models in which the microscopic mechanisms are modeled. They mainly rely on Gurson's famous approach on void growth micromechanical modeling. The second category deals with phenomenological models (or empirical models), in which models are established by considering the macroscopic phenomena involved in ductile fracture. This category contains various types of models that are always based on macroscopic empirical considerations. They can be derived from a rigorous thermodynamical framework or not.

However, these categories are not fully closed and some models may belong to both categories, besides it is possible to categorize these models using a different heuristic. For example, Li and Wierzbicki [START_REF] Li | Prediction of plane strain fracture of AHSS sheets with post-initiation softening[END_REF] proposed to categorize ductile fracture models with respect to their approach of fracture, whether fracture is considered to occur spontaneously when some critical situation is reached with no dependence on the loading history, or whether it considers damage as a cumulative process, this time depending on the loading history. In this situation, fracture occurs when this cumulative damage measure reaches a critical value.

Micromechanical models

The micromechanical approach relies on the study at the micro-scale (the scale of the void responsible for fracture) of the physical process that leads to fracture. The first micromechanical developments for modeling ductile fracture were led by McClintock [START_REF] Mcclintock | A criterion for ductile fracture by the growth of holes[END_REF], and Rice and Tracey [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF]. They isolated a cylindrical void and a spherical void, respectively, and described its growing in a rigid perfectly plastic matrix. They emphasized the important role of the plastic strain and the triaxiality ratio in their modeling. A simplified version of the void growth model for high stress triaxialities is expressed as:

Ṙ R 0 = αexp 3 2 σ h σ y εeq (1.1)
where R is the current void radius, R 0 is the initial void radius, σ h is the hydrostatic pressure, σ y is the matrix yield stress and ε eq is the von Mises equivalent strain. α initially set to 0.283 was then modified by Huang [START_REF] Huang | Accurate dilatation rates for spherical voids in triaxial stress fields[END_REF] for better accuracy. Other authors have extended this model to account for strain hardening effects on ductility (see [START_REF] Budiansky | Void growth and collapse in viscous solids[END_REF][START_REF] Becker | The effect of void shape on the development of damage and fracture in plane-strain tension[END_REF]). By integrating equation (1.1) and assuming a critical void radius above which fracture occurs, this model has led to a simple criterion [START_REF] Marini | Experimental study of cavity growth in ductile rupture[END_REF]:

R R 0 = R R 0 c (1.2)
This criterion does not take into account the three stages of the ductile fracture process, even though the criterion expressed above may be interpreted as a critical void growth above which coalescence occurs.

Besides, the model does not take into account the interactions between voids which implies that coalescence is not represented in the model, as for the void nucleation stage. Another critic is that this void growth has no influence on the material behavior, and no softening can occur with such a modeling.

Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: part I-yield criteria and flow rules for porous ductile media[END_REF] proposed a model that takes this critic into account. He proposed that the distribution of voids in the matrix is idealized as a unique spherical void containing all the volume of void, and surrounded by a rigid perfectly plastic sphere. He performed an upper bound analysis of this model to obtain the yielding condition expressed below in (1.3). The model uses the void volume fraction f which is equivalent to porosity, and represents the ratio of the volume of void over the volume of the sphere. In that case, the growing void has an effect on the material behavior through the void volume fraction that is coupled with the yielding function, originally expessed by Gurson as:

Φ = σ eq σ y 2 + 2f cosh 1 2 σ h σ y -1 -f 2 (1.3)
where σ eq is the von Mises equivalent stress. It has to be noticed that without any void, the Gurson model is exactly the same as the conventional von Mises model, and that damage increases until there is a complete loss of load-carrying capacity only when f = 1.

Since this is an unrealistic behavior because the load-carrying capacity is completely lost only when there is no material anymore, Tvergaard [START_REF] Tvergaard | Influence of voids on shear band instabilities under plane strain conditions[END_REF][START_REF] Tvergaard | On localization in ductile materials containing spherical voids[END_REF] introduced new parameters q 1 and q 2 in his modified version of the Gurson model, expressed as:

Φ = σ eq σ y 2 + 2q 1 f cosh 3 2 q 2 σ h σ y -1 -(q 1 f ) 2 (1.4)
This model is often reffered to as the Gurson-Tvergaard-Needleman model (GTN). By setting q 1 = 1.5 and q 2 = 1 based on the bifurcation analysis of Gurson's model and his numerical studies for a material containing periodic distributions of voids, he obtained a loss of load-carrying capacity at f = 1/q 1 . This value is still too large based on experimental observations carried out by Brown and Embury [START_REF] Brown | The initiation and growth of voids at second phase particles[END_REF] and on unit cell analysis performed by Koplik and Needleman [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF]. They showed that the critical volume fraction above which coalescence occurs is usually less than 15%.

Faleskog et al. [START_REF] Faleskog | Cell model for nonlinear fracture analysis -I[END_REF] showed that q 1 and q 2 depend on the plastic hardening exponent and on the ratio of the yield stress over Young's modulus. Koplik and Needleman highlighted the fact that Gurson's model was not able to predict void coalescence and that a specific criterion should be introduced along with this void growth model.

Consequently, authors have used a so-called critical void volume fraction for calescence in the literature. It simply assumes that coalescence occurs when a critical volume fraction f c is reached, without assuming any stress state dependency or any initial void volume fraction. This was assumed to be a material constant that can be determined experimentally or numerically, a common value is f c = 0.15 [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] or f c = 0.1 [START_REF] Xia | A computational approach to ductile crack growth under large scale yielding conditions[END_REF][START_REF] Ruggieri | Numerical modeling of ductile crack growth in 3-D using computational cell elements[END_REF].

Tvergaard then addressed the problem that the critical void volume fraction does not depend on the initial void volume fraction, and therefore it depends on the material and can not be a simple constant.

Koplik and Needleman also commented on the influence of stress triaxiality on f c and concluded that if the initial void volume fraction was small (ie < 1%), then the influence of stress triaxiality could be neglected. However for larger values, f c at high stress triaxiality was significantly lower than for the case at low stress triaxiality.

Tvergaard and Needleman then proposed to simulate coalescence with the following function:

f * = f for f f c f c + f -fc f R -fc ( 1 q 1 -f c ) for f > f c (1.5)
where f R is the fracture porosity that can be freely chosen. f * can be seen as an effective porosity. This function simply states that if the void volume fraction is below the coalescence criterion, then the void volume fraction follows Gurson's standard void growing rule, and if it reaches a critical value, the void growing is amplified to represent the sudden loss of load-carrying capacity. Figure 1.1 was established by Zhang [START_REF] Zhang | A complete Gurson model[END_REF] and displays the evolution of the effective void volume fraction along with its influence on the load-carrying capacity within a FE computation.

Besson highlighted the fact that any other function of f * could be suitable [START_REF] Besson | Continuum models of ductile fracture: a review[END_REF]. and that low values of f R and f c can lead to convergence problems. It is to be noted that there is a large difference of void volume fraction between the begining and the end of coalescence and that it is un unstable and very quick process, whereas void growth is more stable and produces much less void. Koplik and Needleman [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF] and Zhang [START_REF] Zhang | A complete Gurson model[END_REF] emphasized the accuracy of the function in (1.5) with respect to void coalescence simulation although it is proven to be less accurate for high stress triaxiality values.

Nucleation in a micromechanical-based analysis is much less studied than void growth or coalescence, even though its role in the process of ductile fracture is [START_REF] Zhang | A complete Gurson model[END_REF]) undoubtedly significant. This aspect of ductile fracture is the less known stage of the process and authors have reported that it is a highly material-dependent phenomenon. Indeed, nucleation of voids depends on the distribution of inclusions, their sizes and shapes, their strength, and on the hardening exponent of the plastic law of the matrix. It is often studied in a phenomenological way in the following form:

ḟn = A n ṗ (1.6)
where ḟn is the change of porosity due to nucleation, A n is the strain driven nucleation rate and p the accumulated plastic strain. Then the total change of porosity is expressed as:

ḟ = ḟn + ḟg + ḟc (1.7)
where ḟg and ḟc are the change of porosity due to void growth and coalescence, respectively.

Chu and Needleman proposed a model [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF] of this strain driven nucleation rate that is often coupled with the GTN model. It is expressed as a Gaussian function as:

Models for ductile fracture

A n = f N s N √ 2π exp - 1 2 p -ε N s N 2 (1.8)
where f N is the volume fraction of inclusions at which voids can be nucleated, ε N is the strain for which 50% of the inclusions are broken, s N is the standard deviation of the nucleation strain. Chu and Needleman [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF] calibrated their model to find ε N = 0.3 and s N , and many authors have used the same values without taking into account that they highly depend on the material microstructural properties. Other authors have proposed other forms for A n (e.g. [START_REF] Zhang | A complete Gurson model approach for ductile fracture[END_REF][START_REF] Besson | Modeling of scatter and size effect in ductile fracture: application to thermal embrittlement of duplex stainless steels[END_REF]).

Although it is not often used, a similar expression can be obtained for a stress controlled nucleation intensity, and the total void nucleation rate is then written as:

ḟn = A n ṗ + B n ( σh + σeq ) (1.9)
This allows the hydrostatic stress in equation (1.9) to have an effect on nucleation which has a physical sense. This law is a statistical void nucleation law.

Zhang [START_REF] Zhang | A complete Gurson model[END_REF] highlighted that there exist two other nucleation mechanisms, the cluster of voids nucleation mechanism in which voids are nucleated when some critical condition has been reached, and it is often assumed that voids nucleate at the beginning of plastic deformation, setting a so-called initial void volume fraction f 0 .

The second other type of nucleation model is called a continuous nucleation model, in which it is considered that the number of voids increases linearly with plastic deformation as it was shown by Gurland [START_REF] Gurland | Observation on the fracture of cementie particles in spheroidized 1.05 C steel deformed at room temperature[END_REF] for some specific material. Zhang in [START_REF] Zhang | A complete Gurson model approach for ductile fracture[END_REF] did not recommend to use Chu and Needleman's statistical model due to the difficulty that arises from the determination of the three parameters of the model, and due to the fact that it may lead to non-uniqueness issues, in a numerical point of view.

The model proposed by Chu and Needleman is actually a statistical version of the cluster of voids model, it states that most of the void will be nucleated at a critical value of plastic deformation for the strain-driven void nucleation model.

Gurson's model and low stress triaxiality handling

When stress triaxialities are low, for which shear stress becomes prevalent with regards to hydrostatic stress (η < 0.33), it is proven that Gurson's original model is not efficient to simulate the so-called void sheeting mechanism (see section 1.2.2). Shear damage is actually caused by void rotation, when shear deformation increases, voids are driven towards each other and the ligaments between voids reduce which induces a loss of carrying capacity as it does when necking occurs at high stress triaxialities. This mechanism prevails in the context of metal forming in which stress triaxialities remain low.

Chapter 1. Ductile fracture modeling: a review

Recent works (e.g. [START_REF] Bao | On fracture locus in the equivalent strain and stress triaxiality space[END_REF]) have shown that void coalescence at low stress triaxialities is related to the critical strain, which can be approximated with a function of the stress triaxiality and the amount of shear loading. There exist shear damage models that require the precise information of the microstructure geometry. For instance Bacha's model [START_REF] Bacha | Metal ductility at low stress triaxiality application to sheet trimming[END_REF] requires lots of geometrical parameters to define the microstructure which induces a prohibitive calibration method, considering that this calibration is even possible. Model parameters should remain as low as possible so an identification procedure remains possible in an industrialization perspective of the method. Jackiewicz proposed a review [START_REF] Jackiewicz | Use of a modified Gurson model approach for the simulation of ductile fracture by growth and coalescence of microvoids under low, medium and high stress triaxiality loadings[END_REF] of models that incorporate the shear-driven coalescence mechanism, and which do not require a lot of parameters. In the following a quick review of these approaches is proposed.

Mc Veigh et al. proposed to include the shear-driven coalescence mechanism considering the assumption stated above, namely the relation between this mechanism with a critical strain. He considered that the simplest way to include the shear mechanism was to include a dependence of the porosity rate on the equivalent plastic strain rate, since it describes quite simply the amount of shear deformation rate. He proposed the following expression for the porosity rate, to be compared with equation (1.7):

ḟ = ḟn + ḟg + ḟc + K ṗ (1.10)
The term K ṗ does not have a physical sense but it enables Gurson's model to account for the shear-driven coalescence mechanism that prevails at low stress triaxialities.

Nahshon and Hutchinson proposed a slightly different approach, considering that the variable f is not only related to void growth, but should be interpreted as a universal damage measure, which describes the loss of carrying-capacity in a phenomenological way. It was emphasized that in a pure shear stress state, no void growth occurs with the actual definition of Gurson's model, however there is a loss of load-carrying capacity due to void distortion and ligament reduction between voids.

The authors proposed to add a term in the definition of ḟ to account for damage increase in pure shear stress states, without modifying the accurate definition of this rate in the case of axisymmetric stress states. They expressed the damage rate as follows, in the finite deformation framework, without accounting for damage increase due to nucleation:

ḟ = (1 -f )tr( Ėpl ) + k ω f ω(S) S dev ij Ėpl ij S eq (1.11)
where E pl is the plastic part of Green-Lagrange's deformation tensor, S is the 1.2. Models for ductile fracture the second Piola-Khirchoff stress tensor, k ω is the unique material parameter that sets the magnitude of the damage growth rate due to pure shear stress states, and ω(S) is expressed as:

ω(S) = 1 -[27J 3 (S)/2S 3 eq ] 2 (1.12)
where J 3 is a measure of the third invariant. This definition of ω(S) allows to not modify the axisymmetric cases since ω(S) = 0 in this case, and ω(S) = 1 for pure shear stress states, while it lies in the range 0 < ω(S) < 1 fo all intermediate stress states. In equation (1.11), the term (1 -f )tr( Ėpl ) refers to the damage increase rate due to void coalescence for the unmodified Gurson model.

Xue [START_REF] Xue | Constitutive modeling of void shearing effect in ductile fracture of porous materials[END_REF] performed analytical calculations on a unit cell containing a void and subjected to simple shear, in order to define damage increase due to shear. He analysed the reduction of the distance between the void free surface and the cell boundary during shear straining and defined a shear damage measure considering this assumption. He then obtained an expression of this shear damage rate as:

Ḋshear = qf n g θ E pl eq Ėpl eq (1.13)
where q and n are constants that only depend on the number of dimensions of the model (2D or 3D), and g θ is the Lode parameter dependent function, also related to the third invariant as in the model proposed by Nahshon and Hutchinson described above. This function is defined as:

g θ = 1 -6|θ L |/π and θ L = tan -1 1 √ 3 (2χ -1) (1.14)
where χ = (S 2 -S 3 )/(S 1 -S 3 ) and S 1 , S 2 and S 3 are the maximum, intermediate and minimum principal stress components.

Jackiewicz [START_REF] Jackiewicz | Use of a modified Gurson model approach for the simulation of ductile fracture by growth and coalescence of microvoids under low, medium and high stress triaxiality loadings[END_REF] in his review of these two additional terms proposed by Nahshon and Hutchinson on the one hand [START_REF] Nahshon | Modification of the Gurson Model for shear failure[END_REF], and by Xue on the other hand [START_REF] Xue | Constitutive modeling of void shearing effect in ductile fracture of porous materials[END_REF], highlighted the fact that the model proposed by Xue gives better accuracy of the strain localization, induced by the reduction of inter-particle spacing during large material rotations, than the one proposed by Nahshon and Hutchinson. In his paper, Jackiewicz also introduced a new simple coalescence criterion that states that coalescence occurs when the effective stress (that accounts for damage, i.e. section reduction) becomes singular.

Void shape and rotation influence

One of the major limitations of Gurson's model is that it cannot handle non-spherical voids which is a problem when low stress triaxiality cases are considered. It is possible with Gurson's model to tune parameters q 1 and q 2 by performing unit cell analysis corresponding to the actual microstructure. However, in order to describe more precisely this situation, Gologanu et al. [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-Case of axisymmetric prolate ellipsoidal cavities[END_REF][START_REF] Gologanu | Approximate models for ductile metals containing nonspherical voids-Case of axisymmetric oblate ellipsoidal cavities[END_REF] proposed to completely revise the analysis performed by Gurson by considering ellipsoidal voids (oblate and prolate) in their analysis. They obtained a closed-form expression for the yield surface for both oblate and prolate voids, loaded only axisymmetrically. It is reffered to as the GLD model. The only added parameter, in comparison with Gurson's model, is the shape aspect ratio of the void, namely the ratio of the radius in the direction of the symmetry axis over the radius in the radial direction.

The model is not expressed here but it has to be noticed that it becomes equivalent to Gurson's model when the aspect ratio is equal to 1. Besides, it enables to model flat voids and an equivalent porosity is computed that is equal to the porosity of a cell containing a spherical void with the same projected area as the 'penny-shaped' flat void. Finally, it is possible to consider cylindrical voids with this model, in the case of an infinite void aspect ratio.

As stated above, this model only allows axisymmetric loadings and therefore cannot handle void straining when the principal strain axes are not aligned with the ellipsoidal void axes, and therefore cannot induce void rotation and shearing. This means that this model does not properly fall in the category of models that can handle low stress triaxilities. [START_REF] Castañeda | Constitutive models for porous materials with evolving microstructure[END_REF], and Kailasam and Ponte Castañeda [START_REF] Kailasam | A general constitutive theory for linear and nonlinear particulate media with microstructure evolution[END_REF] proposed a model for 3D ellipsoidal pores in which loadings are not limited to the axisymmetric case contrary to the GLD model. In this model, the evolution of the aspect ratios of the cavities is deduced from the cavity mean deformation rate, which is itself obtained from the macroscopic plastic strain rate, using Eshelby theory [START_REF] Eshelby | The Determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF] considering plastic incompressibility. The formalism developed by Eshelby also provides the cavity mean rotation rate.

Ponte Castañeda and Zaidman

However, Besson [START_REF] Besson | Continuum models of ductile fracture: a review[END_REF] higlighted the fact that this formalism is probably more suitable for linear viscous solids than for plastic solids. Bordreuil [START_REF] Bordreuil | On modelling the growth and the orientation changes of ellipsoidal voids in a rigid plastic matrix[END_REF] considered this problem with a description of cavity rotation rate by using the representative theory developed by Wang [START_REF] Wang | A new representation theorem for isotropic functions: An answer to Professor G. F. Smith's criticism of my papers on representations for isotropic functions[END_REF] for skew-symmetric tensor-valued functions. He developed a model that contains parameters that need to be fitted by unit cell calculations. He showed that rotation rates strongly differ if the matrix is considered as a plastic or viscous solid. Before Gologanu's work on ellipsoidal void growth model [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-Case of axisymmetric prolate ellipsoidal cavities[END_REF][START_REF] Gologanu | Approximate models for ductile metals containing nonspherical voids-Case of axisymmetric oblate ellipsoidal cavities[END_REF], Thomason proposed a coalescence model using the same geometrical considerations [START_REF] Thomason | A theory for ductile fracture by internal necking of cavities[END_REF][START_REF] Thomason | Three-dimensional models for the plastic limit-loads at incipient failure of the intervoid matrix in ductile porous solids[END_REF], an axisymmetric unit cell containing an ellipsoidal void. He assumed that coalescence occurs when the inter-void ligament reaches its plastic limit load. Such an hypothesis tends to be confirmed by observations [START_REF] Benzerga | Coalescence-controlled anisotropic ductile fracture[END_REF] and unit cell calculations [START_REF] Koplik | Void growth and coalescence in porous plastic solids[END_REF]. The model uses a new material parameter, the cell aspect ratio L z /L x . He obtained 1.2. Models for ductile fracture the coalescence condition by writing the mechanical equilibrium of the ligament:

π(L 2 x -R 2 x )C f σ 0 = πL 2 x σ 33 (1.15)
where R x is the radius of the cavity in the radial direction, L x is the length of the cell in the radial direction, σ 0 is the elastic limit stress and C f is a constraint factor that describes the stress concentration in the ligament due to the cavity. C f was computed by Thomason by using the upper-bound theorem for limit-load analysis. Its approximation in the case of an axisymmetric cell is expressed as [START_REF] Thomason | Three-dimensional models for the plastic limit-loads at incipient failure of the intervoid matrix in ductile porous solids[END_REF]:

C f = 0.1 L x /R x -1 S 2 + 1.2 L x /R x (1.16)
Using Besson's notations and simplifications [START_REF] Besson | Continuum models of ductile fracture: a review[END_REF], the coalescence criterion is written as:

1 - R 2 x L 2 x C f = 2 3 + η X(f, η) (1.17)
where η is the macroscopic stress triaxiality ratio, X is defined as the ratio σ eq /σ 0 which depends on the void volume fraction f .

The evolution of the cell aspect ratio is related to the extensional plastic deformation through the relation:

(L z /L x ) (L z /L x ) 0 = exp 3 2 ε pl zz (1.18)
The main drawback of this criterion is that it doesn't account for void growth, and it is only based on geometrical considerations. The consequence is that X in (1.17) is overestimated. Microstructural parameters (void volume fraction f and void aspect ratio R z /R x ) should be computed using a dedicated model, such as Gurson's void growth model which was performed by [START_REF] Zhang | Analyzing ductile fracture using dual dilatational constitutive equations[END_REF][START_REF] Zhang | A new failure criterion for the Gurson-Tvergaard dilational constitutive model[END_REF]. In their works, void shape changes were neglected as it was assumed that voids remained spherical. Thomason's coalescence criterion was also used with the GLD model [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-Case of axisymmetric prolate ellipsoidal cavities[END_REF] to predict X, f and S = R z /R x [START_REF] Benzerga | Coalescence-controlled anisotropic ductile fracture[END_REF][START_REF] Pardoen | An extended model for void growth and coalescence[END_REF][START_REF] Pardoen | Micromechanics-based model for trends in toughness of ductile metals[END_REF] which is confirmed by unit cell calculations.

In their works Pardoen and Hutchinson [START_REF] Pardoen | An extended model for void growth and coalescence[END_REF][START_REF] Pardoen | Micromechanics-based model for trends in toughness of ductile metals[END_REF] combined Thomason's model for coalescence with the GLD model to model void growth. It consits of modeling two different modes of plastic deformation, the void growth model is dedicated to model diffuse plasticity and any other void growth model could be used (e.g. the GTN model), while the second one corresponds to localized plasticity in the intervoid ligament that occurs during coalescence. It means that both yield functions compete with each other, thus there exists a situation in which both yield functions intersect and plastic flow is defined by the cone of normals defined at the intersection point (see figure 1.2).

Φ = σ σ * 2 + 2q 1 f cosh q 2 2 σ h σ * -1 -(q 1 f ) 2 = 0 (1.19)
where σ * is the GTN effective stress and σ is an anisotropic stress measure. Various stress measures can be used. For example authors have often used the Hill model for plastic anisotropy (see e.g. [START_REF] Doege | Prediction of necking and wrinkling in sheetmetal forming[END_REF][START_REF] Rivalin | Ductile tearing of pipeline-steel wide plates: II. Modeling of in-plane crack propagation[END_REF][START_REF] Grange | Anisotropic behavior and rupture of hydrided ZIRCALOY-4 sheets[END_REF][START_REF] Brunet | Failure analysis of anisotropic sheet-metals using a non-local plastic damage model[END_REF]) or more advanced stress measures (see e.g. [START_REF] Barlat | A six-component yield function for anisotropic materials[END_REF][START_REF] Karafillis | A general anisotropic yield criterion using bounds and a transformation weighting tensor[END_REF]). Besson et al. have rigourously derived the Gurson-like yield equation in the case of Hill's anisotropic stress measure [START_REF] Benzerga | Plastic potentials for anisotropic porous solids[END_REF] in which voids are considered to remain spherical.

In cases of non-proportional or cycling loadings, it is important to account for kinematic hardening. Even though Gurson's micromechanical approach cannot account for the micromechanisms involved in a fatigue process, since its approach does not correspond to fatigue mechanisms, authors have intended to account for a predeformation or low-cyle fatigue within a Gurson-like approach (see e.g. [START_REF] Mear | Influence of yield surface curvature on flow localization in dilatant plasticity[END_REF][START_REF] Leblond | An improved Gurson-type model for hardenable ductile metals[END_REF][START_REF] Besson | An extension of the Green and Gurson models to kinematic hardening[END_REF]). However, experimental tests are still missing.

Calculations comparing kinematic hardening and isotropic hardening show that accounting for kinematic hardening accelerates the failure occurence. For all the models described above, it is assumed that once a void has nucleated, the particle that was the site of nucleation no longer influences the void growth or coalescence phenomenon. This is proven to be true as long as the stress triaxiality remains high so that void deformation is not influenced by particle interaction. However, at low stress triaxialities, voids tend to elongate in the direction of maximum deformation and shrink in the direction perpendicular to the direction of maximum deformation. Particles can prevent this shrinkage and thus increase porosity in this specific direction. A model based on this assumption and on the GLD model [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-Case of axisymmetric prolate ellipsoidal cavities[END_REF] was established by Siruguet et al. [START_REF] Siruguet | Effect of void locking by inclusions upon the plastic behavior of parous ductile solids-part II: theoretical modeling and numerical study of void coalescence[END_REF]. It allows to represent ductile fracture under pure shear loadings.

An alternative to Gurson's model

Alternatively to Gurson's approach on ductile fracture modeling on porous media, Cheng et al. [START_REF] Cheng | A stress-based variational model for ductile porous materials[END_REF] proposed a model applied to a rigid-plastic hollow sphere (as in Gurson's approach) that is based on a statical limit analysis combined with a homogenization technique in order to obtain a macroscopic yield criterion that accounts for the three invariants of the macroscopic stress field, while Gurson's model only takes into account the first two invariants of the stress tensor. In the authors' formulation, the internal boundary condition on the stress field is relaxed, leading to a so-called "quasi-lower bound" yield criterion. They obtained this "Stress Variational Model" by considering Hill's variational principal [START_REF] Hill | Mathematical Theory of Plasticity[END_REF] together with the introduction of a Lagrange multiplier to solve the minimization problem.

Phenomenological approaches

Contrary to micromechanical approaches to fracture, the phenomenological approach deals with empirical considerations. In this context, authors try to relate mechanical variables defined at the microscopic level to phenomena that occur at the macroscopic scale, such as the loss of carrying capacity in the context of ductile fracture. In this section and without the aim of proposing an exhaustive review, we first review the most famous fracture criteria that rely on a phenomenological approach, and then we review some of the most notorious models (mostly damage models) proposed by authors that also rely on phenomenological considerations.

Among the ductile fracture scientific community, in many works authors use in their models the notion of the so-called fracture strain, the strain at which fracture occurs. They postulate that fracture occurs when the accumulated plastic strain reaches a critical value, depending on a combination of the hydrostatic pressure, the von Mises equivalent stress and the principal tensile stress.

For all of these models, authors have intended to model the dependence of the fracture strain to different variables that characterizes the stress state, such as the hydrostatic pressure, the von Mises equivalent stress, the principal tensile stress, the principal shear stress, the third invariant of the stress tensor (or the Lode angle). The models are then formulated in the following sense :

ε f 0 f (stress state)dε = C (1.20)
where ε f is the equivalent fracture strain, and C the critical value. Some of these works rely on existing plasticity theories, or intend to establish new theories with different yield surfaces.

Within these criteria, one of the first so-called empirical model is the Cockroft-Latham-Oh criterion [START_REF] Cockcroft | Ductility and the workability of metals[END_REF] in which fracture is driven by the maximum principal stress over the equivalent stress integrated over the strain path. This model was modified by Brozzo et al. [START_REF] Brozzo | A new method for the prediction of formability in metal sheet, sheet Metal forming and formability[END_REF] to include a hydrostatic stress dependence, and Clift et al. [START_REF] Clift | Fracture prediction in plastic deformation processes[END_REF] replaced the maximum principal stress by the equivalent stress.

Another class of models relates fracture to the stress triaxiality state. For example, the Rice and Tracey criterion [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF] states that fracture exponentially depends on the stress triaxiality ratio.

Bao and Wierzbicki [START_REF] Bao | A comparative study on various ductile crack formation criteria[END_REF] have tested a few of these criteria and they concluded that none of these standard criteria could describe the aforementioned two different fracture mechanisms ('internal necking' and 'void sheeting'). They were accurate only in a specific range of triaxiality ratios, and specifically not in the middle range where the 'void sheeting' mechanism prevails.

These authors also emphasized [START_REF] Bao | On fracture locus in the equivalent strain and stress triaxiality space[END_REF] that the aforementioned authors did not carry sufficient tests in the intermediate stress triaxiality range, and that they could not rigorously assume that within this middle range the fracture locus follows the same rule as in the extremas.

They addressed the issue that the elaboration of tests within this range was particularly difficult, this is why the description of the fracture locus within this range is very challenging. They eventually obtained a complete fracture locus for the Al 2024-T351 aluminium alloy used in the aeronautical industry, which highlights their previsions on the middle range triaxiality ratio fracture locus. They also showed the existence of a cut-off value of the stress triaxiality ratio below which fracture never occurs, which was not emphasized by Rice and Tracey [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF] or Johnson and Cook [START_REF] Johnson | Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures[END_REF].

Recently, with a similar approach Bai and Wierzbicki [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF] have extended the Mohr-Coulomb yield criterion, that has been extensively used to describe fracture in the geomaterial community, to the range of ductile fracture. They reformulated this criterion, originally formulated in the principal stress space, in the space of equivalent fracture strain, stress triaxiality and Lode angle (ε f , η, θ). and included a pressure and a Lode angle dependence that were missing in the original formulation.

Figure 1.3 presents the fracture locus for the Johnson-Cook model and for the Bai-Wierzbicki model in a plane stress state, it illustrates the aforementioned differences between the two models. In the general stress state case, the fracture locus for Bai-Wierzbicki's model becomes a 3D surface.

Xue [START_REF] Xue | Damage accumulation and fracture initiation in uncracked ductile solids subject to triaxial loading[END_REF] along with other authors [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF] also emphasized with these studies, that This shows the importance of the Lode angle parameter in the ductile fracture process and the need to account for it in models.

As a consequence, authors have recently intended to describe the fracture strain locus with respect to both the triaxiality ratio and the Lode angle parameter, by several means. Some [START_REF] Bai | A new model of metal plasticity and fracture with pressure and Lode dependence[END_REF] completely reestablished standard J2 plasticity models by incorporating the third invariant and pressure dependencies on the yield condition, or they included this dependency in the damage accumulation process [START_REF] Xue | Damage accumulation and fracture initiation in uncracked ductile solids subject to triaxial loading[END_REF], both intend to account for the Lode angle dependence at different stages of the modeling.

We recall that at high triaxiality ratios, the prevalent mechanism is primary void growth which implies that coalescence happens sooner at high stress triaxiality than at low stress triaxiality level. However a shift occurs at stress triaxialities below 0.3 -0.4 where the void sheeting mechanism prevails. This induces a different shape of the fracture locus below this range, as it is shown in figure 1.3.

Lemaitre's model

Lemaitre [START_REF] Lemaitre | A continuous damage mechanics model for ductile fracture[END_REF] developed a theory of 'Continuum Damage Mechanics' based on the early works of Kachanov [START_REF] Kachanov | Rupture time under creep conditions[END_REF]. Some authors [START_REF] Besson | Continuum models of ductile fracture: a review[END_REF] criticize this terminology knowing that micromechanical models also use the notion of continuum damage. This is why they prefer to refer to this theory as "phenomenological" because they are based on macroscopic considerations.

Lemaitre's notion of damage relies on the macroscopic phenomenon of loss of carrying capacity when plastic strain is increasing during a loading. It is related to the weakening of Young's modulus which is equivalent to the loss of resisting section, as follows in the scalar version:

D = S 0 -S S 0 = E 0 -E E 0 (1.21)
where S 0 is the resisting section of the undamaged material, S is the current resisting section, E 0 is the Young modulus of the undamaged material and E the Young modulus of the current material. This definition of damage allows to introduce the important notion of effective stress, which corresponds to the stress across the undamaged section after having removed from the total section the part that has been damaged, which corresponds to voids and microcracks. It is defined as follows:

σ = σ 1 -D (1.22)
Lemaitre's damage variable can be scalar for isotropic materials in proportional loadings configurations, as in the previous definitions, or it can be tensorial for anisotropic materials and/or for non-proportional loadings.

Lemaitre's damage-plasticity model has been widely extended to different applications such as creep, fatigue or brittle fracture. These different applications are detailed in the book [START_REF] Lemaitre | Ductile, creep, fatigue and brittle failures[END_REF], and are written in a rigorous and consistent thermodynamical framework, in a way that dissipation always remains positive. Considering Lemaitre's isotropic damage plasticity model with isotropic and kinematic hardening, the model has four internal variables: the plastic strain tensor ε p , the scalar kinematic hardening variable α and the scalar damage variable D. The model is actually formulated like von Mises standard plasticity model, but replacing the true stress tensor σ by the effective stress tensor σ defined in (1.22). It relies on Lemaitre's principle of strain equivalence which states:

"The principle of strain equivalence states that any strain constitutive equation of a damaged material is written exactly as for a virgin material except that the stress is replaced by the effective stress."

The yield function associated to this model is then defined as:

Φ = 3 2 dev( σ -X) : dev( σ -X) -(σ y + q(r)) (1.23)
where dev is the deviatoric part of a second order tensor, X is the stress associated to kinematic hardening, q(r) is the hardening law (which can either be chosen as a power law or a saturated law), and σ y is the elastic limit.

Lemaitre's model is a non-associated model so the internal variables are defined by deriving the following dissipation potential:

G D = Φ + S (1 + s)(1 -D) Y S 1+s (1.24)
where Y is a thermodynamic force from which the damage variable derives, also known as the elastic energy release rate, s ans S are Lemaitre's material parameters.

The damage variable evolution law is defined by Lemaitre as:

Ḋ = λ 1 -D Y S s (1.25)
Since this model is used in the present work, its theoretical developments as well as its numerical implementation are detailed in section 2.2 in the following chapter.

Remark 1.2.1 (Lemaitre's model and anisotropic damage). As mentioned above, the damage variable can be tensorial in the case of anisotropic materials such as laminated plates or in the case of non-proportional loadings since in this case the normal to the resisting section is not constant with respect to the loading path. This comes from the fact that Lemaitre's damage variable, contrary to Gurson's porosity definition, is defined from areas and not volumes.

Thus, authors have tried to rebuild Lemaitre's model by considering a tensorial measure of damage. There exist different ways of defining a tensorial damage, it can either be a second order tensor (see e.g. [START_REF] Murakami | A continuum theory of creep and creep damage[END_REF][START_REF] Cordebois | Damage induced elastic anisotropy[END_REF][START_REF] Chow | An anisotropic theory of continuum damage mechanics for ductile fracture[END_REF]) or a fourth-order tensor (see e.g. [START_REF] Chaboche | Continuous damage mechanics -A tool to describe phenomena before crack initiation[END_REF][START_REF] Krajcinovic | Continuous damage mechanics revisited: basic concepts and definitions[END_REF]). The second order version is often prefered due to the complexity induced by the fourth-order version, and the definition of the effective stress tensor in the case of the second order damage tensor is actually an approximation of the fourth-order case (which corresponds to the most general case).

In the case of a second order damage tensor, the effective stress tensor may be written [START_REF] Lemaitre | Ductile, creep, fatigue and brittle failures[END_REF]:

σ = dev(H.dev(σ).H) + σ H 1 -ηD H 1 (1.26)
where dev(a) is the deviatoric part of tensor a, H is defined as H = (1-D) -1/2 , η is a material parameter that drives the effect of damage on the hydrostatic stress and on the deviatoric stress (called the hydrostatic sensitivity of damage), σ H and D H are the hydrostatic stress and hydrostatic damage, respectively. η is claimed to be easy to identify. As mentioned in [START_REF] Lemaitre | Ductile, creep, fatigue and brittle failures[END_REF], this is the only definition of the effective stress tensor using a second-order damage tensor that fulfills the conditions of symmetry of the effective stress, that is independent of the strain behavior, and that is compatible with the thermodynamics framework.

In this case of anisotropic damage, Lemaitre et al. [START_REF] Lemaitre | Anisotropic damage law of evolution[END_REF] have proposed to generalize the definition of the tensorial damage evolution law as:

Ḋ = Y S s | εp | (1.27)
where εp is the absolute value of the plastic strain tensor rate.

Remark 1.2.2 (Damage closure). Lemaitre's has been modified [START_REF] Lemaitre | A course on damage mechanics[END_REF][START_REF] Andrade Pires | Numerical modelling of ductile plastic damage in bulk metal forming[END_REF] to account for compressive behavior under cyclic loading or complex loading paths. This version takes into account the partial recovering of elastic properties when voids or microcracks close, which is not straightforward in the standard version. To this purpose, the strain energy release rate has been modified as:

Y = 1 + ν 2E < σ > + :< σ > + (1 -D) 2 + h < σ > -:< σ > - (1 -hD) 2 - ν E tr(< σ > + ) 1 -D 2 + h tr(< σ > -) (1 -hD) 2 (1.28)
where < a > + and < a > -represent the positive part and the negative part of tensor a, respectively. h is a parameter (such that h 1) that describes the effect of damage closure. If h = 1 there is no damage closure and the original model is recovered, and if h = 0 there is no damage growth since Y = 0 under compression.

The above definition of the energy release rate considering damage closure has been extended to anisotropy by Ju [START_REF] Ju | On energy-based coupled elastoplastic damage theories: Constitutive modeling and computational aspects[END_REF] and Chaboche [START_REF] Chaboche | Development and continuum damage mechanics for elastic solids sustaining anisotropic and unilateral condition[END_REF] for quasi-brittle materials and by Lemaitre [START_REF] Lemaitre | Anisotropic damage law of evolution[END_REF] and Desmorat et al. [START_REF] Desmorat | Modeling microdefects closure effect with isotropic/anisotropic damage[END_REF] for ductile materials.

Remark 1.2.3 (Other extensions of Lemaitre's model). Several extensions of the CDM (Continuum Damage Mechanics) theory have been proposed. Among them, we can cite the works of Saanouni [START_REF] Saanouni | Damage Modeling in Metal Forming Processes[END_REF][START_REF] Saanouni | Virtual metal forming including the ductile damage occurrence: Actual state of the art and main perspectives[END_REF][START_REF] Saanouni | On the numerical prediction of the ductile fracture in metal forming[END_REF] in which the authors enhanced the model to account for anisotropy with respect to several mechanical variables, such as the kinematic hardening, the elastic stiffness, the plastic yield function as well as the thermal effects, in the context of metal forming using a finite strain hypothesis. For example, the new yield function is written in Hill's sense [START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF] as follows:

Φ(σ, X, R, D, T ) = (σ -X) : H : (σ -X) -R √ 1 -D -σ y = 0 (1.29)
where X is the anisotropic kinematic hardening, T the temperature, and H is a fourth order tensor that accounts for an initial plastic anisotropy. The stress tensor is, in this context, defined as:

σ = (1 -D)Γ : ε e -(T -T 0 ) √ 1 -Dk (1.30)
where Γ is a symmetric fourth order tensor that takes into account the elastic stiffness anisotropy and k is a symmetric second order tensor that accounts for the anisotropic thermal elongation.

In the same context of metal forming simulations, Labergere et al. [START_REF] Labergere | Improving numerical simulation of metal forming processes using adaptive remeshing technique[END_REF] have worked on a remeshing technique based on various types of error estimators in order to circumvent the issue of having a "process zone" smaller than the characteristic length (that is the mesh size).

They proposed an error estimator based on the derivatives of the plastic strain (the Hessian matrix of the plastic strain). The computation of the Hessian of the plastic strain is rendered possible by using the diffuse approximation technique for which the local derivatives are immediately available since the approximation is polynomial.

They also proposed an error estimator based on the damage rate Ḋ. This remeshing method was implemented in Abaqus explicit for the solver together with the software DIAMESH2D for the remeshing process (see [START_REF] Labergere | 2D adaptive mesh methodology for the simulation of metal forming processes[END_REF]).

The Rousselier model

Rousselier introduced a model [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF] based on Lemaitre and Chaboche thermodynamical framework, that describes void growth for ductile materials which leads to a quite similar yield surface to Gurson's model. He developed his model with the aim of obtaining a simple model in a consistent thermodynamic framework which depends on the two first invariants of the stress tensor. Rousselier uses the porosity f as an internal variable in such a way that it satisfies the equation of mass conservation stated as:

ḟ = (1 -f )tr( εp ) (1.31)
The yield function is expressed as:

Φ = σ eq 1 -f + σ 1 f D exp tr(σ) (1 -f )σ 1 -R (1.32)
where σ 1 is a material parameter, D is a parameter that is supposed to be independent of the material and R is the isotropic hardening which depends on the accumulated plastic strain p. Although the Rousselier model looks like the Gurson model, they do not enable the same growth mechanisms for different stress states.

For instance, if one considers a pure shear stress state in which tr(σ) = 0 damage growth still occurs with the Rousselier model whereas it does not evolve for the Gurson model.

One of the main advantage of Rousselier's model is that, contrary to Lemaitre's model, it enables to model volume change and density change through the use of porosity as a damage measure, within a consistent thermodynamical framework such as the one proposed by Lemaitre [START_REF] Lemaitre | A continuous damage mechanics model for ductile fracture[END_REF]. It also provides as good predictions as Gurson's model while it is claimed to better capture strain and damage localization.

The presence of a vertex point on the yield surface also enables to preserve a shear plastic flow even at high stress triaxialities.

The Rousselier model has been extended by Tanguy et al. [START_REF] Tanguy | An extension of the Rousselier model to viscoplastic temperature dependent materials[END_REF] to viscoplastic and temperature dependent materials. Remark 1.2.4 (Links between micro-mechanical and phenomenological modeling). Both types of modeling have their advantages and drawbacks. As discussed above, Lemaitre's phenomenological model doesn't account for volume change and thus is not able to properly model void growth, it is essentially related to void nucleation. The main advantage of such a model is its thermodynamically consistent formalism contrary to Gurson's model. In order to be able to use this formalism and to recover the predictive mechanisms obtained with the GTN model, Lemaitre and Desmorat [START_REF] Lemaitre | Ductile, creep, fatigue and brittle failures[END_REF] showed that it is possible to tune the parameters of Lemaitre's model to simulate void growth. A further attempt has been made by Chaboche to combine Lemaitre's formalism for void nucleation together with the GTN void growth model in order to obtain an elliptic model.

Another advantage of Lemaitre's model is that it enables to simulate the effects of damage on elastic properties, and this coupling induces complex numerical integration procedures. Some authors have neglected this coupling [START_REF] Xue | Damage accumulation and fracture initiation in uncracked ductile solids subject to triaxial loading[END_REF] to simplify the model by defining an independent damage evolution law.

1.3 Mesh dependency issue: a review of the localization limiters

Highlighting the mesh dependency issue

Ductile fracture models often induce softening behaviors since, for instance in the case of Lemaitre's model, the development of damage induces a loss of load-carrying capacity. This type of behavior leads to strain localization and in the context of FE modeling, this leads to a pathological dependence of the solution on the mesh size and type. This comes from the fact that FE discretization uses a characteristic length which is the mesh size. When strain localization occurs, the dissipative mechanisms are localized in a small area that depends on the mesh size. Since the total dissipation depends on the area in which dissipation occurs, if standard FE are used, the mesh size inevitably appears in the total dissipation and in the extreme case of an infinitely small mesh size, fracture would occur without any dissipation.

An example of this phenomenon is presented on figure 1.4, in which a simple tensile test is simulated using Lemaitre's isotropic damage model in a plane strain framework. Three meshes using Q4 elements with full integration and with different mesh sizes are considered, an element at the center of the specimen is "weakened" to simulate a micro-defect in order to drive localization, the yield stress σ 2 y of this element is equal to 0.99σ It is obvious that the post-peak phase of the force/displacement curve strongly depends on the mesh size. The dissipated energy during this phase tends to converge to 0 as stated above. It has to be noticed that the global post-peak phase which corresponds to the force/displacement curve does not necessarily coincide with the local post-peak phase at the material point. Indeed, local softening generally occurs before global softening. At the end of the test, the major part of the strain is localized in the weakened element. This is why the total dissipation of the softening phase reduces to 0 as the element width reduces to 0.

This phenomenon is also related to the fact that standard elements are designed to model volumic dissipative mechanisms and not to model surfacic dissipative mechanisms which are inseparable from ductile fracture modeling. They are mainly divided in two types of techniques, local models and non-local models. In this section, our interest is to review a few of these methods, before introducing our method based on the Strong Discontinuity Approach (SDA).

Local methods

Viscous regularization

One of the techniques consists in introducing an artificial viscosity in the model. This local method is based on the fact that the rate dependency of a material induces a regularizing effect. Indeed, in a mathematical point of view, governing equations involving material rate dependency remain elliptic, even with softening materials. This was proved by the work of Needleman [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF] in which he showed with a one-dimensional model under quasi-static and dynamic loadings that a viscoplastic effect removes the pathological mesh-dependency issue by introducing an artificial length scale in the governing equations, which keeps them elliptic. This method was used by Dube et al. [START_REF] Dubé | Rate dependent damage model for concrete in dynamics[END_REF] in the case of a one-dimensional model for concrete under dynamic loading, he numerically showed that the objectivity of the solution was ensured. This type of regularization was also implemented by Sluys et al. [START_REF] Sluys | Wave propagation and localisation in a rate-dependent cracked medium: Model formulation and one-dimensional examples[END_REF] in the case of dynamic loadings in which they determined an explicit expression of the viscoplastic length scale for regularization.

Recently, Niazi et al. [START_REF] Niazi | Viscoplastic regularization of local damage models: revisited[END_REF] have demonstrated the regularizing capabilities of a viscoplastic regularization applied to a modified version of Lemaitre's model for quasi-static problems, considering anisotropic damage. They showed that a so-called viscoplastic effective regularizing length decreases as damage grows, and if this length becomes shorter than the mesh size, the problem becomes mesh-dependent again. These authors recommend to use very fine meshes to avoid this issue.

The limitation of this method lies in the fact that the introduced viscous parameters' physical significations are not necessarily established in the case of a non-viscous plasticity/damage model, while they are only used for their numerical regularizing effects. Besides, the use of these parameters involves characteristic times considerations along with limited time steps.

Parameters depending to mesh size

Another kind of method, probably the most easy to implement, is those that involve parameters that depend on the mesh size. The philosophy of this approach consists in adjusting, with regards to the mesh size, softening parameters in order to ensure the objectivity of the solution. This method was first introduced by Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt [START_REF] Barenblatt | The mathematical theory of equilibrium crack in the brittle failure[END_REF].

A method based on this approach and on non-linear fracture mechanics was then introduced by Hillerborg et al. [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF] in the framework of quasi-brittle materials. It consists in enforcing the computed fracture energy G f to remain independent from the element size. The so-called "fracture process zone" is determined through a FE computation and material parameters governing the softening phase are tuned to give an equal total fracture energy. This energy depends on the number of elements crossed by a displacement discontinuity, as a consequence these parameters depend on the sizes of the elements crossed by discontinuities.

An obvious drawback of this method is the non-physical sense of these elementsize dependent parameters. Another one is the need of a new set of parameters for each FE computation involving a new mesh. On the other hand, it is quite simple to implement in a FE code. Only the material law requires an adjustement with this method.

Cosserat models

The Cosserat brothers developed a generalized continuum theory [START_REF] Cosserat | Théorie des corps déformables[END_REF], also known as the micro-polar continuum theory, that has then been used to treat localization problems at the beginning of the 90's by De Borst and Sluys [START_REF] Borst | Simulation Of strain localization: a reappraisal of the Cosserat continuum[END_REF] and Yu et al. [START_REF] Yu | The finite element method in the 1990's, An Assessment of the Cosserat Continuum through the Finite Element Simulation of a Strain Localisation Problem[END_REF]. The key idea is to introduce additional rotational degrees of freedom within the equilibrium equations, which are treated within an additional momentum equilibrium equation. This enriched continuum formulation implies that a bending modulus appears in the equations that automatically incorporates a characteristic length essential for regularization. These models have also been used in a certain extent to model granular materials by Mulhaus and Vardoulakis [START_REF] Mühlhaus | The thickness of shear bands in granular materials[END_REF] in which Cosserat's theory is used to evaluate the thicknesses and orientations of shear bands for materials with internal friction. More recently, Khoei and Karimi [START_REF] Khoei | An enriched-FEM model for simulation of localization phenomenon in Cosserat continuum theory[END_REF] reemployed Cosserat's theory to model shear bands in ductile materials together with a X-FEM formulation. Cramer et al. [START_REF] Cramer | An approach to the adaptive finite element analysis in associated and non-associated plasticity considering localization phenomena[END_REF] proposed an adaptative mesh method in the context of associated and non-associated plasticity using Cosserat's continuum theory.

The generalized continuum theory developed by the Cosserat brothers have brought a certain attention in the strain localization modeling for ductile fracture since the beginning of the 90's, and are essentially used to model shear bands.

One of the drawbacks is that it is necessary to reformulate the material constitutive equations since couple stresses terms appear along with the traditional Cauchy stress tensor.

Non-local methods

Gradient models

Gradient models can be of different types, those who involve a gradient formulation directly within the equilibrium equation, and those using a gradient of internal variables instead of the variable itself in the material model formulation.

The first kind of gradient models uses a high order velocity gradient, in which the equilibrium equations are modified. The founding principles of these methods were introduced by Toupin (see [START_REF] Toupin | Elastic materials with couple-stresses[END_REF]) and Mindlin (see [START_REF] Mindlin | Microstructure in linear elasticity[END_REF]). The key point is to postulate that the principle of virtual power includes the use of the second gradient of velocity ∇∇v, in addition to the deformation rate ε = ∇ s v (which corresponds to the first gradient of velocity).

The standard formulation of the principle of virtual power is actually a restriction of the second gradient theory, where only the first gradient of velocity appears in the formulation. Although this method relies on a very rigourous thermodynamical framework and is quite elegant, its numerical implementation is complicated by the requirement of the use of shape functions that are able to describe a second gradient of velocity, on the one hand. On the other hand, writing the boundary conditions is quite difficult because they need to be written for the dual variable of the second gradient of velocity which is a third-order stress tensor, and this is not an easy task.

The other type of gradient models use a gradient of internal variables. The introduction of the gradient of a certain internal variable within the constitutive equations or in the evolution equations allows the regularization. While integral models use a convolution of a variable over a given neighborhood instead of the variable itself, these models use their spatial derivatives.

Among these models, we can cite Mulhaus et al. [START_REF] Mulhaus | A variational principle for gradient plasticity[END_REF] which was developed for a softening plasticity model. The Laplacian of the hardening variable is used to compute its dual variable that appears in the plasticity yield function. A generalization of this assumption was proposed by Lorentz et al. [START_REF] Lorentz | A variational formulation for nonlocal damage models[END_REF] in which the free energy is rewritten to account for the spatial derivative dependency of the dual variables to the state variables. Then the differentiation of the free energy to obtain the state equations naturally links the dual variables to the state variables and their gradients. In their works [START_REF] Lorentz | A variational formulation for nonlocal damage models[END_REF] Lorentz et al. showed that the gradient plasticity model proposed by Mulhaus et al. [START_REF] Mulhaus | A variational principle for gradient plasticity[END_REF] is a particular case of their general framework. They also proposed in their paper a damage homogenization and a porosity model based on their general gradient framework.

In a similar way Reusch et al. [START_REF] Reusch | A non-local extension of Gurson-based ductile damage modeling[END_REF] proposed a modification of the GTN model to account for non-local damage. They replaced the porosity f that appears in the yield function by a so-called time-dependent continuum damage field v, assuming that the effect of the damage field v is the same as the one of f :

Φ = σ eq σ y 2 + 2q 1 vcosh 3 2 q 2 σ h σ y -1 -(q 1 v) 2 (1.34)
An additional term related to the gradient of this damage field is also added to the free energy:

Ψ = Ψ s (ε, ε p ) + 1 2 α D l 2 D ∇v.∇v (1.35)
where Ψ s is the standard free elastic energy for the GTN model, l D is the characteristic length of the non-local damage field and α D is a characteristic time.

Recently, Linse et al. [START_REF] Linse | Simulation of crack propagation using a gradient-enriched ductile damage model based on dilatational strain[END_REF] proposed a different modification of the GTN model to account for non-local damage. The non-local aspect appears in the formulation of the porosity rate, in which the dilatiational part of the plastic strain rate is replaced by its average over a neighborhood, whereas in the previous model a new damage field variable was introduced:

ḟ nl g = (1 -f ) εp (1.36)
where εp is defined implicitly as the solution of the diffusion differential equation expressed as:

εp -c∇ 2 εp : I = ε p (1.37)
where ∇ 2 () : I reduces to the Laplace operator, √ c has a dimension of length and plays the role of the characteristic length responsible for regularization. It is often more numerically convenient to use the gradient of plastic strain instead of the gradient of damage since damage is limited to the range [0, 1] contrary to the plastic strain which does not have a superior limit.

Other extensions of non-local damage modeling along with the Gurson model have been proposed by [START_REF] Tvergaard | Effects of nonlocal damage in porous plastic solids[END_REF][START_REF] Leblond | Bifurcation effects in ductile metals with nonlocal damage[END_REF][START_REF] Jackiewicz | Non-local regularization for FE simulation of damage in ductile materials[END_REF], and with the Rousselier model by [START_REF] Samal | Finite element formulation of a new nonlocal damage model[END_REF].

A similar implicit gradient of damage formulation has been formulated as an extension of the Lemaitre damage model by César de Sá et al. [START_REF] César De Sá | Damage modelling in metal forming problems using an implicit non-local gradient model[END_REF]. A non-constitutive damage variable D is implicitly related to the local constitutive damage variable D through the diffusion differential equation (as in (1.37)):

(D -D) + l 2 0 ∇ 2 0 D = 0 (1.38)
where l 0 is the characteristic length scale that drives the area of damage diffusion. Then this non-local damage variable D replaces the standard local damage D variable in the Lemaitre constitutive equation. In their works César de Sá et al. implemented this extended model in a mixed-enhanced strain FE framework, which is convenient when non-local field variables are considered.

One of the drawbacks of these models is the difficulty to implement them in a finite element program since the gradients of internal variables are usually treated as additional degrees of freedom, which implies the use of non-standard FE procedures since the added "dofs" cannot be condensed at the element level due to the non-local property of these variables. The use of gradient of internal variables also implies that high-order finite element are used in order to properly represent these gradients. Another drawback is the difficulty induced with regards to boundary conditions since the gradients should naturally vanish at boundaries which has to be enforced in the model. Usually a simple Neumann boundary condition is introduced as an additional condition:

∇ D.n = 0 on ∂Ω (1.39)
and is weakly implemented in the FE procedure.

Integral models

We now focus on models that involve the use of non-local variables through the use of an integral over a given neighborhood. The use of non-local integral models involves the replacement of some scalar variables by their mean over a neighborhood V , using a weighting function that drives the size of the neighborhood. They were firstly introduced by Pijaudier-Cabot et al. (see [START_REF] Pijaudier-Cabot | Non-local damage theory[END_REF]).

For example, if one considers a local variable f (x), the corresponding non-local averaged field is defined as:

1.3. Mesh dependency issue: a review of the localization limiters

f (x) = V α(x, s)f (s)ds (1.40)
where s defines the "source" term and α(x, •) is the weighting function around point x. In order to reproduce exactly uniform fonctions, the following condition has to be satisfied:

V α(x, s)ds = 1 ∀x ∈ V (1.41)
Finally, the standard average function is defined as:

α(x, s) = α ∞ (||x -s||) V α ∞ (||x -s||)ds (1.42)
where α ∞ (r) is usually chosen as a bell-shaped function that tends to zero when r tends to the boundary of the volume V .

The first immediate limitation of this method is the choice of the weighting function that appears in the integral over the neighborhood, Planas et. al showed that the model is very sensitive to this function (see [START_REF] Planas | Cohesive cracks versus non-local models : closing the gap[END_REF]). Another one is the choice of the characteristic length associated to the weighting function, which also has a very important influence on the response of the model. Jirasek [START_REF] Jirásek | Rotating crack model with transition to scalar damage[END_REF] performed an exhaustive overview of non-local integral damage models based on Lemaitre's damage theory [START_REF] Lemaitre | A continuous damage mechanics model for ductile fracture[END_REF], and compared results with regards to the type of averaging that is undertaken and the formulation of damage model itself. Indeed, it is possible to perform the averaging on different variables: the damage variable D itself, the energy density release rate Y , the equivalent strain (considering the non-local extension [START_REF] Saouridis | Prediction of the failure and size effect in concrete via bi-scale damage approach[END_REF] of the damage model for concrete developed by Mazars [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure[END_REF]), the inelastic stress s = D e ε e -σ...etc... He eventually showed with a simple tensile test that averaging quantities like the damage variable, or the inelastic stress induced a non-negligeable locking effect due to the fact that, when averaging such quantities, the stress was not reduced to zero even for very large imposed displacements, which is not realistic. On the contrary averaging quantities like the energy density release rate or the inelastic stress where consistent and did not induced such unrealistic behavior.

Non-local plasticity models were first introduced by Eringen [START_REF] Eringen | On nonlocal plasticity[END_REF] with no aim of proposing a localization limiter. It then inspired authors like Bazant and Lin who proposed a non-local plasticity model with strain softening [START_REF] Bazant | Non-local yield limit degradation[END_REF] in which the averaged variable is the plastic multiplier computed as:

λ = V α(x, s) λ(s)dV (1.43)
and replaces the local plastic multiplier, which is computationally more efficient than averaging the plastic strain tensor ε p itself since it is scalar. However, this formulation still gives rise to the aforementioned locking effects since it is equivalent to averaging the inelastic stress.

Many other non-local plasticity models have been developed using the same averaging philosophy, some of these models were reported by Bazant and Jirásek in [START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: survey of progress[END_REF]. The main diffences between these models are in the choice of the averaged variables.

Finally, similar to gradient models, as these models are difficult to implement and require information on the neighborhood's variables, important modifications of the global architecture of the finite element resolution need to be carried out.

Discontinuity models

A last type of regularizing method reviewed in this document is the one involving discontinuities of the displacement or the strain field. The aim of these models is not to precisely describe at the micro-scale the localization phenomena, but to account for them at the macro-scale.

A major difference between these methods and the aforementioned ones is that the finite element enhancement only involves modifications at the element level, contrary to non-local methods for example that require the knowledge of some quantities at the vicinity of the considered point, that is the surrounding elements.

The main philosophy is to introduce an incompatible strain field within an element, that can describe a discontinuous strain field (in the case of weak discontinuities) or a discontinuous displacement field (in the case of strong discontinuities). They require the use of incompatible shape functions. The weak discontinuity models are briefly described in the following section while the strong discontinuity formulation (which is the method used in this work) is detailed in a subsequent section.

Weak discontinuities

The weak discontinuities models were initially introduced by Pietrusczak et al. [START_REF] Pietrusczak | Finite element analysis of deformation of strainsoftening materials[END_REF], and then by Ortiz et al. [START_REF] Ortiz | A finite element method for localized failure analysis[END_REF]. The idea is to build elements capable of representing discontinuous strain fields within an element. When a localization mode is detected, the shape function basis is enhanced with incompatible shape functions that can take into account a discontinuous strain field with an orientation (n, m) given by the localization analysis (n is the normal to the discontinuity surface while m is the tangential vector to the discontinuity). For every localized mode (I or II) the strain jump is represented by a single parameter which is considered as an additional degree of freedom.

One drawback of this method is that the additional shape functions do not verify the C 0 continuity at the element boundaries. However, since the additional degrees of freedom are specific to the localized element, using the three field Hu-Washizu variational principle it is possible to eliminate these dofs at the global level by static condensation at the element level. Thus, this method can be implemented in a standard FE industrial code. Another limitation of the model proposed by Ortiz [START_REF] Ortiz | A finite element method for localized failure analysis[END_REF] is that it does not properly introduce a characteristic length necessary to regularize since a localization band is characterized by two strain discontinuity surfaces. Thus, it cannot be considered as a localization limiter.

Belytschko et al. [START_REF] Belytschko | A finite element with embedded localization zones[END_REF] addressed this problem and enhanced this principle by allowing the introduction of a full localization band in a three node element, which means that two surfaces of strain discontinuity are introduced within the element. In this model, the plastic deformation is entirely localized in this band, and the out-of-band material is supposed to remain elastic. Contrary to the work of Ortiz [START_REF] Ortiz | A finite element method for localized failure analysis[END_REF], no additional shape function is introduced in the element, however the effect of the localization band behavior modifies the element behavior which is taken into account by the formulation of the three field Hu-Washizu variational principle.

An important parameter of such a model is the relative area of the localization band and the non-localized zone which implies that a minimum element size is required, in order to be fully crossed by a localization band with a given width. This represents a limitation of this method, indeed a minimum element size is required. Another drawback is that the position of the localization band cannot be properly described.

Recently and based on the work of Ortiz and Belytschko, Huespe et al. [START_REF] Huespe | A finite thickness band method for ductile fracture analysis[END_REF] proposed to introduce a full localization band in a four node element whose kinematics are described by an additional shape function capable of representing a discontinuous strain field across the band. The new shape function is plotted on figure 1.5(a) in the case of a one dimensional element. Considering a one-dimensional finite element discretization, the total displacement is expressed as:

u(x, t) = ū(x, t) + M Ω d (x)α(t) (1.44)
and in an element e:

u e (x, t) = n nodes i N e i (x)d i (t) ū(x,t) + (H Ω d (x) - i∈Ω + N e i (x)) M e Ω d (x) α e (t) (1.45)
where α e is the displacement jump vector across the band. H Ω d can be seen as a penalized Heaviside function. Indeed, instead of jumping from 0 to 1 discontin-uously, it linearly increases in the band (see figure 1.5(a)), which keeps the total displacement field continuous in the element.

The equilibrium equations are written following a Petrov-Galerkin formulation and the additional equilibrium equation weakly imposes the traction continuity across the localization band.

Contrary to the Ortiz's work [START_REF] Ortiz | A finite element method for localized failure analysis[END_REF], the new shape function satisfies the C 0 continuity condition since it vanishes at the element boundary (see figure 1.5(a)).

The model falls into the framework of ductile fracture modeling since the outof-band material is described by the GTN model described in section 1.2.1. Huespe et al. performed a rigorous bifurcation analysis (based on the singularity of the acoustic tensor) to obtain a criterion that triggers the localization band introduction and that defines its orientation. They performed calculations in which mode I and II are activated in axisymmetric notched-bar tensile test. However, they did not perform mesh dependency tests with their model, the aim of this work was to finely describe the localization band behavior.

Mesh dependency issue: a review of the localization limiters

The philosophy of weak discontinuities introduced by Pietrusczak and Ortiz represents the starting point of the more recently developed strong discontinuity methods, which is discussed in the following paragraph.

Strong discontinuities

The foundings of strong discontinuity analysis were proposed by Dvorkin [START_REF] Dvorkin | Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions[END_REF] and Klisinsky and Simo et al. [START_REF] Simo | An analysis of strong discontinuity induced by strain softening solutions in rate-independent solids[END_REF], and developments of this method have been proposed since then [START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures. Part I : theoretical formulation and numerical implementation[END_REF]. The key idea is to introduce a surface of discontinuity of displacement that has a zero measure, and to concentrate the whole dissipation on this surface. Hence, the regularization is ensured since a characteristic length is automatically introduced (and is equal to zero). An immediate consequence is that there is no neighboring issue as in the non-local models.

The main purpose being the regularization of the softening behavior, there is no purpose in finely describing the micro-scale mechanisms (by fully describing the localization band for example). The idea is to enhance the local model with surfacic dissipative mechanisms at the element scale in a FE framework, and to account for their effects at the macro-scale in order to regularize the solution. A straightforward consequence is the necessity to enhance the finite element formulation with shape functions capable of representing discontinuous displacements across the newly introduced surface.

Considering a one-dimensional finite element discretization, the total displacement is written (to be compared with equation (1.44)):

u(x, t) = ū(x, t) + M Γ (x)α(t) (1.46)
where in an element e:

u e (x, t) = n nodes i N e i (x)d i (t) ū(x,t) + (H Γ (x) - i∈Ω + N e i (x)) M e Γ (x) α e (t) (1.47)
In this case, the Heaviside function H Γ is not regularized and is discontinuous (see figure 1.5(b)). The formulation is actually the same as proposed by Huespe in [START_REF] Huespe | A finite thickness band method for ductile fracture analysis[END_REF] but with Ω d being of zero width. Derivation of the Heaviside function implies that the total strain field becomes singular on the discontinuity surface within the element. This total strain field is then expressed as:

ε e (x, t) = n nodes i B e i (x)d i (t) ε(x,t) - i∈Ω + B e i (x)α e (t) + (α e (t) ⊗ n)δ Γ (x) ε (x,t) (1.48)
where B e i is the standard elementary gradient of shape functions matrix for node i, n is the unit normal vector to the discontinuity surface, and δ Γ is the dirac distribution associated to the discontinuity surface.

Different philosophies have been proposed to treat the strain singularity. Some authors proposed that a strain singularity is compatible with standard continuum mechanics, and so the displacement jump and the singular strain directly appear in the stress evolution equations, which finally gives a generalized constitutive equation that includes both the bulk material law and the traction/displacement jump law. The compatibility equations also allow to write a straightforward localization condition (the condition that allows the introduction of a discontinuity).

Another philosophy, proposed by Wells and Sluys [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF][START_REF] Wells | Three-dimensional embedded discontinuity model for brittle fracture[END_REF] and Ibrahimbegovic and Brancherie [START_REF] Ibrahimbegovic | Combined hardening and softening constitutive model of plasticity: precursor to shear slip line failure[END_REF], states that a distinct constitutive law related to the displacement jump allows to deal with the singularity, and allows to model the desired dissipative mechanisms. It is typically chosen as a cohesive law (softening traction/displacement jump law).

Jirásek proposed a review [START_REF] Jirásek | Comparative study on finite elements with embedded discontinuities[END_REF] on the different methods used by authors to treat the singularity of the strain field. In all cases, the finite element formulation relies on the three-field Hu-Washizu variational principle. The differences between those formulation mainly lies on the element formulation, and more precisely on the hypothesis concerning the matrix of incompatible shape functions, for example whether it can represent the kinematics of the displacement discontinuity properly (formulation named KOS: Kinematically Optimal Symmetric) or whether it can satisfy the traction continuity condition across the discontinuity (formulation called SOS: Statically Optimal Symmetric).

In the present work, the formulation used is the one proposed by Simo and Oliver in [START_REF] Simo | A new approach to the analysis and simulation of strain softening in solids, Fracture and Damage in Quasibrittle Structures[END_REF] and is called SKON (Statically and Kinematically Optimal Non-symmetric). This formulation states that the enhanced elements are capable of properly representing the kinematics of a displacement discontinuity, as well as satisfying the traction continuity condition on the discontinuity surface, which leads to a nonsymmetric formulation of the problem. The formulation is detailed further in this document.

Originalities of the present work

It has to be noticed that most of the works relying on the strong discontinuity methods imply that all the dissipative mechanisms are gathered on the discontinuity surface cohesive law, while the behavior in the bulk remains elastic. This means that this method has been mainly used to model brittle fracture. The works on the Strong Discontinuity Approach adapted to inelastic behavior mainly apply to the standard 1.3. Mesh dependency issue: a review of the localization limiters von Mises plasticity material law for the bulk. Hence, the main objective and one originality of this work is to account for plasticity as well as damage, using Lemaitre's fully coupled model for the pre-peak phase, along with the Strong Discontinuity Approach to account for surfacic dissipative mechanisms occuring in the context of ductile fracture.

This means that two different dissipative mechanisms interact with each other, a volumic dissipative mechanism that accounts for plasticity and damage growth, and a surfacic dissipative mechanism that allows to take into account strain localization as well as the development of a macro-crack and its propagation, the philosophy being that we do not intend to precisely describe the apparition of localization bands and their behavior, but to account for their effects at the macro-scale.

As the Strong Discontinuity Method can be seen as a localization limiter technique, we propose in chapter 2 a regularized version of Lemaitre's elasto-plastic damageable model. In chapter 3, we propose an enhancement of the Strong Discontinuity Method, initially restricted to a small strain hypothesis, to the context of finite strain modeling, this hypothesis being more suitable to model ductile fracture that involve large strains and displacements.

Chapter 2

Ductile fracture modeling with the strong discontinuity method in small strain hypothesis

In this chapter, the SDA theoretical formulation and its numerical implementation are described in the context of Lemaitre's phenomenological coupled damage-plasticity model. A four-node quadrangular element enhanced with a linear displacement jump is described, and numerical examples showing the capabilities of the method are presented. A critical damage value is used as a criterion to introduce the discontinuity surface. All numerical developments have been implemented in FEAP, an open source Finite Element software developed by R.L. Taylor at the beginning of the 90's (see [START_REF] Taylor | FEAP website[END_REF] for FEAP's website and [START_REF] Zienkiewicz | The Finite Element Method: Its Basis and Fundamentals[END_REF][START_REF] Zienkiewicz | The Finite Element Method For Solid and Structural Mechanics[END_REF] for the books that the program relies on).

Theoretical formulation

The Strong Discontinuity Approach (also called "Embedded Discontinuity Finite Element Method") was initiated and widely developed since the beginning of the 90's (see [START_REF] Armero | An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids[END_REF], [START_REF] Larsson | Embedded localization band in undrained soil based on regularized strong discontinuity-theory and FE-analysis[END_REF], [START_REF] Larsson | Finite element embedded localization band for finite strain plasticity based on a regularized strong discontinuity[END_REF], [START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 1: Fundamentals[END_REF], [START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 2: Numerical simulation[END_REF], [START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures. Part I : theoretical formulation and numerical implementation[END_REF]) mostly in the framework of brittle fracture, and results attesting of its performance were shown, in terms of regularization ability.

Recent works have implemented the method in the context of oxide layer fracture [START_REF] Li | Prediction of plane strain fracture of AHSS sheets with post-initiation softening[END_REF] or in the context of reinforced-concrete frames bending failure [START_REF] Pham | Stress-resultant models for ultimate load design of reinforced concrete frames and multi-scale parameter estimates[END_REF]. In the following section, a detailed presentation of the theoretical formulation is carried out, while we limit the study to the small deformation case in this chapter.

Let us consider a surface Γ s of unit normal vector n and tangent vector m crossing entirely a domain Ω and dividing this domain into two subdomains Ω + and Ω -(see figure 2.1) such that:

Ω + = { x | x.n > 0 } Ω -= { x | x.n < 0 } (2.1)
Considering the introduction of this surface of discontinuity of displacement, a discontinuous part of the displacement ū(x, t)H Γs (x) is added to the standard displacement ū(x, t) field in (2.2) (see figure 2.2(a) for a schematic representation in a 1D case).

u(x, t) = ū(x, t) + ū(x, t)H Γs (x) (2.2)
where H Γs (x) is the Heaviside function:

H Γs (x) = 1 if x ∈ Ω + 0 if x ∈ Ω - (2.3)
With such a decomposition, the essential boundary conditions involve both the continuous and the discontinuous part of the displacement, (see figure 2.2(a)). In order to limit the contribution of the discontinuous displacement field to the neighborhood of the discontinuity Γ s , a new decomposition of the total displacement field is proposed (see figure 2.2(b) for the 1D case).

u(x, t) = ũ(x, t) + ū(x, t)[H Γs (x) -φ(x)] (2.4) with ũ(x, t) = ū(x, t) + ū(x, t)φ(x) (2.5)
and φ a smooth, at least C 1 function verifying:

φ(x) = 0 if x ∈ Ω -\ Ω - d 1 if x ∈ Ω + \ Ω + d (2.6)
where Ω d is a subdomain of Ω containing the discontinuity Γ s .

For such a decomposition, the contribution of the displacement jump ū(x, t) vanishes on the boundary of the domain so that the essential boundary conditions written on the total displacement u(x, t) can be prescribed to the regular field ũ(x, t). This point is of great importance for the local aspect of the method as it will be detailed in the following paragraph.

Results on derivation of distributions (see [START_REF] Stakgold | Green's functions and boundary value problems[END_REF]) give the strain field as:

ε(x, t) = ∇ s u(x, t) = ∇ s ū(x, t) + ∇ s [ ū(x, t)H Γs (x)] (2.7)
And finally,

ε(x, t) = ∇ s ū(x, t) + ∇ s ū(x, t)H Γs (x) ε(x,t) + ū(x, t) ⊗ n s δ Γs (2.8)
where (•) s denotes the symmetric part of (•). It has to be noticed that the total strain field is now decomposed into a regular part ε and a singular part ū(x, t) ⊗ n s δ Γs , where δ Γs is the Dirac distribution associated with the discontinuity surface.

A four-node quadrangular element enhancement

Introducing a displacement field discontinuity requires an enhancement of the standard shape functions. We choose to rely on the incompatible modes method which enables the modifications to be limited to the element level. Other methods such as X-FEM ( [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF], [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], ...) use a different type of enhancement, which involves a specific treatment at global level. We consider a quadrangular element with four nodes and a linear displacement jump.

On the basis of the works of [START_REF] Linder | Finite elements with embedded strong discontinuities for the modeling of failure in solids[END_REF] and [START_REF] Dujc | Quadrilateral finite element with embedded strong discontinuity for failure analysis of solids[END_REF], the displacement jump is projected on a basis (n 0 , n 1 , m 0 , m 1 ) that corresponds to the different modes of separation of a crack, in the linear case, namely: a constant normal separation mode n 0 , a linear normal mode n 1 equivalent to the rotation of the discontinuity around the discontinuity midpoint x Γs , a constant tangential separation mode m 0 and a linear tangential separation mode m 1 . The amplitude of each separation mode is gathered into a vector α =<

α n 0 α n 1 α m 0 α m 1 > T .
The kinematics proposed (described in figures 2.3, 2.4, 2.5 et 2.6) allow to represent every possible separation modes in two dimensions, in the case of a linear displacement jump.

Remark 2.1.1. All the considered separation modes produce zero strain in the bulk except the one along "m 1 " which produces a deformation in Ω + .

αn0 αn0 m n Γ s Ω - Ω + 1 2 3 4
x Γs The total displacement field is then interpolated by adding shape functions associated to each of these separation modes, in the following form:

u h (x) = a=4 a=1 N a (x)u a +α n0 H Γs (x).n - a∈Ω + N a (x).n +α m0 H Γs (x).m - a∈Ω + N a (x).m +α n1 - a∈Ω + N a (x)(n • m T -m • n T ) • (x a -x Γs )+ (n • m T -m • n T ) • (x -x Γs )H Γs +α m1 - a∈Ω + N a (x) (x a -x Γs ) • m m+ (x -x Γs ) • m mH Γs (2.9)
where N a is the standard shape function associated to node "a", x Γs and x a are the coordinate vectors of the midpoint of the element and node "a" respectively. Remark 2.1.2. In the case of mode n 0 alone, the function φ introduced in section 2.1.1 takes the following form:

φ(x) = a∈Ω + N a (x)n (2.10)
This expression in (2.9) can be put in the general following form:

u h (x) = a=4 a=1 N a (x)u a + Mα (2.11)
where M is a discretized version of the continuous function H Γs (x) -φ(x) introduced in equation (2.4). This function allows the contribution of the separation modes to be limited to a localized element with Ω d (see section 2.1.1) chosen equal to the considered quadrangular element. This guarantees the local feature of the proposed method and no additional nodal degrees of freedom have to be introduced like for X-FEM for instance [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF] [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF].

On the basis of this displacement field, the real strain interpolation is given by:

ε h (x) = a=4 a=1 B a (x)u a +α n0 - a∈Ω + B a (x).n + N.n δ Γs +α m0 - a∈Ω + B a (x).m + N.m δ Γs +α n1 - a∈Ω + B a (x) 0 1 -1 0 (x a -x Γs ) + ξ Γs N • n δ Γs +α m1 - a∈Ω + B a (x) (x a -x Γs ) • m m + MmH Γs + (x -x Γs ) • m Nmδ Γs (2.12)
where ξ Γs = (x -x Γs ) • m (x ∈ Γ s ) represents the curvilinear abscissa along the discontinuity surface Γ s (which is here a straight line),

N =   n x 0 0 n y n y n x   (M is
identically defined in terms of the components of m) and δ Γs represents the Dirac distribution on the discontinuity surface Γ s .

A plot of the incompatible shape functions M (see equation (2.11) for both constant and linear normal modes (n 0 and n 1 )is presented on figure 2.7. It has to be noticed that the shape functions vanish on the boundaries of the quadrangular domain keeping the contribution of the displacement jump reduced to the corresponding element. The displacement jump ū introduced in section 2.1.1 and the current vector α of separation modes are linked by the following relation:

ū • n ū • m = 1 ξ Γs 0 0 0 0 1 ξ Γs     α n 0 α n 1 α m 0 α m 1     (2.13)
where ξ Γs is the curvilinear abscissa along the discontinuity such that ξ Γs ∈ [-Γs 2 , Γs 2 ], Γs is the discontinuity length, ū • n and ū • m are respectively the normal part and tangential part of the displacement jump on the discontinuity surface.

The strain interpolation can be put in the general form: with

ε(x, t) = Bu + Ḡr α + N. Ḡr αδ Γs where α = < α n 0 , α n 1 , α m 0 , α m 1 > T Ḡr =
                           Ḡr n0 = - a∈Ω + B a (x).n Ḡr m0 = - a∈Ω + B a (x).m Ḡr n1 = - a∈Ω + B a (x) 0 1 -1 0 (x a -x Γs ) Ḡr m1 = - a∈Ω + B a (x) ((x a -x Γs ) • m) m + MmH Γs (2.14)
and

           Ḡr n0 = n Ḡr m0 = m Ḡr n1 = ξ Γs .n Ḡr m1 = ξ Γs .m (2.15)
The virtual strain is interpolated in the same form:

δε(x, t) = Bδu + Ḡv δα + N. Ḡv δαδ Γs (2.16)
where δu and δα denote the virtual nodal displacement and virtual displacement jump respectively.

Ḡv and Ḡv are computed from Ḡr and Ḡr by imposing the patch-test [START_REF] Ibrahimbegovic | A modified method of incompatible modes[END_REF] guaranteed by:

G v = G r - 1 A Ωe Ωe G r dΩ (2.17)
where A Ωe = Ωe dΩ is the area of the considered element Ω e .

Taking into account (2.14), (2.15), and (2.17), we have then:

                         Ḡv n0 = Ḡr n0 - 1 A Ω e Ω e Ḡr n0 dΩ - l Γs A Ω e N.n Ḡv m0 = Ḡr m0 - 1 A Ω e Ω e Ḡr m0 dΩ - l Γs A Ω e N.m Ḡv n1 = Ḡr n1 - 1 A Ω e Ω e Ḡr n1 dΩ Ḡv m1 = Ḡr m1 - 1 A Ω e Ω e Ḡr m1 dΩ (2.18)
and

           Ḡv n0 = n Ḡv m0 = m Ḡv n1 = ξ Γs .n Ḡv m1 = ξ Γs .m (2.19)
Remark 2.1.3. The total displacement field introduced in (2.9) and (2.11) are formed only to deduce the interpolation of the deformation fields. By using the incompatible modes method, only the deformation field is modified.

Problem resolution

Having presented the interpolation used for displacements and strains, we introduce in this section the numerical resolution procedure used in this work.

Since we introduce an incompatible strain field, the three-field Hu-Washizu principle is an appropriate formulation of the problem. A two-field formulation could also be suitable, however using the three-field Hu-Washizu principle allows to rewrite the equations in the most general framework possible. The following expression of the non-linear discrete problem comes from developments established by [START_REF] Wilson | Use of incompatible displacement modes for the calculation of element stiffnesses or stresses[END_REF] and [START_REF] Ibrahimbegovic | A modified method of incompatible modes[END_REF] in the framework of the so-called "incompatible modes method". Detailed developments that lead to the further formulation can be found in appendix A.

Find (d, α) such that:

     N el A e=1 [f int,e (x) -f ext,e (x)] = 0 h e (x) = 0, ∀e ∈ [1, N el ] (2.20a) (2.20b) with                        f int,e (x) = Ω e B T (x)σ(ε(x)) dΩ f ext,e (x) = Ω e N T b(x) dΩ + ∂tΩ e N T h(x) dΓ h e (x) = Ω e \Γs ḠT v (x)σ(ε(x))dΩ + Γs Ḡv T (x)t Γs (x)dΓ ε | Ω e (x) = B(x)d + G r (x)α (2.21a) (2.21b) (2.21c) (2.21d)
where d are the nodal displacements, α are the displacement jumps projected along the different modes described in section 2.1.1, and N el is the number of elements. N(x) = [N 1 (x)...N Nn (x)] is the matrix of shape functions for a standard Q4 element, B(x) = [B 1 (x)...B Nn (x)] = LN(x) with L the matrix form of the ∇ s operator, and δu T (x, t) = [δu 1 (x, t)...δu Nn (x, t)], b(x) is the imposed volumic force field and h(x) the external force vector imposed on surface ∂ t Ω e . Ḡv and Ḡv are defined in (2.18) and (2.19).

σ(ε(x)) is the stress tensor in the bulk defined from the continuous part of the strain tensor ε(x). In the following work, the bulk is considered to have a non-linear behavior, considering both plasticity and damage. For further information on the material model, please refer to section 2.2.1 in which the model used for the bulk in this work is detailed.

t Γs (x) refers to the driving traction on the discontinuity surface Γ s . It is linked by a non-linear cohesive law (detailed in section 2.2.2) to the displacement jump ū (and so the separation modes α defined in section 2.1.1 above).

We note that both the bulk and the discontinuity surface account for non-linear behaviors, which is a matter of interest in order to solve the problem defined in (2.20) above. This is the purpose of the following section.

Resolution of the discretized problem

Due to the non-linear behavior of the bulk and discontinuity, system (2.20) is nonlinear. The solution is obtained by linearization and introduction of a pseudo-time "t". Let us consider a global iteration i of a time step n + 1. We denote:

d(t n+1 ) (i) = d (i) n+1 ; α(t n+1 ) (i) = α (i) n+1 (2.22) and ∆d (i) n+1 = d (i+1) n+1 -d (i) n+1
; ∆α

(i) n+1 = α (i+1) n+1 -α (i) n+1 (2.23)
With such notations, the full linearization of system (2.20) can be written as:

                                   N el A e=1 Ω e \Γs B T ∂σ ∂d (i) n+1 dΩ ∆d (i) n+1 + Ω e \Γs B T ∂σ ∂α (i) n+1 dΩ ∆α (i) n+1 = N el A e=1 f ext,e(i) n+1 -f int,e(i) n+1 h e(i) + Ω e \Γs ḠT v ∂σ ∂d (i) n+1 dΩ + Γs Ḡv T ∂t Γs ∂d (i) n+1 dΓ ∆d (i) n+1 + Ω e \Γs ḠT v ∂σ ∂α (i) n+1 dΩ + Γs Ḡv T ∂t Γs ∂α (i) n+1 dΓ ∆α (i) n+1 = 0 (2.24)
Furthermore, constitutive laws for both the bulk and the discontinuity are written in an incremental form as:

   ∆σ (i) n+1 = C an (i) n+1 B∆d (i) n+1 + Ḡr ∆α (i) n+1 ∆t (i) Γsn+1 = Can (i) n+1 ∆ ū(i) n+1 (2.25a) (2.25b)
where we have considered the numerical integration of the material laws for the bulk and the discontinuity surface (see sections 2.2.1 and 2.2.2). C an (i) n+1 and Can (i) n+1 are the tangent moduli at time n + 1, iteration (i) for the bulk and for the discontinuity, respectively.

We have then:

                               ∂t Γs ∂α (i) n+1 = ∂t Γs ∂ ū (i) n+1 • ∂ ū ∂α = Can (i) n+1 1 ξ Γs 0 0 0 0 1 ξ Γs ∂σ | Ω\Γs ∂α (i) n+1 = C an (i) n+1 Ḡr ∂σ | Ω\Γs ∂d (i) n+1 = C an (i) n+1 B ∂t Γs ∂d (i) n+1 = 0 (2.26a) (2.26b) (2.26c) (2.26d)
Due to the expression of the constitutive law on the discontinuity surface in (2.25b), the values of tractions only depend on the displacement jumps and not on nodal displacements, this is why we obtain equation (2.26d).

With these expressions in (2.26), the system to be solved may be written:

       N el A e=1 K e(i) f d,n+1 ∆d (i) n+1 + K e(i) f α,n+1 ∆α (i) n+1 = N el A e=1 f ext,e(i) n+1 -f int,e(i) n+1 h e(i) n+1 + K e(i) hd,n+1 ∆d (i) n+1 + K e(i) hα,n+1 + F e(i) hα,n+1 ∆α (i) n+1 = 0 (2.27a) (2.27b)
where

K e(i) f d,n+1 = Ω e \Γs B T C an (i) n+1 BdΩ ; K e(i) f α,n+1 = Ω e \Γs B T C an (i) n+1 Ḡr dΩ K e(i) hd,n+1 = Ω e \Γs ḠT v C an (i) n+1 BdΩ ; K e(i) hα,n+1 = Ω e \Γs ḠT v C an (i) n+1 Ḡr dΩ F e(i) hα,n+1 = Γs Ḡv T Can (i) n+1 ∂ ū ∂α dΓ (2.28)
The resolution is then sequential. The nodal displacements are obtained by solving at the global level equation (2.27a) while the separation modes ∆α n+1 . This procedure is called the "operator split method".

First the problem is solved at the element level with a standard Newton-Raphson procedure. The local equation is written for local iteration j as follows (i denotes the global interation):

Knowing d (i) n+1 , for each element e ∈ [1, ..., N e ], find α (i,j) n+1 such that : h e(i,j) n+1 + K e(i,j) hα,n+1 + F e(i,j) hα,n+1 ∆α (i,j) n+1 = 0 (2.29)
At the end of the procedure, we have

h e(i,j) n+1 = 0 (2.30)
And the linearized system (2.27) is written in matrix form:

   K e(i) f d,n+1 K e(i) f α,n+1 K e(i) hd,n+1 K e(i) hα,n+1 + F e(i) hα,n+1       ∆d (i) n+1 ∆α (i) n+1    =    f ext,e(i) n+1 -f int,e(i) n+1 0    (2.31)
Using the static condensation of the second equation in (2.31), we have :

∆α (i) n+1 = -K e(i) hα,n+1 + F e(i) hα,n+1 -1 K e(i) hd,n+1 ∆d (i) n+1
(2.32)

And then the global problem is written:

N el A e=1 K e(i) f d,n+1 -K e(i) f α,n+1 K e(i) hα,n+1 + F e(i) hα,n+1 -1 K e(i) hd,n+1 ∆d (i) n+1 = N el A e=1 f ext,e(i) n+1 -f int,e (i) n+1 (2.33)
which is equivalent to:

N el A e=1 Ke(i) n+1 ∆d (i) n+1 = N el A e=1 f ext,e(i) n+1 -f
int,e(i) n+1

(2.34)

with

Ke(i) n+1 = K e(i) f d,n+1 -K e(i) f α,n+1 K e(i) hα,n+1 + F e(i) hα,n+1 -1 K e(i) hd,n+1 (2.35)
System (2.34) is equivalent to the standard non-linear Finite Element problem, the only difference is the formulation of the element tangent stiffness matrix Ke(i)

n+1
given in (2.35). This matrix is not symmetric which involves an appropriate resolution of system (2.34). The fact that the global Finite Element procedure remains unchanged highlights the main advantage of this method. Only a new element needs to be implemented, the rest of the FE resolution architecture remains unchanged, while the introduction of new degrees of freedom, as for X-FEM for example, involves a specific treatment at the global level.

The full enhanced Finite Element resolution procedure is sumed up in the following algorithm:

Data: d n , α n Result: d n+1 , α n+1 i = 0; while A N el e=1 Ke(i) n+1 ∆d (i) n+1 -A N el e=1 f ext,e n+1 -f int,e n+1
> tol do for e = 1 → N el do ∆σ

(i) n+1 = C an (i) n+1 B∆d (i) n+1 if localization not yet detected then Ke(i) n+1 = K e(i) f d,n+1 else j = 0 α (i,j) n+1 = α n while h e(i,j) n+1 > tol do ∆σ (i,j) n+1 = C an (i,j) n+1 B∆d (i) n+1 + Ḡr ∆α (i,j) n+1 ∆t (i,j) Γsn+1 = Can (i,j) n+1 ∂ ū ∂α ∆α (i,j) n+1 h e(i,j) n+1 = Ω e \Γs ḠT v σ (i,j) n+1 dΩ + Γs Ḡv T t (i,j) Γsn+1 dΓ α (i,j+1) n+1 = α (i,j) n+1 -K e(i,j) hα,n+1 + F e(i,j) hα,n+1 -1 h e(i,j) n+1 j = j + 1 end Ke(i) n+1 = K e(i) f d,n+1 -K e(i) f α,n+1 K e(i) hα,n+1 + F e(i) hα,n+1 -1 K e(i) hd,n+1
end f

int,e(i) n+1

= Ω e B T σ

(i) n+1 dΩ f ext,e(i) n+1 = Ω e N T b (i) n+1 dΩ + ∂tΩ e N T h (i) n+1 dΓ end Solve A N el e=1 Ke(i) n+1 ∆d (i+1) n+1 = A N el e=1 f ext,e(i) n+1 -f int,e(i) n+1 to find d (i+1) n+1 i = i + 1 end
Algorithm 1: Full resolution procedure where ∂ t Ω e is the boundary of the domain where external forces are applied.

In the framework of ductile fracture, the plastic strain tends to prevail over the amount of elastic strain, which implies an almost inscompressible strain field. In a numerical point of view, this can induce locking issues for some cases, especially when using a Q4 element in our case (as detailed in section 2.1.1 above). In the following, we briefly summarize the main ingredients of the B-bar method and comment on the modifications needed to couple such a method to the embedded discontinuity method.

B-bar method implementation along with SDA

B-bar method -main ingredients

The B-bar method was designed as an alternative method to mixed formulation to overcome locking obtained for incompressible or near incompressible problems. The basic idea of the B-bar method is to decompose the strain field into its deviatoric and dilational parts and choose different interpolations for those two parts.

In order to fix the notations, we recall for the case of standard finite element the strain interpolation:

ε(x) = B(x)d (2.36)
where d denotes the nodal displacement vector and with

B(x) = B 1 (x) B 2 (x) . . . B n (x) (2.37)
where n is the number of nodes of the considered problem.

For a 2D problem under plane strain hypothesis, we have:

B i (x) =      B 1 (x) = ∂N i (x) ∂x 1 0 0 B 2 (x) = ∂N i (x) ∂x 1 0 0 B 2 (x) = ∂N i (x) ∂x 2 B 1 (x) = ∂N i (x) ∂x 1      (2.38)
where N i is the interpolation shape function associated to node i.

The dilational part of B i is then obtained by:

B sph i (x) = 1 3     B 1 (x) B 2 (x) B 1 (x) B 2 (x) 0 0 B 1 (x) B 2 (x)     (2.39)
The B-bar method consists in an appropriate definition of the dilatational part of the displacement to strain operator B. Different choices are possible (see [START_REF] Hughes | The Finite Element Methods[END_REF]), we've chosen here to replace the dialational part of B i by its mean value over an element:

Bi (x) = Ω e B i (x)
dΩ

Ω e dΩ (2.40)
We obtain then:

Bi (x) =      2B 1 (x)+ B1 (x)) 3 -B 2 (x)+ B2 (x) 3 -B 1 (x)+ B1 (x)) 3 2B 2 (x)+ B2 (x)) 3 -B 1 (x)+ B1 (x)) 3 -B 2 (x)+ B2 (x) 3 B 2 (x) B 1 (x)      (2.41)
Remark 2.1.4. This method was also implemented for axisymmetric cases. The formulation is very close to the one presented above, the main difference is that all integrals are computed as follows :

Ω e • dΩ = • rdrdθdz (2.42)
where (r, θ, z) are the cylindrical coordinates.

In the axisymmetric case, the problem does not depend on θ and the integral may be written

Ω e • dΩ = 2π
• rdrdz or

Ω e • dΩ = • rdrdz (2.43)
if a unit angle θ is considered.

Another difference with the plain strain case is the formulation of matrix B i . Indeed, since the gradient is now formulated in cylindrical coordinates, the third component of the strain tensor is not zero anymore (which accounts for direction θ). B i is now written :

B i (x) =      B 1 (x) = ∂N i (x) ∂x 1 0 0 B 2 (x) = ∂N i (x) ∂x 2 N i (x) x 1 0 B 2 (x) = ∂N i (x) ∂x 2 B 1 (x) = ∂N i (x) ∂x 1      x 1 being r.
(2.44)

B-bar method -numerical examples

To illustrate the modifications induced in the structure response when considering B-Bar method, we present in the following a couple of numerical results . A comparison with standard finite element is also provided.

A three-point bending test on a notched specimen is studied for two meshes with different characteristic sizes in the "process zone" which are presented along with the boundary conditions in figure 2.8. The material model used in this example is Lemaitre's coupled damage-plasticity model (see section 2.2.1 for details), which falls in the framework of near incompressible models. It justifies the use of the B-bar method as shown in the following results. As it is obvious on figure 2.9, the so-called locking phenomenon occurs for this test, for both mesh sizes. The computation is stopped when convergence stops occuring for each test.

The B-bar method allows to "unlock" the problem. The benefits of the method are obvious for this test as it is shown on figure 2.9. Indeed, the structure's response is softer than with the standard finite element computation, and the maximum displacement is much higher in this case. The convergence rate is also better with the B-bar method, when the damage variable of the most deformed element becomes close to 1.

Embedded discontinuities with B-bar method

In section 2.1.1, we presented the kinematics and the implementation of SDA in a four-node quadrangular element, considering a linear jump displacement, for a standard finite element formulation. We now focus on the modifications induced by the use of the B-bar method in a SDA computation framework.

There are only a few modifications to undertake to obtain a coupled formulation of SDA using a B-bar formulation of the strain tensor. As for a non-SDA formulation, the only change induced by a B-bar formulation is the use of matrix B instead of B in all the computed matrices that depend on matrix B.

Then, the total strain with B-bar coupled with the SDA formulation takes the following form : with

ε(x, t) = Bu + Ḡr α + N. Ḡr αδ Γs where α = < α n 0 , α n 1 , α m 0 , α m 1 > T Ḡr =
                           Ḡr n0 = - a∈Ω + Ba (x).n Ḡr m0 = - a∈Ω + Ba (x).m Ḡr n1 = - a∈Ω + Ba (x) 0 1 -1 0 (x a -x Γs ) Ḡr m1 = - a∈Ω + Ba (x) ((x a -x Γs ) • m) m + MmH Γs (2.46)
and

           Ḡr n0 = n Ḡr m0 = m Ḡr n1 = ξ Γs .n Ḡr m1 = ξ Γs .m (2.47)
Numerical examples with SDA coupled (or not) with B-bar will be presented in section 2.3. We will show the interest of introducing a coupling of SDA with B-bar for the case of the three-point bending test, presented above without SDA.

Material model formulation and implementation

In this section, we present the consitutive laws for both the bulk and the surface of discontinuity involved in this model.

Bulk material model : Lemaitre damage model

With the aim of representing the pre-peak phase of the material behavior, we consider the coupled damage-plasticity Lemaitre model [START_REF] Lemaitre | A continuous damage mechanics model for ductile fracture[END_REF] widely used in the literature to model ductile metal failure. We limit the study to the consideration of isotropic damage (the scalar variable D is the isotropic damage) and isotropic hardening ( ξ denotes the accumulated plastic strain and q its dual variable). In this case, the yield function is very close to the von Mises yield function taking into account the effective stress σ = σ/(1 -D), and is written as follows :

φ(σ, D, q) = 3 2 dev( σ) : dev( σ) -σ y 0 + q( ξ) (2.48)
where σ y 0 is the elastic limit. q is the hardening variable associated to the accumulated plastic strain ξ defined in (2.49), and dev( σ) is the deviatoric part of the effective stress tensor.

q( ξ) = (σ y 0 -σ y inf ) 1 -exp(β ξ) (2.49)
In the Lemaitre damage model, the damage evolution law is written as follows:

Ḋ = γ 1 -D -Y r s (2.50)
where Y is the dual variable of the damage variable D (also called damage energy release rate), γ is the plastic multiplier, and s and r are the two Lemaitre parameters that characterize damage evolution.

The consistent tangent operator that appears in equation (2.25) takes the following form:

C an n+1 = dσ n+1 dε n+1 = 1 1 -D n+1 a I -1 3 1 ⊗ 1 + b sn+1 ⊗ sn+1 + c sn+1 ⊗ 1 + d 1 ⊗ sn+1 + e 1 ⊗ 1 (2.51)
where sn+1 = dev σ n+1 dev σ n+1 , I is the fourth-order identity tensor, 1 is the secondorder identity tensor.

Expressions of a, b, c and d along with details on the numerical implementation of this model can be found in appendix B.

Remark 2.2.1. It has to be noticed that we use a fully implicit integration scheme for the elastoplastic material coupled with damage, even though other techniques are used in the literature. We can mention the works of Bergheau et al. [START_REF] Bergheau | Intégration numérique de lois de comportement élastoviscoplastique endommageable et applications[END_REF] who use a semi-implicit integration scheme over time for the elastoplastic damageable material. In their works, the authors take into account interactions between creep phenomena and damage of the material in the context of fracture of pressure vessels. In particular, they compared the efficiencies of several tangent stiffness matrices.

To properly define the discontinuity, several ingredients are needed. We have to define a material model for the discontinuity surface (a cohesive law). We also need to define a criterion to introduce the discontinuity since it is not taken into account at the beginning of the computation. Finally, an orientation of the discontinuity has to be chosen. Note that it is considered to remain fixed during the computation. These three aspects are dealt with in the folowing sections.

Cohesive law on the discontinuity surface

Once a discontinuity suface is introduced within an element, a new surfacic dissipative phenomenon occurs and is added to the bulk dissipation. Bulk and surfacic dissipative phenomena can coexist so that damage can still evolve in the bulk while surfacic dissipation occurs.

To model this surfacic dissipation, we choose a damageable interface model, within the framework of interface thermodynamics. In order to be consistent with the perfectly rigid behavior of the discontiniuty surface that is needed before the crack opens, we choose to model the surfacic behavior in terms of the surfacic compliance, unlike most damageable cohesive models which are defined in terms of the surfacic stiffness and a damage variable (see [START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF], [START_REF] Ortiz | Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis[END_REF], [START_REF] Charlotte | Initiation of cracks with cohesive force models: a variational approach[END_REF], ...).

We choose to use a damage exponential cohesive law (see figure 2.10). No plasticity occurs during the loading, so when the discontinuity surface is completely unloaded, there is no residual displacement jump. As for a damage law for bulk models, the material stiffness decreases with damage evolution, as it is shown on figure 2.10.

The following deals with the thermodynamics of the chosen discrete law. This model uses a surface energy Ψ defined on a discontinuity surface, and depends on the displacement jump between the crack lips ū, the surfacic compliance Q and an internal variable ξ that deals with the cracking irreversibility. The surface energy is given by:

Ψ( ū, Q nn , Q mm , ξn , ξm ) = 1 2 ū.(Q nn n ⊗ n + Q mm m ⊗ m Q ) -1 . ū + Ξ( ξn , ξm ) (2.52)
where n and m are the normal and tangential vectors to the discontinuity, respectively. ( ξn , ξm ) are the internal variable associated to softening in directions n and m, respectively, and their dual variables are ( qn , qm ).

The model takes into account both I and II fracture modes, where Q nn and Q mm are the scalar compliances in directions n and m respectively. Only the terms involving the normal components of ū and t Γs will be defined in the following expressions. Indeed, the cohesive laws are similar in both directions, the only change being that negative jumps are permitted in the tangential direction while it is not in the normal direction (see section 2.2.5 for further details). The dual variables (t n , Ȳn , qn ) associated to the state variable ūn and the internal variables Q nn and ξn , respectively, are given by the state equations:

                   t n = t Γs .n = ∂ Ψ ∂ ūn = ūn Q nn Ȳn = - ∂ Ψ ∂Q nn = 1 2 t 2 n qn = - ∂ Ψ ∂ ξn = σ c [1 -exp(-σc Gc ξn )] (2.53a) (2.53b) (2.53c)
σ c is the limit traction in mode I. It is the maximum principal stress at the onset of localization. When a maximum principal stress criterion is considered, the value of σ c is a user parameter to be set, while this value is computed as the maximum principal stress when a damage criterion is considered (see the following section for further aspects on the introduction criteria).

G c is the fracture energy density of the surface of discontinuity, it corresponds to the energy that is needed to fully open a unit surface of discontinuity.

The proposed model belongs to the framework of non-associated laws, where the yield function φn and the dissipation potential F n , in the case of a mode I fracture, are given by:

   φn (t n , qn ) = |t n | -(σ c -qn ) F n (t n , qn ) = |t n | 2 -ln |tn| σc -(σ c -qn ) (2.54a) (2.54b)
The evolution of the state variables are given in terms of the dissipation potential as:

         ξn = γn ∂F n ∂ qn = γn Qnn = γn ∂F n ∂ Ȳn = γn t n ∂F n ∂t n (2.55a) (2.55b)
where ξn is the Lagrange multiplier for the opening mode in direction n. Together with the Kuhn-Tucker conditions, if φn = 0 then γn > 0 and we have:

t n = σ c exp - σ c G c ξn (2.56)
And finally, equation (2.55b) may be written:

Qnn = γn t n 1 + σ c G c ξn (2.57)
Integration of equation (2.57) gives :

Q nn (t) = t 0 Qnn dt = ξ(t) 0 1 σ c exp σ c G c ξn 1 + σ c G c ξn d ξ = ξn σ c exp σ c G c ξn (2.58)
And finally, the consistent tangent operator, for both linear loading/unloading and dissipative phases, is given by:

Can nn = ∂t n ∂ ūn =        σc ξn exp -σc Gc ξn if φn < 0 σc ξn 1 -( ūn ) 2 σc Gc + 1 ξn exp -σc Gc ξn if φn = 0 (2.59)
Remark 2.2.2. Considering that the normal and tangent behaviors are uncoupled, the consistent tangent modulus corresponding to the tangent behavior Can mm takes the same form as the one corresponding to the normal behavior ( Can nn in (2.59) above). These two consistent tangent moduli are used in the resolution of the local equilibrium equation (2.27) with the element matrices defined in (2.28).

Criteria to introduce the surface of discontinuity

As specified in section 2.2.2, the formulation allows to implement, in a same framework, two different possibilities for the criteria of introduction of the discontinuity:

• A maximum principal stress criterion: the average of the maximum principal stress over the Gauss points of the bulk in the current element is computed at each time step. If this averaged value exceeds the critical value specified as a parameter, then this current element is considered to be localized and a local iterative procedure is undertaken.

standard element if σ < σ c localized element if σ ≥ σ c (2.60)
where σ represents the average of the maximum principal stress over the four Gauss points, and σ c is the critical maximum principal stress.

• A damage criterion: the same procedure as previosuly presented is considered, we compute the average damage value over the Gauss points of the bulk in the current element at each time step, and if this value exceeds the specified critical damage value, the element is localized.

standard element if D < D c localized element if D ≥ D c (2.61)
Differences between the two criteria are discussed in section 2.3.4.

Remark 2.2.3. The integration of the discrete material model requires the information of the trial tractions on the discontinuity. This does not imply any issue if the discontinuity surface is already opened (see the state equation (2.53a)). However, at the onset of localization, that is when the criterion has just been verified, this information is not available and we use the following procedure to compute the trial tractions.

In order to compute these trial tractions in the case of rigid beahavior of the discontinuity (at the onset of localization), we use equation (2.20b) which has to be verified for each localized element, in particular at the very moment when the discontinuity is introduced. Projecting equation (2.20b) on the different modes of opening of the discontinuity < n 0 , n 1 , m 0 , m 1 >, we obtain a set of four equations:

                                                 h e n 0 = h e Ω,n0 Ω e ḠT v n0 σdΩ e + Γs ḠT v n0 t Γs dΓ = 0 h e n 1 = h e Ω,n1 Ω e ḠT v n1 σdΩ e + Γs ḠT v n0 ξ Γ t Γs dΓ = 0 h e m 0 = h e Ω,m0 Ω e ḠT v m0 σdΩ e + Γs ḠT v m0 t Γs dΓ = 0 h e m 1 = h e Ω,m1 Ω e ḠT v m0 σdΩ e + Γs ḠT v m1 ξ Γ t Γs dΓ = 0 (2.62a) (2.62b) (2.62c) (2.62d)
When computing Γs • dΓ by a two points Gauss integration technique, we obtain:

                         h e Ω,n0 + l Γs 2 t 1 n,Γs + t 2 n,Γs = 0 h e Ω,n1 + l Γs 2 - l Γs 2 √ 3 t 1 n,Γs + l Γs 2 √ 3 t 2 n,Γs = 0 h e Ω,m0 + l Γs 2 t 1 m,Γs + t 2 m,Γs = 0 h e Ω,m1 + l Γs 2 - l Γs 2 √ 3 t 1 m,Γs + l Γs 2 √ 3 t 2 m,Γs = 0 (2.63a) (2.63b) (2.63c) (2.63d)
where the Gauss points on the discontinuity are such that ξ Γs = ± l Γs 2 √

3 and where we have denoted as t i n,Γs and t i m,Γs the normal and tangential components of the traction at integration point number i (i = 1, 2).

Orientation of the discontinuity surface

The orientation of the newly introduced discontinuity surface is governed by the type of opening we wish to model. Hence, we have the two following different possibilities for the orientation of the discontinuity surface:

• In the case of a mode I fracture modeling framework, the direction is defined by the direction of maximum principal stress, whether a maximum principal stress or a damage criterion is considered. In particular, the normal to the discontinuity is defined as an average over the Gauss points in the bulk of the current element.

• In the case of a mode II fracture modeling framework, the direction of the discontinuity surface is defined by the direction of maximum shear stress. In the following example dealing with mode II fracture (see section 2.3.5) a damage criterion is used. In this case, the tangent to the discontinuity is oriented along the average of the directions of maximum shear stress over the Gauss points of the element, that is 45 • to the maximum principal stress direction.

Remark 2.2.4. As for the position of the discontinuity surface, we introduce it at the center of the localized element.This consideration is probably a limitation of this work. Indeed, with such a consideration, the continuity of the crack path cannot be ensured, which is not a strong limitation if we consider that the aim of this modeling is to regularize Lemaitre's model in the framework of ductile fracture, and so, to be of concern essentially to the global response of the structure's failure. Authors such as [START_REF] Linder | Finite elements with embedded strong discontinuities for the modeling of failure in solids[END_REF] worked on more flexible positionings of the discontinuity within a fournode element. They also consider branching discontinuities (see [START_REF] Linder | Finite elements with embedded branching[END_REF]) which imply a treatment at global level, this reduces the local aspect of the method.

Interpenetration processing

In fact the normal displacement is not allowed to be negative to prevent interpenetration. Numerically, the element acts as if no localization had occured and equation h e = 0 is considered to be solved. It means that the local iterative procedure is stopped and two cases are to be considered for the computation of the consistent tangent operator Can given in (2.59):

• If there was no displacement jump at the previous time step ( ūn •n = 0) and the discontinuity has never opened ( ξn = 0), which means the discontinuity still has a perfectly rigid behavior, then the consistent tangent operator takes a very large value (∼ 10 12 ) to simulate the perfectly rigid behavior of the discontinuity (penalty method to treat rigid behavior at initiation of localization).

• If the discontinuity has already opened normally at least once in the loading history ( ξn > 0) then, the consistent tangent operator takes the value of the linear case specified in equation (2.59) (with ξn+1 = ξn ).

In both cases, the stress is updated considering no displacement jump.

Numerical results and discussion

In order to test the capabilities of the proposed strategy to regularize Lemaitre damage model, three numerical tests are considered: a simple traction test, a threepoint bending test and a tensile test on a cylindrical notched specimen.

As specified previously, we use in the following examples a B-bar extension of the method (see section 2.1.3 for the details of the implementation).

Uniaxial tensile test

Results presented in this section deal with a uniaxial tensile test of a 10 cm x 5 cm x 5 cm plate under plane strain hypothesis (see figure 2.11). A displacement is incrementally prescribed on the right side of the specimen. Parameters for Lemaitre's model were calibrated for an AISI 1010 low carbon steel in a rolled state by Benallal et al. [START_REF] Benallal | Crack initiation and propagation analyses taking into account initial strain hardening and damage fields[END_REF]. These parameters are presented in table 2.1, and the studied meshes along with the normal opening at the end of the test are presented in figure 2.12 An elastic unloading is carried out around 17 mm of prescribed displacement and we note that damage already occurred in the bulk (See figures 2.13(a) and 2.13(b)). The discontinuity surface was introduced when the maximal principal stress reached σ d = 920 MPa for every element, whereas an element located at the bottom of the specimen was weakened with σ d = 919 MPa in order to control the onset of the failure process. Then the discontinuity surface is propagated straight from this weakened element and orthogonally to the loading direction, which means a mode I type of behavior was considered even though it is well-known that such a test involves a mode II fracture mode. This test was performed to attest the regularizing capability of the method and not the precise simulation of ductile fracture, thus it should only be considered as a first step towards fracture modeling.

As it is seen on figure 2.13(a), the global softening behavior is quite brittle which involves that the global softening phase begins when the first element is localized. Besides, it is shown that the solution remains mesh independent. In figure 2.13(b), the damage evolution in the weakened element is presented, damage evolves until it reaches a value of 0.13 that corresponds to the onset of localization in the element. Thus, for this test, the volumic and surfacic mechanisms are never activated at the same time within an element.

Embedded discontinuities with the B-bar method

In order to attest the regularization capabilities of the SDA method combined with a B-bar formulation, let us recall the three-point bending test to which the mesh dependency issue was demonstrated above in section 2.1.3. It was shown that fracture (which means the most deformed element has reached a damage value of 1) occured sooner for the fine mesh, in terms of maximum displacement, than for the fine mesh, in quite large proportions. This is a well-knwon result when using a softening material behavior such as Lemaitre's model.

Three different meshes are studied, the element size differ in the cracking area such that there are two elements in front of the notch for the "coarse" mesh while there are five elements in this area for the "fine" mesh and four elements for the "intermediate" mesh.

The same test with the same meshes was carried out, but this time taking into account the surfacic mechanisms. The criterion of introduction of the discontinuity surface is a damage criterion. When the average of the damage variable within an element reaches 0.1, the critical traction σ c is computed as the average of maximum principal stresses over the Gauss points of the element at this very moment of localization. As the model takes into account a constant fracture energy G c , the dissipated energy is the same for every opening discontinuity, whatever this critical traction value is. Hence, the regularization is ensured.

The SDA regularization is undertaken for this test using a fracture energy quite large (G c = 1000 MPa.mm), enabling the material to remain ductile for both the bulk and the surfacic beahavior, which means both dissipative mechanisms can be activated at the same time. Only mode I fracture is considered for this test. Results in terms of reaction/displacement at nodes where the displacement is applied are shown in figure 2.14.

It is obvious that results with SDA are much less mesh-dependent for this test. The benefits of the use of this method for this test are the equivalence of the computed ultimate loads for both meshes, along with the similar prediction of both the onset of localization and dissipated energy, while they highly differ without the use of SDA. Besides, there is a good convergence rate even in the global post-peak phase, while it induces bad convergence rates without SDA when passing the global peak.

Even though these results are not completely mesh-independent while the propagation goes further, it is necessary to keep in mind that the considered meshes have characteristic lenghts in the process zone that are very different. Besides, even though the "fine" mesh curve and the "coarse" mesh curves are not totally superposed, the "intermediate" mesh curve shows that results tend to converge to a response very close to the "fine" one in figure 2 The model and the corresponding meshes are presented in figure 2.15 along with the crack opening at the end. Results in terms of force/displacement are presented in figure 2.16 (for a unit thickness). The influence of the activation of the linear mode on the global response has been studied for two topologically identical meshes, but with different element sizes. The relatively coarse mesh considers 1015 Q4 elements while the fine mesh considers 4060 elements.

We note that results differ slightly when we take into consideration the linear mode, and that the peak is higher when only a constant displacement jump is considered. This may be due to the well-known "locking phenomenon" that occurs in the case of the constant jump alone. The benefits of the consideration of a linear opening mode to avoid this stress-locking were enlightened by [START_REF] Manzoli | A general technique to embed non-uniform discontinuities into standard solid finite elements[END_REF] in the case of bending tests. The interest in the linear mode is obvious for this test, since we observe that a stress concentration occurs near the notch that induces an heterogeneous stress state in the neighboring elements. Thus, a linear opening of the crack allows for a more correct representation of the traction state along the interface for these elements than with only a simple constant crack opening. mesh a mesh b deformed mesh a (x20)

Damage criterion

The same test was performed, this time using a damage criterion in order to show the versatility of the method with respect to the criterion used to introduce the discontinuity surface.

Three different meshes are studied, the element size differs in the cracking area such that there are two elements in front of the notch for the "coarse" mesh while there are five elements in this area for the "fine" mesh and four elements for the "intermediate" mesh.

A first set of tests was undertaken to show the mesh dependency if no SDA regularization is used, meaning only Lemaitre's damage model is responsible for the global softening behavior. Results show that for a smaller mesh size, convergence stops occuring sooner in the test than with the larger mesh, and the peak appears also sooner. This is due to damage localization occuring faster for smaller elements in the process zone.

Results in terms of reaction/displacement at the top of the specimen (where the displacement is prescribed) along with the considered meshes (coarse and fine) and boundary conditions are shown on figure 2.18 and figure 2.17, respectively. The same test with the same meshes was carried out, but this time taking into account the surfacic mechanisms. The criterion of introduction of the discontinuity surface is a damage criterion. When the average of the damage variable within an element reaches D c , the critical traction σ c (in the traction/displacement jump law) is computed as the average of maximum principal stresses over the Gauss points of As the model takes into account a constant fracture energy G c , the dissipated energy is the same for every opening discontinuity, whatever this critical traction value is. Hence, the regularization is ensured.

In order to preclude the loss of uniqueness at the material point level, D c was set to 0.1 which is a quite low value that ensures that the standard localization condition is not yet reached. However we plan to use the standard localization criterion within this framework in future works. This first implementation should be considered to be as a first step to localization prediction The SDA regularization is undertaken for this test using a fracture energy quite large (G c = 1000 M P a.mm), enabling the material to remain ductile for both the bulk and the surfacic beahavior, which means both dissipative mechanisms can be activated at the same time. Only mode I fracture is considered for this test. Results in terms of reaction/displacement at nodes where the displacement is applied are shown in figure 2.19. It is obvious that results with SDA are much less mesh-dependent for this test. The benefits of the use of this method for this test are the equivalence of the computed ultimate loads for both meshes, along with the similar prediction of both the onset of localization and dissipated energy, while they highly differ without the use of SDA. Besides, there is a good convergence rate even in the global post-peak phase, while it induces bad convergence rates without SDA when passing the global peak.

Even though these results are not completely mesh-independent while the propagation goes further, it is necessary to keep in mind that the considered meshes have characteristic lengths in the process zone that are very different. Besides, even though the "fine" mesh curve and the "coarse" mesh curves are not totally superposed, the "intermediate" mesh curve shows that results tend to converge to a response very close to the "fine" one in figure 2.19.

Tensile test on a notched bar

A tensile test on a cylindrical notched specimen is performed in this section, considering axisymmetric strain hypothesis. An adjustment of the numerical implementation formulation was carried out to take into account such an hypothesis. In particular, the incompatible shape functions appearing in (2.14), (2.15), (2.18) and (2.19), were adjusted . Adjustments in the computations of the enhanced elementary tangent stiffness Ke(i) n+1 and the residual f int,e(i) n+1

were also undertaken. The same set of parameters used in the previous test is considered for this test. A damage criterion (D c = 0.1) was used considering two different mesh sizes, as shown in figure 2.20. The fracture energy density is G c = 100 M P a.mm, which corresponds to a more brittle behavior than the one studied in 2.3.3. For this set of parameters, even though the accumulated plastic strain is maximum near the notch (see figure 2.22(c)), the dependence of Lemaitre's damage model on the triaxiality ratio (see figure 2.22(b)) allows the damage (see figure 2.22(a)) to be maximum at the center, and so the cracking to be initiated at this location. Then the crack is propagated straight to the notch. tion of the damage variable as well as the mode I crack opening in the first localized element, that is the element at the center of the specimen. It is shown that damage is still evolving while the discontinuity begins opening, and that its evolution is progressively slowing down with respect to the crack opening. This phenomenon does not occur with the uniaxial tensile test presented in section 2.3.1, in which damage stops evolving exactly when the damage threshold is reached, which means that the bulk is elastically unloading at this moment.

Another test was undertaken with a different set of parameters. A fracture energy density of G c = 200 M P a.mm and a damage criterion of D c = 0.1 was used and this time it is shown on figure 2.25(a) (where the crack opening is represented) that there is a bifurcation in the propagation path, which corresponds to the wellknown "cup cone" propagation phenomenon. This result was numerically studied first by Tvergaard and Needleman (see [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]) with the GTN model, and a detailed study of this phenomenon was presented by Besson et al. [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] comparing the prediction capabilities of the models proposed by Gurson and Rousselier by the mean of localization analyses.

This result illustrates the link between the fracture mode and the surfacic model's parameters, more precisely, it is shown that higher values of the damage criterion and the fracture energy (which means higher ductility) induce the onset of the socalled cup cone. Additionally, this is a first step towards a sensitivity analysis of these parameters on the cup cone appearance, and also towards an identification of these parameters for a specific material, as we know that these parameters can be effectively identified.

Shear test

A final test, exhibiting mode II fracture is performed on a notched specimen. Even though this test is originally to be performed in a 3D framework, it has been extended to the 2D case in order to show the capabilities of the 2D enhanced model presented above in a mode II fracture framework. The geometry of the specimen is shown on figure 2.26 with the boundary conditions. We assume a plane strain hypothesis and the thickness of the specimen is 20 mm. Several meshes are tested with different mesh sizes, but with the same bias. They are shown in the final deformed configuration on figure 2.27(b) along with the mode II crack opening in the undeformed configuration on figure 2.27(a).

The criterion used for this test is a damage criterion, which is set to D c = 0.12, and the fracture energy density G c is set to 50 M P a.mm. As for the previous tests, this criterion ensures that the localization condition is not reached at the gauss point.

This value guarantees that the loss of ellipticity condition is not reached for this level of damage. Again and as for all the tests presented, any criterion (like the loss of ellipticity criterion) could be used to initiate localization. For this test, the discontinuity surface is introduced with the tangent to the discontinuity in the same direction as the maximum shear direction. This guarantess that the shear stress in the bulk is oriented in the same direction as the tangent traction on the discontinuity, which will enable the mode II opening mode on the discontinuity surface along with shear elastic unloading in the bulk. A well-known result appears in figure (a), indeed, we note that the post-peak responses differ with respect to mesh size, giving different dissipated energies after the global peak is reached, the slope of the softening response increases with the mesh size to a point that the energy becomes almost zero for a very refined mesh, inducing a loss of convergence of the global implicit iterative Newton procedure for even finer meshes. However, we note that with our method (see figure 2.28(b)), this phenomenon does not appear at all, the slopes of the softening responses are equivalent for all mesh sizes, giving equivalent dissipated energies after the peak. We also note that with our method the peak is not reached at the same prescribed displacement for the coarse meshes a and b and the fine meshes c and d, this is due to the fact that the coarse meshes a and b are not converged, even in the pre-peak phase of the test. This implies that the damage field in the 'process zone' is not well represented for these meshes, and thus the damage criterion becomes irrelevant.

On the contrary, results for meshes c and d are perfectly superimposed, which tends to prove the mesh-size objectivity of the proposed method.

Extension to explicit formulation of the FE problem with SDA

In this section, we present the strong discontinuity method's extension to the explicit global integration scheme, while remaining in the field of the quasi-static problems.

The purpose is to overcome convergence issues with the implicit formulation, which can be difficult to deal with considering a strong discontinuity enhancement in which we deal with three levels of non-linearities (at the Gauss point level, the element level and the global level). Another purpose is to reduce computational costs. We will show with an example that the total computational costs can be drastically reduced, if the time increment is carefully chosen, along with the mass density parameter.

The only difference that exists with the implicit scheme is the integration of dynamic effects in the formulation, which does not change the elementary stiffness matrix, but adds acceleration terms in the residue formulation.

The global equilibrium equation is written:

Mü(t) + K(u(t))u(t) = f (t) (2.64)
A centered difference scheme is used to compute the acceleration and velocity terms:

   ü(t) 1 ∆t 2 (u(t + ∆t) -2u(t) + u(t + ∆t)) u(t) 1 ∆t (u(t + ∆t) -u(t -∆t)) (2.65)
where ∆t is the time-increment.

Critical time step

The use of an explicit integration scheme induces a potential instability if the maximum eigenvalue of the amplification matrix A is greater than 1, where A is defined in a recursive form as:

u(t + ∆t) = Au(t) (2.66)
Let us compute the amplification matrix A. From (2.64) and (2.65), we may write:

M ∆t 2 (u(t + ∆t) -2u(t) + u(t -∆t)) + K(u(t))u(t) = f (t) (2.67)
which after some algebraic manipulations gives:

u(t + ∆t) = (2I -∆t 2 M -1 K)u(t) + ∆t 2 M -1 f (t) -u(t -∆t) (2.68)
Thus we have

A = 2I -∆t 2 M -1 K (2.69)
Then the stability condition of this explicit integration scheme may be written:

2 -∆t 2 λ max ≤ 1 (2.70)
where λ max is the maximum eigenvalue of M -1 K. Equation (2.70) is equivalent to the following condition:

∆t ≥ 2 λ max (2.71)
This means that if the time step exceeds this value, the stability of the scheme is not ensured anymore. Thus an automatization of the time increment has to be implemented, which has been done in this work.

Remark 2.4.1. Matrices M and K are the assembled mass and stiffness matrices, respectively. This means that condition (2.71) deals with these global matrices, and that the maximum eigenvalue of M -1 K has to be computed at the global level. This induces prohibitive additional costs. To overcome this issue, a standard approximation has been considered. Condition (2.71) is actually written at the elementary level, which involves an additional computation for each element but dramatically decreases the total cost of the condition computation since we now deal with matrices that have the size of the elementary matrices. The computation of the maximum eigenvalue at the elementary level then becomes acceptable. This approximation is common in industrial FE codes.

In a practical point of view, in this work, condition (2.71) is computed for every element at each time step, and if the condition is violated, a new time increment is computed and the computation of the global equilibrium is started again considering this new time increment.

Remark 2.4.2. In this work the mass matrix is lumped, which gives a diagonal elementary mass matrix, easy to invert in order to compute the maximum eigenvalues of M e -1 K e Remark 2.4.3. We recall that we lie in the framework of quasi-static problems. This means that the dynamic effects induced by the use of the new scheme have to be lowered as much as possible, since we do not lie in the dynamics computation framework. In other terms and in a practical point of view, it means that the mass density parameter is not viewed as a physical parameter, but more as a numerical parameter that enables the total cost to be lowered, as long as the dynamic effects remain low enough. Two parameters are responsible for theses effects, the mass density and the time increment.

If the mass density increases, then the stability condition (2.71) gives larger critical time increments, while increasing the undesirable dynamic effects and vice versa. Increasing the total time of the test enables to lower the dynamic effects. As a consequence, there is no unique time increment/mass density ideal couple that gives the desired computational cost. A choice has to be made, one can choose to take the physical mass density of the specimen and then increase the total time of the test, or one can do the other way around: choose a total time and then decrease the mass density to lower the dynamic effects.

Comparison between explicit and implicit formulation

The tensile test on a round notched axisymmetric specimen was performed comparing the explicit and implicit formulations. The material is an AISI 120 steel, with E = 210000 M P a and with Poisson's ratio ν = 0.3 and the constitutive law used for the bulk is a von Mises plasticity law with an exponential hardening law of the form σ y = 620 + 3300 * (1 -exp(-βR)) with β = 0.4. Only mode I opening is taken into account for this test. The discrete law on the discontinuity surface is an exponential rigid-damageable law defined with the fracture energy density G c = 200 MPa.mm. The criterion used for the onset of localization (that is the introduction of the discontinuity surface) is a Rice and Tracey criterion with ( R R 0 ) c = 1.3. The mass density for the explicit case is set to the physical mass density, that is ρ = 7800 kg/m 3 , and the total time is set to 200 s for both cases.

The geometry of the specimen along with the mesh and the discontinuity opening at the end of the test are presented in figure 2.29. Remark 2.4.4. For this test, the geometry of the specimen along with the Rice and Tracey criterion used imply that the triaxiality ratio is maximum at the center of the specimen which also imply that the onset of localization appears at the center of the specimen (since the Rice and Tracey criterion strongly depends on the triaxiality ratio), while the maximum accumulated plastic strain remains near the notch. This test has already been performed using Lemaitre plasticity-damage constitutive law As shown in figure 2.30, results are quite similar which proves that the explicit formulation does not exhibit any substantial dynamic effects, with the chosen mass density and total test time.

In table 2.3, the computational costs variables are presented for both formulations.

From table 2.3, it is clear that the interest of the use of an explicit formulation is established. Indeed, the total time computation is 22.4 times smaller with the explicit formulation, and the risk of a loss of convergence is also lowered in this case. However, as explained in remark 2.4.3, the choice of the total time, time increment and mass density parameters are not straightforward and one can easily choose these parameters in a way that the computational cost becomes prohibitive, or in a way that the dynamic effects are no longer insignificant.

Conclusion

In the small hypothesis framework, a regularization method using the strong discontinuity method was presented to deal with the loss of uniqueness of the postpeak phase when dealing with Lemaitre's damage model. The strong discontinuity method allows the consideration of both surfacic and volumic dissipative mechanisms while a standard finite element method only permits to model volumic dissipative mechanisms.

We note that the added incompatible shape functions in localized elements vanish at the elements' nodes, and that the resolution procedure is carried out by an "operator split method", which implies that the local equilibrium is solved sequentially with the global equilibrium resolution. Together with the static condensation of the added degrees of freedom at the global level, no degree of freedom is added in the global equilibrium resolution. Thus, there is no modification carried out in the Newton iterative procedure, unlike in the X-FEM method for example, in which degrees of freedom are added and taken into account at the global level.

The enhanced element presented in this chapter is a quadrangular bilinear element (first presented by [START_REF] Linder | Finite elements with embedded strong discontinuities for the modeling of failure in solids[END_REF]) enhanced with a linear displacement jump which is an extension of the more commonly used three-nodes triangular element that only describe constant strain fields. This new element was never implemented in the context of ductile failure considering a coupled damage-plasticity model in the bulk, and particularly together with Lemaitre's damage model. In the context of ductile failure, only models that consider J2-plasticity theory have been used (see e.g. [START_REF] Armero | New finite elements with embedded strong discontinuities in the finite deformation Range[END_REF]). The traction/displacement jump law used to model the cohesive behavior of the discontinuity surface is a damageable behavior without plasticity consideration.

Lemaitre's damage variable was used as an introduction criterion of the discontinuity surface within a localized element. Both mode I and mode II were considered in the examples presented. Results demonstrating the regularization capabilities of the method were shown in the cases of a three-point bending test, a shear test exhibiting mode II failure mode, and a tensile test on a cylindrical notched bar, in the frameworks of plane strain and axisymmetric computations. Indeed, a good regularization was observed in terms of ultimate load prediction and dissipated energies while conserving good convergence rates. It was shown in the case of the shear test, that it is needed to work with fine enough meshes in order to obtain a good representation of the damage field in the 'process zone' that triggers the introduction of strong discontinuities. This well-known feature is related to standard finite element topics, rather than the consideration of strong discontinuities.

Besides the regularization abilities of the method, results showing different crack paths driven by both the damage criterion and fracture energy density, were also presented. The so-called cup cone phenomenon was also obtained in the case of the axisymmetric computation. A correlation between the cup cone appearance and both the damage criterion and the fracture energy was demonstrated, which is a first step towards an identification of these material parameters and towards precise crack path prediction. We also observe that these fracture parameters are phenomenological parameters that have a physical sense. Nevertheless, further work needs to be carried out to precisely describe the sensitivity of the cup cone occurence to these fracture parameters. This last result lightens Lemaitre's model ability to depict the cup cone phenomenon while this ability has already been attested for Rousselier and Gurson's models [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF].

We recall that there is no resctriction to the implementation of a different bulk model into this strong discontinuity formulation, as long as it is able to give information of the damage state in the bulk or any other variable that can be used as an introduction criterion of the discontinuity surface. In this perspective, further developments and tests will be performed using different ductile damage models for the volumic behavior.

As presented above, this chapter is limited to the small strain hypothesis, in the following chapter a finite strain hypothesis is considered.

Chapter 3

Ductile modeling with the strong discontinuity method in large strain hypothesis

In the previous chapter, the SDA theoretical formulation and its numerical implementation were described in the context of small strain hypothesis. In this chapter, we aim at extending this formulation to the large strain hypothesis in the ductile fracture framework.

SDA in the finite deformation range: theoretical formulation and FE implementation

In the previous chapter considering small strain hypothesis, we emphasized the necessity to enhance the spaces of both the displacement fields and the deformation fields to take into account a strong discontinuity. The following deals with the extension of the previous work to the finite strain framework.

Theoretical formulation

The enhanced deformation gradient formulation

In the following section, a detailed presentation of the theoretical formulation of the enhanced deformation gradient operator is carried out. All entities refering to the reference (material) configuration are designated with upper case letters and the ones refering to the current (spatial) configuration with lower case letters. Since the enhancement of the displacement field is similar to the small strain hypothesis, the developments of the one-dimensional case of the displacement field enhancement is not expressed in this section. For further details, the reader shall refer to section 2.1.1, from equation (2.1) to equation (2.6).

We recall that the transformation from the initial configuration to the deformed configuration defines the following mapping:

(x, t) → x(t) = ϕ(X, t) (3.1)
Considering the enhancement with a displacement jump, this transformation function ϕ is expressed as:

ϕ(X, t) = φ(X, t) + ūH Γs (X) (3.2)
where φ is the standard mapping function from the initial to the deformed current configuration. We recall that it is expressed as:

φ(X, t) = X + u(X, t) (3.3)
Then, we may express the total displacement vector as:

u(X, t) = φ(X, t) -X + ūH Γs (X) = ū(X, t) + ūH Γs (X) (3.4)
In this context, the deformation gradient is written as follows:

{ Figure 3.1: Motion decomposition scheme F(X, t) = F(X, t) + [ ū(X, t) ⊗ N]δ Γs F (3.5)
where F and F represent the regular part and the singular part of the deformation gradient, respectively, and where N is the unit normal vector to Γ s in the initial configuration. Following the notations used in section 2.1.1, F is expressed as:

F(X, t) = I + ∇ū(X, t) = I + ∇ũ(X, t) -∇[ ū(X, t)φ(X)] + ∇ ū(X, t)H Γs (X) (3.6)
where we recall that φ is an additional function that allows to have the same essential boundary conditions for the total displacement field u and ũ. Results on derivation of distributions were used to express F (see [START_REF] Stakgold | Green's functions and boundary value problems[END_REF]) and δ Γs is the Dirac distribution associated to the discontinuity surface. See figure 3.1 for a schematic representation of the proposed motion decomposition.

In (3.5), ū refers to the spatial jump. We may introduce the material jump using the pull-back transformation:

Ū = F-1 ū (3.7)
for which we propose a normal-tangential decomposition as follows:

Ū = α n N + α m M (3.8)
where M is the unit tangential vector to the discontinuity surface Γ s in the undeformed configuration.

With such a decomposition and using equation (3.7), we may write:

ū = α n FN + α m FM = α n n + α m m (3.9)
where

n = FN m = FM (3.10a) (3.10b)
n and m are the convected vectors corresponding to N and M, respectively. It has to be noticed that n and m are not necessary orthogonal nor unitary vectors. Remark 3.1.1. n is not the normal to γ s = ϕ(Γ s ), the transformed of the discontinuity surface Γ s .

One of the normal vectors to γ s is n = F-T N (using Nanson's formula on the transport of a piece of oriented surface).

It follows that n is not unitary either and that n = n . We have the following relations:

n .m = ( F-T N).( FM) = N.( F-1 FM) = N.M = 0 (3.11) because (N, M) is an orthonormal basis.
It follows that m is tangent to the convected surface γ s , remembering that m is not necessarily unitary.

See figure 3.2 in which a graphical representation of bases (n , m ) and (n , m ) is proposed, as well as the decomposition of the displacement jump in basis (n , m ) as it was defined in (3.9).

We may define a multiplicative decomposition of the deformation gradient, as it was proposed by Armero and Garikipati [START_REF] Armero | An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids[END_REF], as:

F = F I + Ū ⊗ Nδ Γs = F F (3.12)
where F = I + Ū ⊗ Nδ Γs . This new motion decomposition is shown on figure 3.3.

Problem definition

In this section, we properly define the problem that is to be linearized and solved using standard non-linear resolution procedures. The formulation in this context of large deformation is similar to the infinitesimal case. The differences lie in the fact that the problem is defined in terms of the first Piola Kirchhoff stress tensor P if we work in the undeformed configuration (Lagrangian formulation), or in terms of the Kirchhoff stress tensor τ if the problem is formulated in the deformed configuration (Eulerian configuration).

Figure 3.3: Multiplicative motion decomposition scheme

Since we introduce an incompatible strain field, the three-field Hu-Washizu principle is an appropriate formulation of the problem. The following expression of the non-linear discrete problem in finite deformation comes from developments established by [START_REF] Wilson | Use of incompatible displacement modes for the calculation of element stiffnesses or stresses[END_REF] and [START_REF] Ibrahimbegovic | A modified method of incompatible modes[END_REF] in the framework of the so-called "incompatible modes method".

Considering w, θ, and η the variations of displacement fields, first Piola-Kirchhoff stress fields and displacement gradient fields, respectively, the formulation states:

Find (u, P, H) such that:

∀(w, θ, η),                Ω ∇w : P dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω θ : [∇u -H] dΩ = 0 Ω η : [P(H) -P] dΩ = 0 (3.13)
where h and b are the prescribed surfacic and body forces, respectively. H is the real deformation field. In this context, P(H) is a second order tensor that satisfies the material law.

With the incompatible modes method, both variations and real displacement gradients (η and H) are decomposed in a standard compatible part (∇w and ∇u) and an incompatible part (η and H):

η = ∇w + η H = ∇u + H (3.14)
With such a decomposition, system (3.13) can be written:

∀(w, θ, η),                Ω ∇w : P dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω θ : H dΩ = 0 Ω ∇w : [P(H) -P] dΩ + Ω η : [P(H) -P] dΩ = 0 (3.15)
The method then states that the space of stress fields (θ and P) and incompatible deformation fields (η and H) are constructed orthogonally.

Hence,

∀(w, η),                Ω ∇w : P(H) dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω ∇w : [P(H) -P] dΩ = 0 Ω η : P(H) dΩ = 0 (3.17)
Considering that the material law is verified exactly, that is P(H) = P, system (3.17) is written:

∀(w, η),          Ω ∇w : P(H) dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω η : P(H) dΩ = 0 (3.18a) (3.18b)
Then, let us consider that the added parts of deformations are decomposed into a regular and a singular part:

η = η + ηδ Γs H = H + Hδ Γs (3.19)
Then, system (3.18) may be written as:

∀(w, η),        Ω ∇w : P(H) dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω\Γs η : P( H) dΩ + Γs η : P | Γs ( H) dΓ = 0 (3.20)
where Γ s is the discontinuity surface and where classical results on distributions have been used.

Besides, the singular part of the deformation fields η and H are supposed to derive from discontinuous displacement fields w and ū respectively, so they may be written as :

η = w ⊗ N H = ū ⊗ N (3.21)
The final formulation of the Hu-Washizu principle in the framework of incompatible modes methods is then written as : ∀(w, η),

         Ω ∇w : P(H) dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω\Γs η : P(H) dΩ + Γs w.[P | Γs (H).N] dΓ = 0 (3.22a) (3.22b)
Alternatively, as presented above, the previous problem may be written in terms of the Kirchhoff stress field τ :

∀(w, η),        Ω ∇w : τ (H) F-T dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω\Γs η : τ (H) F-T dΩ + Γs w.τ | Γs (H) F-T .N dΓ = 0 (3.23)
where we have used the relation P = τ F-T . Introducing n = F-T N, system 3.23 may finally be written:

∀(w, η),        Ω ∇w F-1 : τ (H) dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω\Γs η F-1 : τ (H) dΩ + Γs w.[τ | Γs (H).n ] dΓ = 0 (3.24)
where we have used the following tensorial identity A : (BC) = AC T : B T (A, B and C are second order tensors).

which can be expressed in terms of the convected components of the variations of the displacement gradient and its regular enhanced part as:

∀(w, η),        Ω ∇ x w : τ (H) dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω\Γs ηx : τ (H) dΩ + Γs w.[τ | Γs (H).n ] dΓ = 0 (3.25) with ∇ x w = ∇w F-1 and ηx = η F-1 .
Working in the deformed configuration and with Kirchhoff stresses allows to use the symmetries of the Khirchhoff tensor field, which implies that:

∀(w, η), ∇ x w : τ (H) = ∇ s x w : τ (H) ηx : τ (H) = ηs x : τ (H) (3.26)
where we recall that ∇ s = 1 2 (∇ + ∇ T ). This consideration will notably simplify the finite element formulation of the problem.

Another reason that we choose to work with quantities referring to the deformed configuration is that the chosen constitutive law (which is used for the bulk model in this work, and will be expressed further) deals with Kirchhoff stresses and with the spatial tangent modulus.

The final continuum formulation, that will be the basis of the finite element problem formulation, is expressed as:

∀(w, η),          Ω ∇ s x w : τ (H) dΩ - Ω w.b dΩ - ∂tΩ w.h dΓ = 0 Ω\Γs ηs x : τ (H) dΩ + Γs w.[τ | Γs (H).n ] dΓ = 0 (3.27a) (3.27b)
We may now introduce the driving traction t Γs on the discontinuity surface which is defined as:

t Γs = P | Γs .N = τ | Γs .n = t n n + t m m (3.28)
for which we have projected this traction on the convected basis (n , m ). And equation (3.27b) can be rewritten:

Ω\Γs ηs x : τ (H) dΩ + Γs w.t Γs dΓ = 0 (3.29)

Finite element implementation

The following section deals with the discretization of the continuum problem defined above.

The three-node triangular element enhancement

As specified above, we consider a triangular element with three nodes and one regular integration point. Given that strain and stress fields are constant over the element, the enhanced mode shall also be interpolated as a constant displacement jump over the element to ensure consistency.

Starting from equation (3.4), an interpolation of the displacement can be written as (see e.g. [START_REF] Dvorkin | Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions[END_REF][START_REF] Simo | An analysis of strong discontinuity induced by strain softening solutions in rate-independent solids[END_REF]): The total displacement field is interpolated with a standard part, corresponding to the standard nodal interpolation, and with additional non-standard shape functions corresponding to the added separation modes discussed in 2.1.1, in the following form, if we consider the decomposition ū = α n n + α m m on (n , m ) for the spatial jump ū:

u h (X) = a=3 a=1 N a (X)d a + H Γs (X) - a=3 a=1 N a (X)] ū (3.30) 1 2 3 X Γ s α n N M α m
u h (X) = a=3 a=1 N a (X)d a +α n H Γs (X)n - a∈Ω + N a (X)n +α m H Γs (X)m - a∈Ω + N a (X)m (3.31)
where N a is the standard shape function associated to node "a", and d a are standard nodal displacements for node "a".

The non-conforming shape function appearing in (3.31) is plotted in figure 3.5.

u h (X) = a=4 a=1 N a (X)d a + M(X)[n m ]α (3.32)
where M is a discretized version of the continuous function H Γs (x) -φ(x) introduced in equation (2.4) in section 2.1.1. It is defined as:

M(X) = H Γs (X) - a∈Ω + N a (X) (3.33)
A representation of M is given on figure 3.5. α is a vector that contains the components of the spatial jump, following the decomposition proposed in (3.31). It is expressed as

α = [α n α m ] T .
This function allows the contribution of the separation modes to be limited to a localized element, Ω d is then chosen equal to the considered triangular element. This guarantees the local feature of the proposed method and no additional nodal degrees of freedom have to be introduced like for X-FEM for instance [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF] [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF].

On the basis of this displacement field, the real deformation gradient interpolation is given by:

F h (x) = I h + a=3 a=1 B ns a (X)d a +α n - a∈Ω + B ns a (X).n + N.n δ Γs +α m - a∈Ω + B ns a (X).m + N.m δ Γs (3.34)
where

I h =     1 1 0 0     is the second order identity tensor in vector form, N =    N X 0 0 N Y N Y 0 0 N X   
where N X and N Y are the components of normal vector N (in the undeformed

configuration), B ns a (X) =     ∂Na ∂X 0 0 ∂Na ∂Y ∂Na ∂Y 0 0 ∂Na ∂X   
 is the unsymmetric gradient matrix referring to the undeformed configuration, and δ Γs represents the Dirac distribution on the discontinuity surface Γ s .

The deformation gradient interpolation can be put in the general form:

F h (X, t) = I h + B ns d + Ḡr α + Ḡr αδ Γs where α = < α n , α m > T Ḡr = Ḡrn Ḡrm and Ḡr = Ḡrn Ḡrm with        Ḡrn = - a∈Ω + B ns a (X).n Ḡrm = - a∈Ω + B ns a (X).m (3.35)
and

Ḡrn = N.n Ḡrm = N.m (3.36)
Noticing that we have ū = α n n + α m m , then we can rewrite the deformation gradient in (3.35) as:

F h (X, t) = I h + B ns d + Ḡ r + Ḡr δ Γs ū (3.37)
where

Ḡ r = - a∈Ω + B ns a (X) and Ḡr = N (3.38)
This decomposition can be related to equation (3.31) in which an interpolation of the displacement is defined in the same manner. And we notice that these shape functions do not depend on the current deformation state.

The virtual deformation gradient is interpolated in the same form:

δF(x, t) = I h + B ns δd + Ḡv δα + Ḡv δαδ Γs (3.39)
where δu and δα denote the virtual nodal displacement and virtual displacement jump respectively. Then, the discretization of the virtual displacement gradient η introduced in (3.14) takes the following form:

η = B ns δd + Ḡv δα + Ḡv δαδ Γs (3.40)
for which the enhanced part introduced in (3.14) is defined as: η = Ḡv δα + Ḡv δαδ Γs (3.41) Ḡv and Ḡv are computed from Ḡr and Ḡr by imposing the patch-test [START_REF] Ibrahimbegovic | A modified method of incompatible modes[END_REF] guaranteed by:

G v = G r - 1 A Ωe Ωe G r dΩ (3.42)
where A Ωe = Ωe dΩ is the area of the considered element Ω e in the initial configuration. This allows equation (3.18b) (or equation (3.22b)) to be verified for a constant Piola-Kirchhoff stress field P.

Taking into account (3.35), (3.36), (3.42) and the fact that the enhanced shape functions are constant over the element, we have then:

               Ḡvn = Ḡrn - 1 A Ω e Ω e Ḡrn dΩ =0 - l Γs A Ω e N.n = - l Γs A Ω e N.n Ḡvm = Ḡrm - 1 A Ω e Ω e Ḡrm dΩ =0 - l Γs A Ω e N.m = - l Γs A Ω e N.m (3.43)
and

Ḡvn = N.n Ḡvm = N.m (3.44)
It has to be noted that, with these shape functions at hand for a T3 element, equations (3.22) and (3.27b) give very simple expressions since the different fields involved in the formulation remain constant for such an element.

We recall that the Hu-Washizu formulation may be expressed in terms of the first Piola-Kirchhoff stress tensor P in (3.22), or in terms of the Kirchhoff stress tensor τ in equation (3.27). The virtual shape functions defined above are the ones that are supposed to be used with P.

Considering that the constitutive law for the bulk is, in our case, given in terms of the Kirchhoff tensor stress, we need to express the virtual shape functions that appear in the spatial version of the local equation (3.27) as:

       ḡvn = - l Γs A Ω e n .n ḡvm = - l Γs A Ω e n .m (3.45)
and

ḡvn = n .n ḡvm = n .m (3.46)
where

n =   n x 0 0 n y n y n x  
, for which we have used the fact that the Kirchhoff stress tensor is symmetric. n x and n y are the components of n defined as n = F-T N in remark 3.1.1.

Finite element problem formulation

With the interpolation functions defined above in section 3.1.2, the following expressions are verified:

         Ω e \Γs ḡT v τ dΩ = Ω e \Γs ḠT v P dΩ Γs ḡv T τ | Γs dΓ = Γs Ḡv T P | Γs dΓ (3.47a) (3.47b) with ḡv = [ḡ vn ḡvm ], Ḡv = [ Ḡvn Ḡvm ], ḡv = [ ḡvn ḡvm ] and Ḡv = [ Ḡvn Ḡvm ].
We can easily show that, from definition (3.46) we have the following relation:

ḡv T τ | Γs = [n m ] T τ | Γs .n = [n m ] T t Γs = t n t m (3.48)
where we have used the definition of the driving traction and its projection on (n , m ), as it was introduced in (3.28).

Indeed we have the following relations:

n .t Γs = t n n .n + t m n .m = t n (3.49)
since, as it was specified in remark 3.1.1, n .m = 0, and n .n = 1.

In a similar manner, we can show that m .t Γs = t m . Finally, we have:

Γs ḡv T τ | Γs dΓ = l Γs t n t m (3.50)
Using the same procedure, we show that :

Ω e \Γs ḡT v τ dΓ = l Γs A Ω e Ω e \Γs [n m ] T τ .n dΓ = l Γs [n m ] T τ .n (3.51)
since in a T3 element, the stress field τ remains constant over the element.

Using equations (3.50) and (3.51), the discrete local equilibrium equation (3.47a) is finally equivalent to:

[n m ] T τ .n = t n t m (3.52)
We recall that t n and t m are respectively the normal and tangential components of the driving traction that appears in the cohesive law. This cohesive law is expressed in incremental form through the following scalar equations:

ṫn = Can (n,n) αn ṫm = Can (m,m) αm (3.53a) (3.53b)
with Can = Can

(n,n) 0 0 Can (m,m)
.

The tangent modulus of the cohesive law is here supposed to be diagonal, which implies that mode I and II do not interact with each other as it may be possible for some special cohesive laws (see [START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures. Part I : theoretical formulation and numerical implementation[END_REF]).

Using the discretization specified in section 3.1.2, the system to be solved may finally be sumed up as follows:

Find (d, α) such that:      N el A e=1 [f int,e -f ext,e ] = 0 h e = 0, ∀e ∈ [1, N el ] (3.54a) (3.54b)
with

                           f int,e = Ω e B T x τ dΩ f ext,e = Ω e N T bJ dΩ + ∂tΩ e N T hJ dΓ h e = Ω e \Γs ḡT v τ | Ω e \Γs dΩ + Γs ḡv T t Γs dΓ = -Γs [n m ] T n ,T τ | Ω e \Γs + Γs t n t m (3.55a) (3.55b) (3.55c) (3.55d)
where d are the nodal displacements, α are the displacement jumps projected along the different modes described in section 3.1.2, and N el is the number of elements. N = [N 1 N 2 N 3 ] is the matrix of shape functions for a standard T3 element,

B x = [B x,1 B x,2 B x,3 ] = L x N with L x the matrix form of the ∇ s
x operator referring to the shape functions derivatives in the spatial configuration. J is the jacobian of the transformation, such that J = det( F). Remark 3.1.2. It has to be noticed that the stresses τ are a regular distribution over Ω e , then we can write the following relation:

f int,e = Ω e B T x τ dΩ = Ω e \Γs B T x τ dΩ (3.56)
Remark 3.1.3. In the following part that deals with the linearization of system (3.54), the external forces f ext,e are not taken into account since the examples presented in section 3.4 do not account for external forces. However, if these external forces were to be accounted for, it has to be noticed that their linearization should not be omitted since the jacobian of the transformation J appears in their formulation in (3.55b).

Linearization of the discretized problem

The non-linear system (3.54) needs to be linearized with respect to (d, α) in order to be solved in a Newton procedure. A detailed presentation of the full linearization of (3.54) is proposed in appendix C. The final linearized system is written as follows:

     N el A e=1 [δf int,e ] = N el A e=1 [f ext,e -f int,e ] h e + δh e = 0, ∀e ∈ [1, N el ] (3.57a) (3.57b)
where

δf int,e = Ω e B ns,T x H T c an T r + τ B ns x dΩδd
-

Ω e B ns,T x H T c an T r + τ a∈Ω + B ns a,x .[n m ]dΩδα (3.58)
and

δh e = -[n m ] T n ,ns,T τ + τ a=3 a=1 B ns a,x + c an T r a=3 a=1 B a,x δd + [n m ] T n ,ns,T τ + τ a∈Ω + B ns a,x + c an T r a∈Ω + B a,x .[n m ] + Can (n,n) 0 0 Can (m,m) δα (3.59)

Resolution of the discretized problem

The numerical solution of the linearized problem is obtained by introducing a pseudotime "t". Let us consider a global iteration i of a time step n + 1. We denote:

d(t n+1 ) (i) = d (i) n+1 α(t n+1 ) (i) = α (i) n+1 (3.60) and ∆d (i) n+1 = d (i+1) n+1 -d (i) n+1 ∆α (i) n+1 = α (i+1) n+1 -α (i) n+1 (3.61)
With this "time" discretization expressions, the linearized system (3.57) to be solved may be written:

             N el A e=1 ∂f int,e ∂d (i) n+1 ∆d (i) n+1 + ∂f int,e ∂α (i) n+1 ∆α (i) n+1 = N el A e=1 f ext,e(i) n+1 -f int,e(i) n+1 h e(i) n+1 + ∂f int,e ∂d (i) n+1 ∆d (i) n+1 + ∂f int,e ∂α (i) n+1 ∆α (i) n+1 = 0 (3.62)
which can be rewritten in matrix form:

       N el A e=1 K e(i) f d,n+1 ∆d (i) n+1 + K e(i) f α,n+1 ∆α (i) n+1 = N el A e=1 f ext,e(i) n+1 -f int,e(i) n+1 h e(i) n+1 + K e(i) hd,n+1 ∆d (i) n+1 + K e(i) hα,n+1 ∆α (i) n+1 = 0 (3.63a) (3.63b)
where the matrices involved in this system rely on the linearization detailed above. They are defined as:

K e(i) f d,n+1 = Ω e \Γs B T x c an T r B x dΩ + Ω e \Γs B ns,T x τ B ns x dΩ K e(i) f α,n+1 = Ω e \Γs B T x c an T r ḡr dΩ + Ω e \Γs B ns,T x τ ḡns r dΩ K e(i) hd,n+1 = ḡv ns,T c an T r B x + τ + τ B ns x K e(i)
hα,n+1 = ḡv ns,T c an T r ḡr + τ + τ ḡns Remark 3.1.4. It has to be noted that, contrary to the small strain formulation in section 2.1.2, and since we choose to formulate the problem in the current configuration, all the quantities involved in the formulation, and particularly the different matrices containing the derivatives of the shape functions as well as the spatial vectors (n , m ) and (n , m ), depend on the current pseudo-time step t n+1 and the global iteration i. This is why their references are omitted in the definition of the matrices in (3.64).

Remark 3.1.5. It also has to be noted that, contrary to the small strain case, the matrices in (3.64) contain a "geometric" part that is not related to the material model, but to the non-linear transformation from the undeformed to the deformed configuration. All these matrices can be put in the following general form:

K e(i) * ,n+1 = K mat,e(i) * ,n+1 + K geo,e(i) * ,n+1 (3.65)
where K

mat,e(i) * ,n+1

contains the material tangent modulus among other quantities, and K geo,e(i) * ,n+1 contains the "geometric" terms.

The resolution of this sytem is carried out using the same "operator split" method as in the infinitesimal case. Thus, the reader is invited to go to section 2.1.2 for further details concerning this matter.

The final problem can then be put in the general following form:

N el A e=1 Ke(i) n+1 ∆d (i) n+1 = N el A e=1 f ext,e (i) n+1 -f int,e(i) n+1 (3.66)
with

Ke(i) n+1 = K e(i) f d,n+1 -K e(i) f α,n+1 K e(i) hα,n+1 -1 K e(i) hd,n+1 (3.67)
System (3.66) is equivalent to the standard non-linear Finite Element problem, the only difference is the formulation of the element tangent stiffness matrix Ke(i)

n+1
given in (3.67). This matrix is not symmetric which involves an appropriate resolution of system (3.66). The fact that the global Finite Element procedure remains unchanged highlights the main advantage of this method. Only a new element needs to be implemented, the rest of the FE resolution architecture remains unchanged, while the introduction of new degrees of freedom, as for X-FEM for example, involves a specific treatment at the global level.

The full enhanced Finite Element resolution procedure is sumed up in the following box:

Data: d n , α n Result: d n+1 , α n+1 i = 0; while A N el e=1 Ke(i) n+1 ∆d (i) n+1 -A N el e=1 f ext,e n+1 -f int,e n+1 > tol do for e = 1 → N el do Compute deformation gradient F(i) n+1 = I h + B ns ∆d (i) n+1
Compute Kirchhoff stress and tangent modulus τ

(i) n+1 = f ( F(i) n+1 , q n ) c an,(i) T rn+1 if localization not yet detected then Ke(i) n+1 = K e(i) f d,n+1 else j = 0 α (i,j) n+1 = α n while h e(i,j) n+1 > tol do Fh,(i,j) n+1 = I h + B ns ∆d (i) n+1 + Ḡr ∆α (i,j) n+1
begin Compute material models for the bulk and disc. surface τ

(i,j) n+1 = f ( F(i,j) n+1 , q n ) c an,(i,j) T rn+1 ∆t (i,j) Γs,(n,m),n+1 = Can (i,j) (n,m),n+1 ∆α (i,j) n+1 Can (i,j) n+1 end [n ,(i,j) n+1 m ,(i,j) n+1 ]
= Fh,(i,j)

n+1 [N M] [n ,(i,j) n+1 m ,(i,j) n+1 ] = Fh,(i,j) -T n+1 [N M] h e,(i,j) n+1 = -Γs [n ,(i,j) n+1 m ,(i,j) n+1 ] T n ,T,(i,j) n+1 τ (i,j) n+1 + Γs t (i,j) Γs,n+1 α (i,j+1) n+1 = α (i,j) n+1 -K e(i,j) hα,n+1
-1 h e,(i,j) n+1 j = j + 1 end Compute updated deformation gradient, Kirchhoff stress and tangent modulus

F(i) n+1 = I h + B ns ∆d (i) n+1 + Ḡr ∆α (i) n+1 , τ (i)
n+1 = f ( F(i) n+1 , q n ) , c an,(i) T rn+1
Compute enhanced tangent stiffness matrix Ke,(i)

n+1 = K e,(i) f d,n+1 -K e,(i) f α,n+1 K e, (i) hα,n+1 -1 K e,(i) hd,n+1 end B (i)
x n+1 = f( Fh,(i) n+1 , B) f int,e(i) n+1 = Ω e B (i),T x n+1 τ (i) n+1 dΩ f ext,e(i) n+1 = Ω e N T b (i) n+1 J dΩ + ∂tΩ e N T h (i) n+1 J dΓ end Solve A N el e=1 Ke,(i) n+1 ∆d (i+1) n+1 = A N el e=1 f ext,e(i) n+1 -f int,e,(i) n+1 to find d (i+1) n+1 i = i + 1 end
Algorithm 2: Full resolution procedure where q n is the set of internal variables for the bulk constitutive model computed at time t n .

Material model formulation and implementation

In this section, the material model for the bulk material is formulated and its numerical implementation is discussed. Furthermore, a criterion different from the small deformation case is formulated and discussed.

Bulk material model: Large strain von Mises plasticity model

We recall that the bulk material has to be able to model the bulk dissipative mechanisms that occur during the pre-peak phase of the material response. This means that, at least the plastic dissipation has to be taken into account in our material model. In order to circumvent the convergence issues in the implicit non-linear global problem resolution induced by Lemaitre's coupled damage-plasticity model, we choose to use an isotropic von Mises plasticity model for the bulk material, combined with a Rice & Tracey void growth model and its associated coalescence criterion, that we will use as the criterion that triggers the introduction of a discontinuity surface. However, combining a coupled damage-plasticity model such as Lemaitre's in a finite strain scenario is in perspective.

Since in this part of the work a large strain hypothesis is considered, the von Mises plasticity model is not straightforward, and a proper adaptation of the small strain case has to be made. The following section deals with the formulation of the von Mises J2-plasticity model in the large strain scenario.

Theoretical formulation

Contrary to the small strain hypothesis in which the additive decomposition of the total strain ε into an elastic part and a plastic part is assumed, in a finite strain hypothesis, we use a multiplicative decomposition of the total deformation gradient F:

F = F e F p (3.68)
where we recall that the deformation gradient defines the transformation between the undeformed and the deformed configurations:

F = ∂x ∂X (3.69)
where X ∈ Ω X refers to the reference configuration (or undeformed) and x ∈ ω x = ϕ(Ω X ) refers to the deformed (or spatial) configuration.

See figure 3.6 in which the motion decomposition is schematically represented. Such a decomposition has been considered by many authors for the past forty years (see e.g. Lee and Liu [START_REF] Lee | Finite strain elastic-plastic theory particularly for plane wave analysis[END_REF], Lee [START_REF] Lee | Elastic-plastic deformations at finite strains[END_REF], Kroner and Teodosiu [START_REF] Kroner | Lattice defect approach to plasticity and viscoplasticity, Problems of Plasticity[END_REF], Mandel [START_REF] Mandel | Contribution théorique a l'etude de l'écrouissage et des lois de l'écoulement Plastique[END_REF][START_REF] Mandel | Thermodynamics and plasticity, Foundations of Continuum Thermodynamics[END_REF], Kratochvil [START_REF] Kratochvil | On a finite strain theory of elastic-inelastic materials[END_REF], Sidoroff [START_REF] Sidoroff | Un modèle viscoélastique non linéaire avec configuration intermediaire[END_REF], Nemat-Nasser [START_REF] Nemat-Nasser | On finite deformation elasto-plasticity[END_REF], ...etc...). This decompostion relies on the principle of an intermediate local stress-free configuration which is now known to lead to a simple class of integration algorithms, and in particular to a straightforward extension of the return-mapping algorithm of the infinitesimal J2-plasticity theory, which is of great interest in our case. This principle relies on micromechanical considerations with regards to the plastic flow occurring in a single-crystal, studied by G.I. Taylor and other authors (see [START_REF] Taylor | Bakerian lecture: the distortion of an aluminum crystal during a tensile test[END_REF][START_REF] Taylor | The plastic extension and fracture of aluminum crystals[END_REF]). They defined the plastic part of the deformation gradient as the dislocation motion of a crystal subjected to a single slip-system, while the elastic part is the deformation caused by stretching and rotation of the crystal lattice. This principle has been extended to crystal subjected to several slip systems, for which a review of articles has been proposed by Asaro [START_REF] Asaro | Micromechanics of crystals and polycrystals[END_REF].

The principle of the local intermediate stress-free configuration implies that we can interpret the tensor F e-1 as the deformation that releases the stresses from the deformed configuration.

As it was presented in section 3.1.1, we chose to write the governing equilibrium equations in the spatial (or Eulerian) configuration since this configuration deals with the Kirchhoff stress tensor τ which remains symmetric as the Cauchy stress tensor σ, which is computationally more convenient. Thus we need to define all the strain and stress tensors in the spatial configuration. Associated to this current configuration, the left Cauchy-Green strain tensor is defined as:

b = FF T (3.70)
Another useful strain measure is the right Cauchy-Green strain tensor (related to the reference configuration) defined as:

C = F T F (3.71)
and its plastic part:

C p = F pT F p (3.72)
which is useful to define the elastic left Cauchy-Green tensor:

b e = F e F eT = FF p-1 F p-T F T = F F pT F p -1 F T (3.73)
which, in view of 3.72 can be rewritten:

b e = FC p-1 F T (3.74)

Model formulation

Assuming an isotropic stress response and an incompressible plastic flow, we can write:

J p = detF p = detC p = 1 (3.75)
which leads to:

J = detF = detF e detF p = detF e = J e (3.76)
Considering these assumptions, we proceed to the model formulation as follows. This model was first proposed by Simo (see [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition: Part I. Continuum formulation[END_REF] for the continuum formulation and [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition[END_REF] for its integration algorithm). Assuming the notion of a stress-free intermediate configuration, the free energy function is of the form: Ψ = U (J e ) + Ψ( be ) (3.77) where be = J e-2/3 b e is the volume preserving elastic left Cauchy-Green tensor and U and Ψ are the volumetric and deviatoric part of the total free energy Ψ, respectively. They take the following form:

U (J e ) = 1
2 κ 1 2 (J e2 -1) -ln(J e ) Ψ( be ) = 1 2 µ tr[ be ] -3

(3.78)
where κ and µ are the bulk modulus and shear modulus, respectively. This uncoupled stored energy function leads to an uncoupled volumetric-deviatoric stressstrain relationship, defined as follows:

τ = J e p1 + s (3.79)
where 1 is the second-order identity tensor, and with:

p = U (J e ) = κ 2 (J e2 -1)/J e s = dev[τ ] = µ dev[ be ] (3.80)
Let us now define the Mises-Huber yield function, which is the finite strain equivalent of the von Mises infinitesimal yield function:

φ(τ , r) = ||dev[τ ]|| -2 3 [σ y + q(r)] (3.81)
where q(r) is the non-linear hardening function and r is the standard hardening parameter.

The associative-flow rule is uniquely defined by the principle of maximum dissipation, and this flow rule can be defined in the spatial description as:

L v b e = -2 3 γtr[b e ]n (3.82)
where L v is the Lie derivative, γ is the plastic multiplier, and where as in the infinitesimal theory, n is defined as:

n = s ||s|| (3.83)
Since in the integration algorithm, detailed in a further paragraph, the flow rule is defined in terms of the plastic strain tensor Cp-1 n+1 at time step n + 1, the flow rule can also be defined in the reference configuration as:

∂ ∂t Cp-1 = -2 3 γ tr[b e ]F -1 nF -T (3.84)
where Cp-1 is the inverse of the volume preserving plastic right Cauchy-Green tensor.

The evolution of the hardening parameter is defined as:

ṙ = 2 3 γ (3.85)
and γ is the plastic multiplier that satisfies the standard Kuhn-Tucker loading/unloading conditions: γ 0, φ(τ , r) 0, γφ(τ , r) = 0 (3.86) with the consistency condition:

γ φ(τ , r) = 0 (3.87)

Integration algorithm

We now deal with the integration algorithm of the aforementioned constitutive model. It was developed by Simo [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition[END_REF][START_REF] Simo | Computational inelasticity[END_REF]. Using a backward Euler difference scheme together with the plastic flow rule (3.84) and the evolution of the hardening variable (3.85), we obtain the following discrete evolution equations:

Cp-1 n+1 - Cp-1 n = -2 3 ∆γ C p-1 n+1 : C n+1 F -1 n+1 n n+1 F -T n+1 r n+1 -r n = 2 3 ∆γ (3.88)
After some algebraic manipulations, the spatial discrete evolution equations are written:

             be n+1 = fn+1 be n fT n+1 -2 3 ∆γ tr[ be n+1 ]n n+1 n n+1 = s n+1 ||s n+1 || s n+1 = dev[τ ] r n+1 = r n + 2 3 ∆γ (3.89)
where fn+1 is the volume preserving part of the deformation gradient between the deformed configuration at time step t n and the deformed configuration at time step t n+1 . It is expressed as:

fn+1 = Fn+1 F-1 n = (J n+1 /J n ) -1/3 f n+1 = (J n+1 /J n ) -1/3 F n+1 F -1 n (3.90)
where f n+1 = F n+1 F -1 n . With these developments at hand, we can now define the return-mapping algorithm, and as in the infinitesimal case, it is a strain driven algorithm. It is based on the elastic predictor and plastic corrector standard procedure, which is summed up in the following paragraph.

Let us consider a trial state in which plastic flow does not occur at time step t n+1 . We have then:

C p-1,trial n+1 = C p-1 n r trial n+1 = r n (3.91)
Then, we can express the trial elastic state of the different variables:

             τ trial n+1 = p n+1 J n+1 1 + s trial n+1 s trial n+1 = µ dev[ be,trial n+1 ] p trial n+1 = U (J n+1 ) be,trial n+1 = fn+1 be n fT n+1 r trial n+1 = r n (3.92)
And the trial yield function is expressed as:

φ trial n+1 = φ(τ trial n+1 , r trial n+1 ) = ||dev[τ trial n+1 ]|| -2 3 [σ y + q(r trial n+1 )] (3.93)
Then, if φ trial n+1 0 the Kuhn-Tucker conditions in (3.86) are satisfied with ∆γ = 0, and this trial elastic state is the solution at time t n+1 . Thus, we can define the updated variables as (•) n+1 = (•) trial n+1 where (•) trial n+1 are defined in equation (3.92).

However, if φ trial n+1 0, then the Kuhn-Tucker conditions are not satisfied and this trial elastic state is not admissible, so that τ n+1 = τ trial n+1 . Time step t n+1 is a plastic step with ∆γ = 0. The radial return procedure corresponding to a plastic step is detailed in the following paragraph.

By taking the trace of be n+1 in expression (3.89), we can write:

tr[ be n+1 ] = tr[ be n+1 ] (3.94)
which, using the expressions of the elastic response in (3.92), leads to:

s n+1 = s trial n+1 -2 3 µ∆γ tr[ be,trial n+1 ]n n+1 (3.95)
After some manipulations in (3.95), we can write:

n trial n+1 = s trial n+1 ||s trial n+1 || = n n+1 (3.96)
which implies the following equation:

||s n+1 || + 2μ∆γ = ||s trial n+1 || (3.97)
where

μ = 1 3 µtr[ be,trial n+1 ] (3.98)
which will be useful to compute the updated yield function φ n+1 . Now we need to compute the plastic multiplier ∆γ, which is done, as in the standard infinitesimal case, by iteratively solving φ(τ n+1 , r n+1 ) = 0 with a Newton scheme. It cannot be directly computed since the hardening function q(r) is non linear in our case. Substituting equation (3.97) in (3.81), we can express the following equation to be solved in order to obtain ∆γ:

φ(∆γ) = ||s trial n+1 || -2 3 σ y - 2 3 q(r n + 2 3 ∆γ) + 2μ∆γ = 0 (3.99)
Then ∆γ is computed using the standard Newton scheme as follows:

∆γ (k+1) = ∆γ (k+1) - φ(∆γ (k) ) φ (∆γ (k) ) (3.100)
which is solved until φ(∆γ (k+1) ) 0.

In equation (3.100), φ (∆γ (k) ) is expressed:

φ (∆γ (k) ) = -2μ 1 + q (r n + 2 3 ∆γ (k) ) 3μ (3.101)
Remark 3.2.1. In the framework of finite strain J2-plasticiy theory, other numerical procedures exist. For example, Simo [START_REF] Simo | Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory[END_REF] developed an elegant procedure that is very similar to the one presented above, the difference being that the procedure relies on the principal stretches defined by the spectral decomposition of the left Cauchy-Green tensor b, which allows to obtain the logarithmic strain ε, the problem then leads to a similar formulation of the infinitesimal problem, which is solved in the principal basis. A straightforward change of basis to the original global basis leads to a similar formulation of the updated Kirchhoff stress tensor, and updated left Cauchy-Green strain tensor. The linearization leads to the formulation of the elastic-plastic modulus in a similar manner as in the infinitesimal case.

We also need to compute the spatial consistent tangent elastic-plastic modulus. For details of the linearization procedure leading to the following formulation of the elasto-plastic modulus, please refer to [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition[END_REF][START_REF] Simo | Computational inelasticity[END_REF].

For an elastic step, the elasticity tensor is expressed as:

c = (JU ) J1 ⊗ 1 -2JU I + c (3.102)
where:

c = 2μ[I -1 3 1 ⊗ 1] -2 3 ||s||[n ⊗ 1 + 1 ⊗ n] (3.103)
and where we recall that:

   s = µdev[ be ] n = s/||s|| μ = 1 3 µtr[ be ] (3.104)
For a plastic step, after the exact linearization proposed by Simo [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition[END_REF][START_REF] Simo | Computational inelasticity[END_REF], we obtain the following expression for the consistent elastic-plastic tangent modulus:

c ep n+1 = c trial n+1 -β 1 ctrial n+1 -2μβ 3 n ⊗ n -2μβ 4 n ⊗ dev[n 2 ] s (3.105)
where [•] s refers to the symmetric part of [•], and where:

                             β 0 = 1 + q (r n+1 ) 3μ
β 1 = 2μ∆γ ||s trial n+1 || β 2 = 2 3 1 - 1 β 0 ||s trial n+1 || μ ∆γ β 3 = 1 β 0 -β 1 + β 2 β 4 = 1 β 0 -β 1 ||s trial n+1 || μ (3.106)
Remark 3.2.2. This consistent elasto-plastic tangent modulus corresponds to the spatial modulus defined as c an T r that appears in the definition of the element tangent stiffness defined in (3.64). We recall that it appears in the rate form of the constitutive equation:

τ T r = L v (τ ) = c an T r : ε = c ep : ε (3.107)
where τ T r is the Trusdell rate of the Kirchhoff stress.

It is actually computed by first "pulling back" the stress tensor s in the reference configuration using the following relation: S = F -1 sF -T , then derivating this material stress tensor with respect to C and finally "pushing forward" this derivative in the spatial configuration with the relation:

c ep ijkl = 2F iI F jJ F kK F lL ∂S IJ ∂C KL (3.108)
where Einstein's summation convention has been adopted. Again, please refer to [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition[END_REF][START_REF] Simo | Computational inelasticity[END_REF] for details regarding the linearization procedure.

Since in the present finite strain formulation the discrete cohesive law has the same formulation as in the small strain case, the reader is invited to refer to section 2.2.2 for details concerning this matter.

Instability issue

There is still a numerical difficulty to circumvent in the presence of strong discontinuities in a non-linear finite element framework. Indeed, when a surface of discontinuity of displacement is introduced, because of the non-smooth transition between the rigid phase and the plastic/damageable phase of the cohesive law, a loss of convergence may occur in some cases during the Newton-Raphson non-linear solving procedure at global level. This is in fact due to the presence of a bifurcation point that requires a special treatment to circumvent this issue. At this bifurcation point, two branches are solutions of the global equilibrium equation. Thus numerically, the solution oscillates between these two solutions until the maximum number of iterations is reached at the time step of interest.

To overcome this problem, a possible way is to use a small perturbation to smooth the different equilibrium paths occurring when crossing a bifurcation point. It remains to make sure that the solution lies on the right equilibrium path, which corresponds to the path with minimal energy. This is done considering that the unique stable path remains continuous, so continuation should be ensured with the use of small enough time steps.

The use of perturbation methods implies that the solution loses its physical sense, the "amount of loss of physical sense" being proportional to the amount of perturbation. Thus, the perturbation should remain small enough to keep the physical sense of the solution.

In the following subsection, the viscosity perturbation method introduced by Oliver et al. [START_REF] Oliver | Theoretical and computational issues in modelling material failure in strong discontinuity scenarios[END_REF] is detailed.

Viscosity regularization

The method introduces an artificial viscous behavior of the discontinuity surface. An additional term appears in the weak form of the traction continuity condition (3.27b). It is in this context expressed as follows:

∀(w, η), Ω\Γs ηs x : τ dΩ + Γs w.[τ | Γs .n + γ u] dΓ = 0 (3.109)
The viscosity parameter γ drives the amount of perturbation of the initial solution.

The constitutive law on the discontinuity surface is expressed in rate form as:

ṫ * Γs = Can u + γ ü (3.110)
where ṫ * Γs is the total traction on the discontinuity which includes the viscosity contribution.

Projecting this equation on basis (n , m ), we obtain the following expression for the time derivatives of the components of the driving traction ṫ * Γs on basis (n , m ): where we have used the decomposition: ū = α n n + α m m , and its time derivation by a "pull-back/push-forward" operation, as it is detailed in appendix C, and particularly in (C. [START_REF] Benallal | Crack initiation and propagation analyses taking into account initial strain hardening and damage fields[END_REF]). Using the time discretization proposed in section 3.1.2, this constitutive equation is written in incremental form as:

∆t * ,(i) (n,m),n+1 = Can,(i) (n,m),n+1 ∆α (i) n+1 + γ∆ α(i) n+1 (3.112)
where ∆t * ,(i)

(n,m),n+1 = ∆t (i) n,n+1 ∆t (i) m,n+1 , C(i) (n,m),n+1 = Can,(i) (n,n),n+1 0 0 Can,(i) (m,m),n+1
is the consistent tangent moduli of the discontinuity surface constitutive law ( Can,(i) (n,n),n+1 for mode I and Can,(i) (m,m),n+1 for mode II), and where ƥ

(i)
n+1 is defined as:

∆• (i) n+1 = • (i+1) n+1 -• (i) n+1 (3.113)
Using an Euler implicit scheme to approximate the discontinuity jump rate α(i) n+1 , we have:

α(i) n+1 α (i) n+1 -α n t n+1 -t n (3.114)
which gives: (3.117)

∆ α(i) n+1 = α(i+1) n+1 - α(i) n+1 = α (i+1) n+1 -α n t n+1 -t n - α (i) n+1 -α n t n+1 -t n = ∆α (i) n+1 ∆t n+1 (
The updated solution procedure in box 3.1.2 is then the following: Oliver and Huespe [START_REF] Oliver | Theoretical and computational issues in modelling material failure in strong discontinuity scenarios[END_REF] highlighted the fact that associating an artificial viscosity behavior to the discontinuity surface induces a critical time step above which following the stable continuous equilibrium path is no longer ensured. Indeed, since in the case of a quasi-static problem the resolution relies on a pseudo-time discretization, for a given time-step with a given perturbation, nothing guarantees that the response will lie on the stable equilibrium branch, and the solution might fall on the unstable branch. The method detailed above only ensures the continuity of the stable equilibrium path, but does not ensure its uniqueness. Thus, an appropriate procedure should be used to ensure that the pseudo-time discretization allows the solution to remain on the stable equilibrium path. In other words, this procedure should control the pseudo-time step discretization with regards to the amount of perturbation introduced.

Find (d, α) such that:      N el A e=1 [f int,e -f ext,e ] = 0 h e = 0, ∀e ∈ [1, N el ] (3.118a) (3.118b) with                            f int,e = Ω e
The following deals with the critical time step evaluation, which directly allows an automatization of the time step to ensure the algorithmic uniqueness of the solution.

We recall that the fully linearized system (3.63) is written :

    N el A e=1 0 0 χ Ω e,loc        K e(i) f d,n+1 K e(i) f α,n+1 K e(i) hd,n+1 Ke(i) hα,n+1       ∆d (i) n+1 ∆α (i) n+1    =    f ext,e(i) n+1 -f int,e(i) n+1 0    (3.120)
where

    N el A e=1 0 0 χ Ω e,loc   
 is an interpretation of the fact that the first set of equations is assembled over the whole structure, while the second set of equations is only written for the localized elements (χ Ω e,loc is the characteristic function of the domain Ω e,loc ). This system can be rewritten:

[ K (i) n+1 ]∆β (i) n+1 = r (i) n+1 (3.121)
with

β (i) n+1 =    d (i) n+1 α (i) n+1    (3.122)
and

r (i) n+1 =    f ext,e(i) n+1 -f int,e(i) n+1 0    (3.123)
If the solution is not unique at a given time step t n+1 , we have:

∆β (1),(i) n+1 -∆β (2),(i) n+1 = 0 (3.124)
and

[ K (i) n+1 ](∆β (1),(i) n+1 -∆β (2),(i) n+1 ) = 0 (3.125)
which implies the singularity of matrix [ K

(i) n+1 ].
Thus, the positive definiteness of matrix [ K (i) n+1 ] will ensure the uniqueness of the solution. Ensuring this positive definiteness will then allow us to define a relation between the viscosity parameter and the critical time step during the discontinuity evolution.

If we consider that the original stiffness matrix [K e(i) f d,n+1 ] defined in (3.64) is positive definite, which is consistent with the fact that we always introduce a discontinuity while the global tangent stiffness remains positive (see remark 3.3.2 below), then a necessary condition for the positive definiteness of matrix [

K (i) n+1 ] is the positive definiteness of matrix [ Ke(i) hα,n+1
]. Oliver and Huespe [START_REF] Oliver | Theoretical and computational issues in modelling material failure in strong discontinuity scenarios[END_REF] proposed that a reasonable approximation of this condition is to state that the symmetric part of this matrix remains positive, which is expressed as:

[K e(i) hα,n+1 + 1 γ ∆t n+1 ] s → positive ∀e ∈ [1, N el ] (3.126)
which is equivalently expressed as

λ e(i) min,n+1 + γ ∆t n+1 > 0 ∀e ∈ [1, N el ] (3.127)
with λ e min being the minimum eigenvalue of elementary matrix [K e(i) hα,n+1 ]. Then, the critical time increment may be expressed as: In a practical point of view, the critical time step should be computed for every iteration of a given time step and for each localized element. Finding the minimum of these values is simply done, at element level, by comparison between the current element's critical time step and the critical time step computed for the last localized element. At the end of the current converged time step, the new value of the time incement is set according to the computed critical time increment.

∆t crit n+1 = min
An alternative to ensure uniqueness is to fix the time increment and to adapt the value of the perturbation, while still ensuring equation (3.126). However, considering that the value of the perturbation evolves independently from the user, nothing guarantees that the perturbation is small enough to keep the physical sense of the solution. Both cases are discussed in section 3.3.1 below. Remark 3.3.2. Working in the finite deformation range implies that even with an elastic behavior, the global response can enter softening due to the geometric nonlinearities. As a consequence, the solution may not be unique and singularity of the global tangent stiffness matrix may occur. In a structural analysis point of view, this means that for a simple tensile test, material striction occurs even for an elastic behavior, and the response becomes mesh-dependent.

The purpose of the strong discontinuity method is to preclude this singularity by introducing a surfacic behavior that is in charge of modeling all the softening behavior of the structure. Thus, it should always be introduced before the global response begins to soften.

In other words, the criterion used to introduce the discontinuity should contain both the "material" condition and the "geometric" condition. Contrary to the small strain case, the positive definiteness of the global stiffness tangent matrix doesn't only involve positive definiteness of the acoustic tensor: det(n.C (an) .n) > 0 ∀n (3.129) but it also contains some geometric information since the global tangent stiffness matrix may now be written:

[ K] = [ Kmat ] + [ Kgeo ] (3.130)
Such a unified localization condition is not established here, but we recall that it should be kept in mind when establishing a localization criterion. Authors have worked on the establishement of a unified "geometric" and material localization condition (see Armero [START_REF] Armero | An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids[END_REF] for example). Having considered this point, in the following, we ensure that the discontinuity is always introduced before the global response begins to soften by considering reasonable strain levels (below the "geometric" condition) at the onset of failure. This is generally true in a ductile metal failure analysis in which equivalent fracture strain levels are in general below 0.3.

These considerations mean that it is reasonable to consider that the matrix K e(i) f d,n+1 is always positive definite.

Numerical examples

In this section, examples attesting the performance of the viscous regularization detailed above are shown. As an introducing example, a simple tensile test considering plane strain hypothesis is shown. As the test is completely homogeneous, nothing in the model drives the crack initiation and the discontinuity is then introduced within every element of the structure exactly at the same time. This is a good example to show the loss of uniqueness problem and the numerical issues that it implies at global level, since several discontinuities are activated at the same time. The considered mesh along with the boundary conditions are shown in figure 3.7.

The material is an elasto-plastic material, following the J2-plasticity model detailed in section 3.2, with an exponential hardening law. The criterion used to introduce the discontinuity is a Rice and Tracey criterion, as detailed in section 3. Fracture energy density G c = 1000 MPa.mm Table 3.1: Model parameters for the volumic and surfacic models If no procedure is undertaken to regularize the bifurcation problem, then at the onset of localization convergence is lost, as shown in table 3.2 which exhibits the residual norm and the energy norm for each iteration of the non-converged step. The time increment is the same for every time step of the computation and is set to 0.01.

It is obvious that the solution oscillates between two non-converged solutions, one exhibiting a negative energy norm, which is typical of loss of "positivity".

The same test is performed with the viscous regularization implemented, but without the auto time-stepping procedure. The value of the viscosity is set to 0.1 which is low enough to keep a physical sense to the solution. Results in terms of convergence are shown in table 3. [START_REF] Armero | New finite elements with embedded strong discontinuities in the finite deformation Range[END_REF] We can see that convergence is still not ensured. The time increment is actually too large and the solution falls in the unstable branch of the equilibrium path. It is obvious that a path continuing procedure such as the auto time stepping procedure detailed above is necessary to preclude this branch switching.

In this last test, the auto-time stepping procedure is implemented and results in terms of convergence are shown in table 3 It is shown that convergence is ensured, and the time step dramatically decreases to a value of ∆t = 6.568 × 10 -5 s, while it was initially set to ∆t = 1 × 10 -2 s at the previous time step.

An alternative procedure to ensure uniqueness would be to set the viscosity parameter automatically with regards to the current time increment, while still ensuring equation (3.127). Then the critical viscosity parameter is computed as:

γ crit = min e |λ e min,n+1 | ∆t n+1 (3.131)
The main advantage of this technique is that the time increment remains fixed and thus does not dramatically decrease as in the previous test (see table 3.5). The main drawback is that nothing drives the amount of viscosity with regards to the physical sense of the solution. In order to show the loss of physical sense of the solution, a final test is performed with exactly the same parameters as the previous test, the unique difference being that the viscosity parameter is set automatically using equation (3.131).

The responses for both cases (automatic time increment and automatic viscosity parameter) are plotted in figure 3 As it is shown in table 3.5, the automatically computed viscosity parameter is quite high and thus explains why both curves in figure 3.8 are not superimposed. Even though with the automatic viscosity computation, the solution remains unique and is computed very quickly, it is obvious that for this test the physical sense of this solution is not well preserved, as we can see in figure 3.8 that a non-negligeable residual force remains at the end of the test.

In the following paragraph, a sensitivity analysis is performed with regards to the viscosity parameter for the same test and for the case of automatic time increment computation. Results in terms of global responses are shown, along with a table gathering results in terms of total CPU time, average automatic time increment and number of time steps (see table 3 We can see from table 3.6 that very low values of the viscosity parameter imply dramatically high computational costs, while they allow to keep a realistic physical sense. On the other hand, high values of the viscosity parameter induce a loss of physical sense while implying very low computation costs. From these results and for this test, we conclude that reasonable values for the viscosity parameter are between 0.1 and 1. This range allows the physical sense to be realistic while keeping reasonable computational costs. However, this conclusion might only be true for this test setting, and other numerical examples should be performed to confirm these conclusions.

Numerical examples

In this section, every mentioned numerical example considers the regularization method mentioned above. Mode I, mode II and mixed mode discontinuity opening are also considered depending on the studied case. As it was performed in earlier numerical examples, a simple tensile test is undertaken to show the regularizing capabilities of the method. An element is "weakened" in the lower region of the mesh to simulate the presence of a micro-defect which drives the onset of localization, which is performed using a Rice and Tracey criterion. For this "weakened" element, the Rice & Tracey criterion is chosen to be slightly lower (( R R 0 ) c = 1.29) than for the rest of the mesh ( R R 0 ) c = 1.30). As for the rest of the material models, the same parameters as in section 3.3.1 are used for this test, and the normal to the discontinuity surface is chosen according to the maximum principal stress direction, as for the examples presented in section 3.3.1.

The considered fine mesh contains identical elements with a size ten times smaller than the coarse mesh (the mesh considered for all the previous tests in section 3.3.1). Results in terms of reaction/displacement are shown in figure 3.11 along with the meshes and boundary conditions. The final deformed meshes are plotted on figure 3.12 where the discontinuity opening in mode I is also plotted. It is obvious that the post-peak phase of the response is regularized for this example. The only difference between the responses is a slight shift of the peak due to difference of the size of the "weakened" element used to onset the localization. Indeed, for the coarse mesh, once the first element is localized, half of the section of the specimen is unloading due to the cohesive behavior in the localized element, while only a tenth of the section specimen is unloading at the onset of localization for the fine mesh. This explains the slight difference between the post-peak phases for the two considered meshes.

Simple tensile test exhibiting mixed mode fracture

A similar test was performed, this time considering mode II and mixed mode fracture. An element is also "weakened" at the bottom of the specimen, using the same values for the Rice & Tracey criterion as in the previous example. We recall that the constitutive model used for the bulk material is a von Mises plasticity model.

For this test, the discontinuity surface is introduced in such a way that its normal is normal to the maximum shear direction, that is 45 • with respect to the maximum principal direction. This direction corresponds approximately to the direction of localization given by the loss of ellipticity mathematical condition, in the case of a von Mises plasticity constitutive model for a simple tensile test in plane strain hypothesis. There are actually two directions that give an even maximum shear, we choose to arbitrarily impose the direction of the discontinuity surface in such a way that the scalar product of m (the tangent vector to the discontinuity surface) with the maximum principal direction n I remains negative. The total fracture energy density is set to G c = 2000 M P a.mm. This fracture energy density actually corresponds to the sum of the fracture energy density in mode I and in mode II, corresponding respectively to the energy needed to fully open the crack and to its full sliding, until there is no residual traction left. As a first approximation, this distribution of fracture energy density between both modes is chosen to be equal, such that G cI = G cII = 0.5G c .

A first set of results presented shows mode I opening and mode II sliding at the end of the test for three different mesh sizes in figures 3.13(a) and 3.13(b), respectively. From these figures, we can see that, on the one hand the crack propagates at 45 • from the weakened element until the other side of the specimen, and on the other hand, mode I opening remains infinitesimal with regards to mode II sliding, which confirms the fact this test exhibits a mode II fracture opening mode, which is attested in many experiments documented in the literature. We can also see that for the intermediate and fine meshes, the discontinuity surface appears not only in the crack path, but also in the neighborhood of the crack path. In this area, the discontinuity surface only exhibits residual mode I opening which remains insignificant with regards to the opening in the actual crack path. This is due to the redistribution of the stress at each time step, indeed, the adjacent element to the crack tip tends to be naturally more loaded than the other neighboring elements. It is also to be noted that the "activation" of these neighboring elements do not influence the convergence rate of the solution.

In figure 3.14, the responses corresponding to these three tests are presented, in which it is obvious that the solution remains size-independent.

We now focus on the influence of the distribution of the total fracture energy density between mode I and mode II. In the previous test, we recall that we chose For all these tests, only the coarse mesh used in the previous example is considered, the purpose being to show which mode is predominant when the distribution of fracture energy densities changes.

Mode I opening and mode II sliding at the end of the tests are presented in figures 3.15 and 3.16.

The responses for these different tests are presented in figure 3.17.

From this figure, we can see that the mixed mode case corresponds to a more "brittle" response of the specimen, exhibiting less dissipation. This was predictable since, for theses cases, the fracture energy density corresponding to the "activated" single mode is lower than for the cases "mode I alone" or "mode II alone". This test actually never exhibits both modes at the same time, one always become predominant, which is exhibited in figures 3.15 and 3.16. From these figures, we can Figure 3.17: Reaction/displacement curves different fracture energy density distributions conclude that the dominant mode is the mode with the lowest fracture energy density, which again was a predictable result since a minimum fracture energy density for a considered mode means a quicker loading of the cohesive surface for this mode.

Compact Tension test (CT)

To attest the regularizing capabilities of the aforementioned method, a compact tension test (CT) is performed for different meshes. This test is usually performed to compute the stress intensity factor in mode I (K 1c ) at the crack tip, since an initial crack is present in the specimen. However, in our case we only use this test to validate the capabilities to predict the crack propagation path, and to ensure that for different mesh sizes, the computed maximum loads and displacements at maximum load are close enough for different mesh sizes. Half of the specimen is modelled since the test is symmetrical. The geometry and boundary conditions are shown in figure 3.18. As it was mentioned, a crack already exists which induces a singular stress state at the crack tip. This means that a fine enough mesh has to be considered to catch the gradients in the vicinity of the crack tip in order to obtain a converged result even in the pre-peak phase.

A vertical displacement is prescribed at the top of the pin, which is considered to be a rigid body. The material model and parameters are the same as for the previous test, the Rice and Tracey criterion is set at R R 0 c = 1.2. The different tested meshes are presented in figure 3.19. For the finest mesh, the elements in the region of interest (along the crack path) are approximately four times smaller than for the coarsest mesh. Results in terms of reaction and displacement at the top of the pin are presented in figure 3.20. The deformed mesh ("coarse" mesh) is also presented in figure 3.21. Results show good convergence, especially in terms of maximum load and dis-placement at maximum load, as well as for the description of the softening phase, even though this test is very mesh-dependent since a crack already exists in the specimen at the beginning of the test. Furthermore, it is well-known that T3 elements converge slower than Q4 elements for example, this is another reason to model this test with small enough T3 elements.

Conclusion

In this chapter, an extension of the strong discontinuity formulation to its finite range was proposed, based on the work of Armero et al. (see [START_REF] Armero | An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids[END_REF]). The theoretical formulation of the method was proposed considering plasticity evolution, and its finite element implementation was described, using a current configuration formulation, that is using Kirchhoff stress tensors. This choice was made to ensure that the constitutive laws usually obtained considering small strain hypothesis can be easily transferred to the finite range, without major additional developments. A Rice and Tracey void coalescence criterion was presented and justified for the extension of SDA in the finite deformation framework.

In a further section, we focused on the well-known issue in the presence of bifurcation points exhibited by the introduction of a strong discontinuity of the displacement field, and its influence on the convergence rates during the global Newton-Raphson resolution procedure. A method using the introduction of an artificial viscous behaviour within the cohesive law was described. We showed that adding this artefact in the formulation exhibited the presence of critical pseudo-time steps. To circumvent this issue, we proposed a method to compute this critical time step, and more particularly to automatically compute the pseudo-time increment in order to keep the physical sense of the solution, based on the work of Oliver and Huespe (see [START_REF] Oliver | Theoretical and computational issues in modelling material failure in strong discontinuity scenarios[END_REF]). Finally, results attesting the performance of the proposed method were shown for a simple tensile test and for a compact tension test, using a Rice and Tracey criterion to drive the appearance of a dissipative surfacic behaviour, thus modeling the failure of the specimen.

Moreover, a Compact Tension test was presented in terms of the capabilities of the formulation to predict the ultimate load and the displacement at ultimate load. Finally, a simple tension test was performed in mode II and mixed mode, in order to show the versatility of the method and the influence of the crack energy density distribution between mode I and mode II on the final result and the crack opening and sliding.

An obvious drawback of the T3 element presented in this chapter is its tendency to exhibit locking when subjected to quasi-incompressible strain, which is the presently the case since we deal with the J2-plasticity material model for the bulk. The work presented in this chapter should be seen as a first step towards a more complex formulation that aims to overcome the locking effect as it was proposed with the enhancement of the standard Q4 element with a "B-bar" formulation. Indeed, a perspective of this work on the T3 element would be its adaptation to a Q4 formulation considering a B-bar enhancement, for instance. Another solution would be to implement SDA in the finite strain context combined with a mixed formulation, particularly suited to overcome locking effects when dealing with quasi-incompressible strain scenarios.

Conclusions and perspectives

The goals of this work were to adapt an existing theoretical framework, that has proven its capability in terms of objective description of brittle and quasi-brittle fracture, to the context of ductile fracture. One requirement was to propose a numerical method that limits the enhancements of existing Finite Element non-linear procedure to the element level, thus enabling the adaptation of this academic work to an industrial context which requires the use of non-intrusive procedures. This work relies on the Strong Discontinuity Analysis framework, initially developed in the brittle-fracture context. One of the main originalities of the presented work is the consideration of two coexisting dissipative mechanisms, a volumic dissipative mechanism and a surfacic dissipative directly resulting from the Strong Discontinuity Analysis framework, since it naturally involves a cohesive behavior. As for the volumic dissipative mechanism, few authors have considered a coupled damage-plasticity behavior such as Lemaitre's material model together with the Strong Discontinuity framework. This represents one of the major contributions of this work.

We first proposed the enhancement of a bilinear quadrangular element with a linear displacement jump (first presented by [START_REF] Linder | Finite elements with embedded strong discontinuities for the modeling of failure in solids[END_REF]) which is an extension of the more commonly used three-nodes triangular element that only describes constant strain fields. This element was never implemented in the context of ductile fracture considering a coupled damage-plasticity model in the bulk, and particularly together with Lemaitre's damage model. We proposed to appeal Lemaitre's damage variable as an introduction criterion of the discontinuity surface within a localized element. The proposed model allows to account for both mode I and II. We also presented a B-bar adaptation of the model to circumvent the locking phenomenon that we emphasized and issued by substantial plastic strain amounts for some cases.

The results presented demonstrated the good regularizing capabilities of the model. Indeed we presented accurate ultimate load predictions as well as consistent ultimate displacements, dissipated energies and crack paths.

Moreover, the so-called "cup-cone" phenomenon was obtained in the case of an axisymmetric computation. In particular, we emphasized that the transition from a plane to slant fracture mode (corresponding to a "cup-cone" in an axisymmetric framework) was, in our modeling context, driven by both the damage criterion D c and the fracture energy density G c . This enables interesting perspectives with regards to a possible identification of these material parameters. Besides, since these parameters have a physical sense and G c is a global and well-known fracture parameter, it most probably renders the identification process easier than for models that involve material parameters that do not have this obvious physical meaning.

Another objective was to propose an adaptation of the latter work to the context of large strain scenarios. This is of great importance when dealing with ductile fracture, which by definition is a phenomenon that involves strain amounts that are no longer consistent with a small strain hypothesis.

In this context, we proposed a numerical model that relies on the enhancement of a three node element with a constant displacement jump. The Finite Element implementation was performed in the current (deformed) configuration. This choice was motivated by the fact that it enables to appeal constitutive laws initially developed in the small strain context without major additional developments since, among other reasons, in this configuration the symmetries of the stress tensors are preserved.

The volumic behavior is based on the J2-plasticity theory adapted to the large strain framework, initially proposed by [START_REF] Simo | Computational inelasticity[END_REF]. As an introduction criterion for the surface of discontinuity, we used the Rice and Tracey criterion in the examples presented. Results emphasizing the regularizing effects of the method were presented, for mode I, mode II and mixed mode. We highlighted the influence, once again, of the fracture energy density parameter G c , this time underlying the fact that different distributions of G c between mode I and II involved different "dominating" modes (see section 3.4.2).

For all the developments involved in chapters 2 and 3, that is the Q4 enhanced element in small strain hypothesis, and the T3 element in finite strain hypothesis, it has to be noted that the structure of the program is such that any type of material model, for the bulk or the discontinuity surface, can be "plugged" in the element. Indeed, as long as a standard displacement Finite Element formulation is adopted, (thus excluding formulations in which the stress, strain and displacement are no longer dependent on each other through the compatibility equation and the material law, which would require a specific Finite Element formulation such as mixed or enhanced strain formulations), the elements that have been programmed can be seen as modules in which we can plug: a material law for the bulk, a cohesive law, and the necessary informations to introduce the discontinuity surface, that is a criterion and an orientation of the surface.

Hence, the perspectives of this work are quite numerous. For instance, it becomes straightforward to study fracture of a hyperelastic material in the finite strain range using the enhanced T3 element proposed in chapter 3. Evidently, the fracture criterion and the orientation of the crack still have to be carefully considered and studied. Even though this work was limited to study ductile fracture, we could have easily studied brittle fracture appealing the mentionned developments, which would have been out of the scope of this work as it was initially defined.

Since in the aeronautical industry, there is great interest of studying thin structures, another interesting perspective of the present work is the adaptation to plane stress hypothesis while only plane strain and axisymmetric problems were considered, the ultimate perspective being the possibility to consider full 3D problems.

We realize there are several limitations with regards to this work. For instance, for the Q4 element, Lemaitre's coupled-damage plasticity model was considered, and a rigorous criterion to introduce the discontinuity surface is the loss of ellipticity of the equilibrium equations, which corresponds to the singularity of the acoustic tensor defined as: ∃n /det(n.C an .n) = 0, where C an is the anelasic modulus (see [START_REF] Rice | A note on some features of the theory of localization of deformation[END_REF]). Such a condition is not straightforward for a complex material model such as Lemaitre's one, and authors have shown that this condition has an analytical formulation only for specific loadings, such as pure shear...etc... A rigorous solution to this problem would have been to compute the acoustic tensor for all the directions of space n at each time step (and each Gauss point) and to state that strain localization occurs when the acoustic tensor becomes singular, thus defining the moment, position and direction of strain localization, which would have induced prohibitive computational costs.

As a final word, we also emphasize the fact that the presented work involved a great amount of programming, at the Gauss point level for the bulk material for the coupled plasticity-damage model and the cohesive law, and especially at the element level for which solving the local equilibrium equation as well as obtaining the tangent stiffness matrix were not straightforward tasks for both elements. These also involved long debugging phases throughout this work period. and ω = 3G + ω(∆γ)H( ξn+1 ) 3 2 dev σtrial n+1 -σ y0 + q( ξn+1 ) (B.12)

For further information on the numerical integration of the material law, please refer to [START_REF] De Souza Neto | A fast, one-equation integration algorithm for the Lemaitre ductile damage model[END_REF].

It is clear that this consistent tangent operator is not symmetric due to the different expressions of c and d. As a consequence, a specific treatment is required to solve the global system, since the element tangent stiffness is also not symmetric. where we have omitted the term Γs since it disappears when writing the linearized equation h e + δh e = 0.

To completely compute this expression, we need to express the following quantities: δ FT , δ F-1 and δτ . Let us focus on the material model for the bulk to compute δτ . The time derivative of the Kirchhoff stress is defined as follows:

τ = c an T r : ε + lτ + τ l T (C.4)
where l is the deformation rate (or deformation gradient rate: l = ∇ u), and where c an T r is the tangent modulus defined in terms of the Trusdell rate of stress. It is equivalent to the Lie derivative of the Kirchhoff stress, and thus it is an objective definition of the Kirchhoff stress rate. Indeed, we have:

τ T r = τ -lτ -τ l T = c an T r : ε (C.5)
where τ refers to an objective stress rate, in this case, the Trusdell rate of the Kirchhoff stress. We could have used another objective stress rate such as the Jaumann stress rate (defined as: τ J = τ -Ωττ ΩT = c an J : ε where c an J is the Jaumann tangent modulus and Ω is the total spin tensor). The reason for the choice of one type of stress rate or another comes from the material model formulation. Since with our model, the output tangent modulus corresponds to the Trusdell rate of the Kirchhoff stress c an T r , we choose to use this one and derive our linearization considering this point. We will deal with this consideration in section 3.2.

The spatial strain rate is expressed as: ε = 1 2 (l + l T ). It follows that the terms lτ and τ l T are not material related but come from the geometry of the deformed configuration, as it is the case for the terms involving δ F-1 and δ F in (C.3).

With straightforward manipulations, a small variation of τ can be expressed in matrix form (for the 2D case): where ∇ x δu and ∇δu refer to the tensor forms of the spatial and material gradients of a small variation of total displacement δu, respectively.

Using equations in (C.7), we can write the following relations that will be useful to properly compute (C. After some matrix manipulations, we obtain the following formulation for the linearization of h e with respect to the total displacement u: We have also used the material model initially defined in (3.53). We recall its rate form expression: where both mode I and II are defined independently, and where the same cohesive model as in section 2.2.2 is used.

It remains to express the small variation of the total displacement δu and its spatial gradient ∇ x δu (and the symmetric part ∇ s x δu) in terms of small variations of nodal displacements δd and enhanced modes δα.

From the definition of the discretized total displacement in (3.30), we may write the small variation of the total displacement as: where ū is the spatial jump defined as : ū = F Ū where Ū is the material jump. Let us recall the multiplicative definition of the total deformation gradient in (3.12):

δu h (X) =
F = F I + Ū ⊗ Nδ Γs = F F (C.14)
The velocity gradient will allow us to define the gradient of a small displacement variation. It is defined as: With these expressions and following standard derivation procedures, we obtain the following expression for the vector form of the gradient of small variation of displacement: and its symmetric part: We have used the fact that the first integral that appears in h e deals with the domain Ω e \ Γ s . This is why the term in factor of δ Γs in (C.17) and (C.18) does not appear in this expression.

l = ḞF
(∇ x ) v δu =
(∇ s x ) v δu =

C.2 Linearization of f int,e

We follow the same procedure to obtain the linearization of the internal forces in (3.54). We can express f int,e as follows: Introducing the formulation of the gradient of a small variation of displacement (∇ x ) v δu defined in (C.17), we finally obtain the following formulation for a small variation of the internal forces with respect to d and α:

f int,e = Ω e

1. 1

 1 Ductile fracture physical mechanisms . . . . . . . . . . . . . . . . . . 1.2 Models for ductile fracture . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Micromechanical models . . . . . . . . . . . . . . . . . . . . . 1.2.2 Phenomenological approaches . . . . . . . . . . . . . . . . . . 1.3 Mesh dependency issue: a review of the localization limiters . . . . . 1.3.1 Highlighting the mesh dependency issue . . . . . . . . . . . . 1.3.2 Local methods . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Non-local methods . . . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Discontinuity models . . . . . . . . . . . . . . . . . . . . . . . Chapter 2 Ductile fracture modeling with the strong discontinuity method in small strain hypothesis 2.1 SDA: theoretical formulation and FE implementation . . . . . . . . . 2.1.1 Kinematics of the strong discontinuity approach (SDA) enhancement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Problem resolution . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3 B-bar method implementation along with SDA . . . . . . . . . 2.2 Material model formulation and implementation . . . . . . . . . . . . 2.2.1 Bulk material model : Lemaitre damage model . . . . . . . . . v Appendices Appendix A Formulation of the discretized problem List of Figures 1.1 (a) schematic plot of equation (1.5), (b) response of a material point in FE analyses with and without the modified function in 1.5 (see [178]) 1.2 Yield surface defined as the intersection of Thomason's model in (1.17) and the GTN model in equation (1.4) (thick line) . . . . . . . 1.3 Johnson-Cook's and Bai-Wierzbicki's models fracture loci . . . . . . . 1.4 (a) meshes and boundary conditions ("weakened" element in gray), (b) Reaction/Displacement curves . . . . . . . . . . . . . . . . . . . . 1.5 (a) compatible and incompatible shape functions, and displacement fields in the weak discontinuity framework of Huespe et al. [68], (b) compatible and incompatible shape functions, and displacement fields in the strong discontinuity framework . . . . . . . . . . . . .

  ), (b) fine mesh, (c) coarse mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.21 Force/displacement curve . . . . . . . . . . . . . . . . . . . . . . . . . 2.22 (a) damage variable, (b) triaxiality ratio, (c) accumulated plastic strain at the onset of plasticity . . . . . . . . . . . . . . . . . . . . . . 2.23 (a) damage variable, (b) triaxiality ratio, (c) accumulated plastic strain at the onset of localization . . . . . . . . . . . . . . . . . . . . 2.24 (a) Damage variable, (b) mode I crack opening in the first localized element, with respect to the prescribed displacement . . . . . . . . . 2.25 (a) Crack opening at the end of the test (b) Reaction/displacement curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.26 Geometry and boundary conditions of the shear specimen . . . . . . . 2.27 Crack opening in mode II at the end of the test in the deformed configuration (scale factor=1) for different mesh sizes . . . . . . . . . 2.28 (a) Vertical reaction/displacement curve without SDA (b) Vertical reaction/displacement curve with SDA . . . . . . . . . . . . . . . . . 2.29 (a) geometry of the specimen and symmetry simplifications, (b) FE mesh and discontinuity opening in mm at the end of the test . . . . . 2.30 Response for explicit formulation (blue cross) and implicit formulation (red square) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Motion decomposition scheme . . . . . . . . . . . . . . . . . . . . . . 3.2 Representation of bases (n , m ) and (n , m ), and projection of ū . . 3.3 Multiplicative motion decomposition scheme . . . . . . . . . . . . . . 3.4 Representation of the enhanced modes for the T3 triangle in its "undeformed" configuration (only the displacement jump is shown) . . . 3.5 Non-conforming shape function for the T3 triangle . . . . . . . . . . . 3.6 Multiplicative decomposition of the deformation gradient . . . . . . . 3.7 Mesh and boundary conditions . . . . . . . . . . . . . . . . . . . . . . 3.8 Responses with two different perturbation regularization procedures . 3.9 Responses for different values of viscosity parameters . . . . . . . . . 3.10 Meshes and boundary conditions . . . . . . . . . . . . . . . . . . . . 3.11 Reaction/displacement curves . . . . . . . . . . . . . . . . . . . . . . 3.12 Discontinuity opening in mode I on deformed meshes (×2) at the end of the computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 3.13 For three different mesh sizes, (a) mode I opening, (b) mode II sliding on deformed mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Reaction/displacement curves . . . . . . . . . . . . . . . . . . . . . . 3.15 Deformed meshes for the cases (a) G cII = G c and mode I is not allowed, (b) G cI = 0.2G c and G cII = 0.8G c . . . . . . . . . . . . . . . 3.16 Deformed meshes for the cases (a) G cI = 0.8G c and G cII = 0.2G c , (b) G cI = G c and mode II is not allowed . . . . . . . . . . . . . . . . 3.17 Reaction/displacement curves different fracture energy density distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18 Geometry and boundary conditions for the CT specimen . . . . . . . 3.19 Different meshes for the CT specimen . . . . . . . . . . . . . . . . . . 3.20 Reaction/displacement curves for the CT test . . . . . . . . . . . . . 3.21 Deformed mesh at the end of the CT test . . . . . . . . . . . . . . . .

Contents 1 . 1 2 1. 2 2 1. 2 . 1 3 1. 2 . 2 13 1. 3 24 1. 3 . 4

 11222213221332434 Ductile fracture physical mechanisms . . . . . . . . . . . Models for ductile fracture . . . . . . . . . . . . . . . . . . Micromechanical models . . . . . . . . . . . . . . . . . . . Phenomenological approaches . . . . . . . . . . . . . . . . Mesh dependency issue: a review of the localization limiters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.3.1 Highlighting the mesh dependency issue . . . . . . . . . . 20 1.3.2 Local methods . . . . . . . . . . . . . . . . . . . . . . . . 22 1.3.3 Non-local methods . . . . . . . . . . . . . . . . . . . . . . Discontinuity models . . . . . . . . . . . . . . . . . . . . . 28

Figure 1 . 1 :

 11 Figure 1.1: (a) schematic plot of equation (1.5), (b) response of a material point in FE analyses with and without the modified function in 1.5 (see [178])

Figure 1 . 2 :

 12 Figure 1.2: Yield surface defined as the intersection of Thomason's model in (1.17) and the GTN model in equation (1.4) (thick line)

Figure 1 . 3 :

 13 Figure 1.3: Johnson-Cook's and Bai-Wierzbicki's models fracture loci

Figure 1 . 4 :

 14 Figure 1.4: (a) meshes and boundary conditions ("weakened" element in gray), (b) Reaction/Displacement curves

Figure 1 . 5 :

 15 Figure 1.5: (a) compatible and incompatible shape functions, and displacement fields in the weak discontinuity framework of Huespe et al. [68], (b) compatible and incompatible shape functions, and displacement fields in the strong discontinuity framework

Figure 2 . 1 :Figure 2 . 2 :

 2122 Figure 2.1: Domain Ω crossed by a discontinuity surface Γ s

Figure 2 . 3 :

 23 Figure 2.3: α n0 parameter on Γ s (constant normal mode)

Figure 2 . 4 :Figure 2 . 5 :Figure 2 . 6 :

 242526 Figure 2.4: α n1 parameter on Γ s (linear normal mode)

Figure 2 . 7 :

 27 Figure 2.7: Incompatible shape function M (a) for constant normal mode, (b) for linear normal mode

  by solving the local (element) equation (2.27b) for fixed nodal displacements d (i)

Figure 2 . 8 :Figure 2 . 9 :

 2829 Figure 2.8: Meshes and boundary conditions for the three-point bending test

Figure 2 . 10 :

 210 Figure 2.10: Traction/normal opening law

Lemaitre 1 -Figure 2 . 11 :Figure 2 .

 12112 Figure 2.11: Geometry of the specimen and boundary conditions for the uniaxial tensile test

Figure 2 . 15 :

 215 Figure 2.15: Considered meshes and deformed mesh at the end of the test (x20)

Figure 2 . 16 :

 216 Figure 2.16: Reaction/displacement curve

Figure 2 . 17 :

 217 Figure 2.17: Considered meshes and boundary conditions for the three-point bending test

Figure 2 . 18 :

 218 Figure 2.18: Force/displacement curve when considering only Lemaitre's damage model

Figure 2 . 19 :

 219 Figure 2.19: Force/displacement curve for two different meshes, with and without SDA

Figure 2 .

 2 Figure 2.20: (a) geometry and boundary conditions (see [48]), (b) fine mesh, (c) coarse mesh

Remark 2 . 3 . 1 .

 231 Figure 2.21: Force/displacement curve

Figure 2 .

 2 Figure 2.22: (a) damage variable, (b) triaxiality ratio, (c) accumulated plastic strain at the onset of plasticity

Figure 2 .

 2 Figure 2.23: (a) damage variable, (b) triaxiality ratio, (c) accumulated plastic strain at the onset of localization

Figure 2 .Figure 2 .

 22 Figure 2.24: (a) Damage variable, (b) mode I crack opening in the first localized element, with respect to the prescribed displacement

Figure 2 . 26 :

 226 Figure 2.26: Geometry and boundary conditions of the shear specimen

Figure 2 . 27 :

 227 Figure 2.27: Crack opening in mode II at the end of the test in the deformed configuration (scale factor=1) for different mesh sizes

Figure 2 . 28 :

 228 Figure 2.28: (a) Vertical reaction/displacement curve without SDA (b) Vertical reaction/displacement curve with SDA

Figure 2 .Figure 2 . 30 :

 2230 Figure 2.29: (a) geometry of the specimen and symmetry simplifications, (b) FE mesh and discontinuity opening in mm at the end of the test

Contents 3 . 1 88 3. 1 . 1 96 3. 2 108 3. 3 116 3. 3 . 1 116 3. 4 126 3. 4 . 1 128 3. 4 . 2 132 3. 4 . 3 138 3. 5

 31881196210831163111641264112842132431385 SDA in the finite deformation range: theoretical formulation and FE implementation . . . . . . . . . . . . . . Theoretical formulation . . . . . . . . . . . . . . . . . . . 88 3.1.2 Finite element implementation . . . . . . . . . . . . . . . Material model formulation and implementation . . . . . 108 3.2.1 Bulk material model: Large strain von Mises plasticity model108 3.2.2 Theoretical formulation . . . . . . . . . . . . . . . . . . . Instability issue . . . . . . . . . . . . . . . . . . . . . . . . . Viscosity regularization . . . . . . . . . . . . . . . . . . . Numerical examples . . . . . . . . . . . . . . . . . . . . . . Simple tensile test in mode I fracture . . . . . . . . . . . . Simple tensile test exhibiting mixed mode fracture . . . . Compact Tension test (CT) . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

Figure 3 . 2 :

 32 Figure 3.2: Representation of bases (n , m ) and (n , m ), and projection of ū

Figure 3 . 4 :

 34 Figure 3.4: Representation of the enhanced modes for the T3 triangle in its "undeformed" configuration (only the displacement jump is shown)

Figure 3 . 5 :

 35 Figure 3.5: Non-conforming shape function for the T3 triangle

τ 11 0 τ 12 0 τ 12 0 τ 22 0 0 τ 11 0 τ 12 0 τ 12 0 τ 22

 22 = -a∈Ω + B a,x .[n m ], ḡns r = -a∈Ω + B ns a,x .[n m ], ḡv ns = -n ,ns .[n m ], τ =

Figure 3 . 6 :

 36 Figure 3.6: Multiplicative decomposition of the deformation gradient

[n m ] T .

Remark 3 . 3 . 1 .

 331 In equation (3.110) Can is different from Can (n,n) equation (3.116), since a projection has been performed between (3.110) and (3.111).

N

  T bJ dΩ + ∂tΩ e N T hJ dΓ h e = Ω e \Γs ḡT v τ | Ω e \Γs dΩ + Γs ḡv T t Γs + γ u)dΓ = -Γs [n m ] T n ,T τ | Ω e \Γs .n + Γs t n Critical time step and automatic time stepping implementation

2 ,UFigure 3 . 7 :

 237 Figure 3.7: Mesh and boundary conditions

Figure 3 . 8 :

 38 Figure 3.8: Responses with two different perturbation regularization procedures

Figure 3 . 9 :

 39 Figure 3.9: Responses for different values of viscosity parameters

Figure 3 . 10 :

 310 Figure 3.10: Meshes and boundary conditions

Figure 3 . 12 :

 312 Figure 3.12: Discontinuity opening in mode I on deformed meshes (×2) at the end of the computation

Figure 3 . 13 :Figure 3 . 14 :

 313314 Figure 3.13: For three different mesh sizes, (a) mode I opening, (b) mode II sliding on deformed mesh

Figure 3 . 15 :Figure 3 . 16 :

 315316 Figure 3.15: Deformed meshes for the cases (a) G cII = G c and mode I is not allowed, (b) G cI = 0.2G c and G cII = 0.8G c

Figure 3 . 18 :Figure 3 . 19 :Figure 3 . 20 :Figure 3 . 21 :

 318319320321 Figure 3.18: Geometry and boundary conditions for the CT specimen

  δh e = -[N M] T δ FT N x N y

C. 1 .

 1 Linearization of h eMoreover, we have the following expressions for δ F-1 and δ F:δ F = ∇δu δ F-1 = -F-1 ∇δu F-1 = -F-1 ∇ x δu (C.7)

F- 1 ∇ x δu 0 0 F- 1 11 2τ 12 τ 12 τ 22 τ 11 τ 12 

 10112 ] T δ FT = [N M] T (∇ x δu F) T = [N M] T FT (∇ x δu) T = [n m ] T (∇ x δu) T N x N y ∇ x δu = -n x n y 0 0 0 0 n x n y ∇ x δu 0 0 ∇ x δu (C.8) Inserting these expressions in (C.3) (only diagonal terms appear in the tangent modulus matrix), we can express the linearization of h e as follows:δh e = -[n n ] T (∇ x δu) T n ,ns,T H T τ + [n m ] T n ,ns,T ∇ x δu 0 0 ∇ x δu H T τ-[n m ] T n ,ns,T H T c an T

τ 11 0 τ 12 0 τ 12 0 τ 22 0 0 τ 11 0 τ 12 0 τ 12 0 τ 22

 22 δh e = -[n n ] T n ,ns,T τ + τ (∇ x ) v δu + c an T r (∇ s x ).9) and (C.10), (∇ x ) v denotes the vector form of the gradient. We have then (∇ x ) v u =

N

  a (X)δd a + H Γs (X) -a∈Ω + N a (X)]δ ū (C.13)

  -T Nδ Γs = l + u ⊗ n δ Γs (C.15) Thus we have identified u with the material jump rate U pushed forward in the current configuration through the relation u = F U. u actually corresponds to the Lie derivative of the material jump, thus it is an objective rate. We can now write the following expression for the small variation of the spatial jump: δ ū = F(δα n N + δα m M) = δα n n + δα m m (C.16)

  δd a + -a∈Ω + B ns a,x .[n m ] + δ Γs n ,ns .[n m ] δα (C.17

  )

BB

  a,x δd a + -a∈Ω + B a,x .[n m ] + δ Γs n .[n m ] δα(expressions (C.17) and (C.18) in (C.10), we obtain:δh e = -[n m ] T n ,ns,T τ + τ a,x δd + [n m ] T n ,ns,T τ + τ

21 )12 2τ 22 τ 12 τ 22 τ 11 τ 12 H

 2112 variation of f int,e is written as:δf int,e =Following the same definition as in (C.6), we have for the second integral in (C.21):(∇ x ) v δudΩ = Ω e B ns,T x H T c an T r + τ + τ (∇ x ) v δudΩ (C.22)where τ is defined in (C.11) and τ = first integral in (C.21), we can express B ns,Tx first integral in (C.21) can be expressed as: some matrix manipulations is finally expressed as:x ) v δudΩ (C.27) Gathering equations (C.27) and (C.22) into (C.21), we now obtain: δf int,e = Ω e B ns,T x H T c an T r + τ + τ (∇ x ) v δudΩ -T c an T r + τ (∇ x ) v δudΩ (C.28)

  .[n m ]dΩ δα (C.29)

  Meshes and boundary conditions for the three-point bending test . . 2.9 Force/displacement curve for two different meshes, with and without the B-bar method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Traction/normal opening law . . . . . . . . . . . . . . . . . . . . . . 2.11 Geometry of the specimen and boundary conditions for the uniaxial tensile test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 (a) meshes considered, (b) mode I opening at the end of the test . . . 2.13 (a) Reaction/displacement curves, (b) damage evolution in the weakened element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.14 Force/displacement curve for two different meshes, with and without SDA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 Considered meshes and deformed mesh at the end of the test (x20) . 2.16 Reaction/displacement curve . . . . . . . . . . . . . . . . . . . . . . . 2.17 Considered meshes and boundary conditions for the three-point bending test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 Force/displacement curve when considering only Lemaitre's damage model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.19 Force/displacement curve for two different meshes, with and without SDA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.20 (a) geometry and boundary conditions (see

. . . . . 2.1 Domain Ω crossed by a discontinuity surface Γ s . . . . . . . . . . . . 2.2 Schematic representation of the displacement fields for a 1D problem 2.3 α n0 parameter on Γ s (constant normal mode) . . . . . . . . . . . . . 2.4 α n1 parameter on Γ s (linear normal mode) . . . . . . . . . . . . . . . 2.5 α m0 parameter on Γ s (constant tangential mode) . . . . . . . . . . . . 2.6 α m1 parameter on Γ s (linear tangential mode) . . . . . . . . . . . . . 2.7 Incompatible shape function M (a) for constant normal mode, (b) for linear normal mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8

  Orientation of the discontinuity surface . . . . . . . . . . 58 2.2.5 Interpenetration processing . . . . . . . . . . . . . . . . . 59 2.3 Numerical results and discussion . . . . . . . . . . . . . . 59 2.3.1 Uniaxial tensile test . . . . . . . . . . . . . . . . . . . . . 60 Tensile test on a notched bar . . . . . . . . . . . . . . . . 72 2.3.5 Shear test . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 2.4 Extension to explicit formulation of the FE problem with SDA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 2.4.1 Critical time step . . . . . . . . . . . . . . . . . . . . . . . 80

	2.2.1 Bulk material model : Lemaitre damage model . . . . . . 53 2.2.2 Cohesive law on the discontinuity surface . . . . . . . . . 54 2.2.3 Criteria to introduce the surface of discontinuity . . . . . 57 2.3.3 Three-point centered bending test on a notched specimen 66 2.2.4 2.3.2 Embedded discontinuities with the B-bar method . . . . . 64 2.3.4

Contents 2.1 SDA: theoretical formulation and FE implementation . 36 2.1.1 Kinematics of the strong discontinuity approach (SDA) enhancement . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.1.2 Problem resolution . . . . . . . . . . . . . . . . . . . . . . 43 2.1.3 B-bar method implementation along with SDA . . . . . . 49 2.2 Material model formulation and implementation . . . . . 53 2.4.2 Comparison between explicit and implicit formulation . . 82 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 2.1 SDA: theoretical formulation and FE implementation 2.1.1 Kinematics of the strong discontinuity approach (SDA) enhancement

  Ḡr n0 Ḡr n1 Ḡr m0 Ḡr m1 and Ḡr = Ḡr n0 Ḡr n1 Ḡr m0 Ḡr m1

  Ḡr n0 Ḡr n1 Ḡr m0 Ḡr m1 and Ḡr = Ḡr n0 Ḡr n1 Ḡr m0 Ḡr m1 (2.45)

Table 2 .

 2 3: Computational cost for implicit and explicit formulations

	Total comput. time (s) Mean time increment (s) Final displ. (mm)
	Implicit case 83171	10 -3	0.78
	Explicit case 3709	∼ 10 -4	0.78

  .4, for the time step that corresponds to

	Iteration Residual norm	Energy norm
	7.759 × 10 3	2.713 × 10 2
	1.194 × 10 3	-1.968 × 10 1
	2.046 × 10 2	1.829 × 10 -1
	7.821 × 10 1	-7.425 × 10 -2
	2.122 × 10 2	1.968 × 10 -1
	7.822 × 10 1	-7.426 × 10 -2
	2.122 × 10 2	1.968 × 10 -1
	7.822 × 10 1	-7.426 × 10 -2
	2.122 × 10 2	1.968 × 10 -1
	7.822 × 10 1	-7.426 × 10 -2

Table 3 .

 3 2: Residuals during the Newton iterative procedure at the non-converged time step corresponding to the onset of localization

	Iteration Residual norm	Energy norm
	7.759 × 10 3	2.713 × 10 2
	1.171 × 10 3	-1.934 × 10 1
	2.061 × 10 2	1.856 × 10 -1
	7.781 × 10 1	-7.501 × 10 -2
	2.153 × 10 2	2.027 × 10 -1
	7.781 × 10 1	-7.501 × 10 -2
	2.153 × 10 2	2.027 × 10 -1
	7.781 × 10 1	-7.501 × 10 -2
	2.153 × 10 2	2.027 × 10 -1
	7.781 × 10 1	-7.501 × 10 -2

Table 3 .

 3

	3: Residuals during the Newton iterative procedure at the non-converged
	time step with viscous regularization

Table 3 .

 3 4: Residuals during the Newton iterative procedure at the time step at which discontinuities are introduced, with viscous regularization and auto-time stepping

  .8.

	Test type	visc. par.	Time inc. (s) Time steps numb. CPU time (s)
	Auto time	1 × 10 -1	∼ 7.65 × 10 -5	1.0515 × 10 4	1.085 × 10 1
	Auto visc. par. ∼ 1.2 × 10 1	1 × 10 -2	5.60 × 10 2	4.5 × 10 -1

Table 3 .

 3

	5: Comparison of automatic time increment procedure and automatic vis-
	cosity procedure

  .6 and figure 3.9). The end of computation is considered when the final displacement reaches 28 mm.

	Test numb. Visc. par. Av. time inc. (s) Time steps num. CPU time (s)
	1	1 × 10 -2	∼ 7.65 × 10 -6	1.045 45 × 10 5	9.955 × 10 1
	2	1 × 10 -1	∼ 7.65 × 10 -5	1.0515 × 10 4	1.085 × 10 1
	3	1	∼ 7.31 × 10 -4	1.107 × 10 3	1.28
	4	1.0 × 10 1	∼ 7.05 × 10 -3	1.61 × 10 2	1.9 × 10 -1
	5	1.00 × 10 2	∼ 6.80 × 10 -2	6.2 × 10 1	6 × 10 -2

Table 3 . 6

 36

: Computational results for different viscosity parameters

Remerciements

is the matrix of shape functions for a standard Q4 element, B(x) = [B 1 (x)...B Nn (x)] = LN(x) with L the matrix form of the ∇ s operator, and δu T (x, t) = [δu 1 (x, t)...δu Nn (x, t)] the vector of virtual displacements (N n being the number of nodes of the discretized problem), and with the discretized expression of virtual strains in (2.16), the discretization of the system gives: ∀(δu, δα),

Noticing that

with t Γs the traction vector on the surface Γ s , system (A.1) can be put into the following form:

Appendix B

Numerical implementation of Lemaitre's model

The numerical implementation of this model is carried out using the one equation return-mapping proposed by de Souza Neto et al. (see [START_REF] De Souza Neto | A fast, one-equation integration algorithm for the Lemaitre ductile damage model[END_REF]). Such an equation is solved using a Newton-Raphson procedure. This equation represents the discretized version of equation (2.50) and recalled in (B.1).

where ∆γ is the plastic multiplier that we seek,

and

where

G denotes the shear modulus, ε e,trial n+1 stands for the elastic trial strain, and ξn+1 = ξn + ∆γ.

In equation (B.1), Y (∆γ) is written as follows:

where K is the compressibility modulus and pn+1 = K 1 3 tr(ε n+1 ) is the hydrostatic part of the effective stress.

De Souza et al. in [START_REF] De Souza Neto | A fast, one-equation integration algorithm for the Lemaitre ductile damage model[END_REF] formulated the expression of the consistent tangent operator as: dev σ n+1 , I is the fourth-order identity tensor, 1 is the secondorder identity tensor, and

In the last expression, H( ξn+1 ) is the updated slope of the hardening curve, and together with (2.49), we obtain:

and a 1 , a 2 , a 3 , a 4 are defined as:

where

Appendix C

Linearization of the finite strain discretized problem

The following developments are based on a method proposed by Crisfield [START_REF] Crisfield | Non-Linear Finite Element Analysis of Solids and Structures[END_REF]. His method relies on the linearization of the already discretized problem (the FE formulation), and particularly the internal forces f int,e . Other methods rely on the linearization of the continuum problem (also described in the book written by Crisfield), and then the linearized continuum formulation is discretized, which is a method that we find less convenient.

C.1 Linearization of h e

Let's focus on the linearization of equation (3.55c), for which we recall the expression:

We may express h e differently as: With these expressions at hand, a small variation of h e is expressed as: