Gaël Micheline Laurent

Riccardo Juan

Benoit, Fabien, François, Julien, Malcolm Antoine Alexandre

Pierre Renaud

Sylvain

Keywords: Plasticity, dislocations, coarse-graining, correlations

REMERCIEMENT

Cette thèse a été pour moi une expérience formidable et incroyablement enrichissante. Cette périlleuse tâche, qui est de comprendre un petit peu plus la nature qui nous entoure, se solde souvent par des incompréhensions et de nouvelles questions. Cependant, cet exercice n'en reste pas moins une expérience incroyablement fascinante et stimulante. Cette passion des sciences, je la dois en particulier à mes encadrants dont ma gratitude ne s'arrête pas aux quelques phrases qui suivent. Ainsi, je remercie mon directeur de thèse Alphonse Finel pour m'avoir suivi attentivement et pour avoir essayé de me transmettre sa rigueur scientifique et ses connaissances. Je remercie également Yann Le Bouar pour ses nombreuses qualités qui ne se résument pas seulement à sa pédagogie, son écoute, sa patience ou encore ses précieux avis. De même, se fut un plaisir de travailler avec Benoît Appolaire à qui je dois beaucoup pour m'avoir éclairé et conseillé dans de nombreuses situations. Je remercie aussi Samuel Forest qui a su suivre à distance l'avancée de mon travail et apporter un oeil extérieur. Enfin, je remercie le jury pour l'intérêt qu'il a porté à mon travail.

Dans cette aventure, j'ai eu la chance de croiser des gens formidables, en particulier au Laboratoire d'Etude des Microstructures, qui m'ont accueilli dans une très bonne ambiance. Je remercie Benoit Devincre et Ladislas Kubin pour les discussions intéressantes sur la plasticité ; Jean-Sebastien pour ses nombreux coups de main ; Frédéric pour ses emails distrayants ; Denis Gratias pour sa passion des sciences qu'il transmet naturellement ; Hakim Amara pour sa sympathie. Je dois aussi beaucoup à Léonard qui a dû subir les discussions sur les dislocations dans le bureau. Je tiens également à remercier Pierre-Antoine et Antoine pour les échanges que l'on a pu avoir sur nos problématiques communes ainsi que Luc pour ses phonons et son soutien lors de la rédaction. Je remercie aussi

INTRODUCTION I.1 Context

-2 dimension (or surface defect): interfaces, grain boundaries, stacking faults, cracks, etc. Each defect is governed by a specific kinetics and dislocation movement has been shown to be the main mechanism of the plasticity. This phenomenon is at the heart of the present thesis.

A dislocation, illustrated in figure I.1, is characterized by two quantities: -The Burgers vector characterizes the plastic shift of a part of the crystal during the movement of the dislocation. For usual perfect dislocations, the Burgers vector is a translation vector of the crystal and is often along a dense direction of the crystal.

-The dislocation line represents the position of the defect. The line is locally oriented along a unit line vector l. 2 ) between two phases with different orientations. On the top, there is a vacancy (dim. 0 ). On the right, there is a portion of dislocation loop (dim. 1 ) enclosing an area sheared by an elementary translation b, the Burgers vector. l is the dislocation line vector. Drawing inspired by [START_REF] Read | Dislocations in Crystals[END_REF].

The movement of a dislocation inside a plane containing the direction of the Burgers vector is called glide. Glide planes (or slip planes) are usually the densest planes of the crystal. For example, in a FCC material illustrated in figure I.2, the four {111} planes are the commonly observed slip planes.

Figure I.2: Lattice structure of a FCC material with slip-systems (dense planes) illustrated by grey planes. Adapted from [START_REF] Forest | Computational Micromechanics of Single Crystals and Polycrystals[END_REF].

The description of crystal plasticity is therefore usually based in slip systems, defined by a slip plane and a slip direction (i.e. the direction of the Burgers vector). At low temperature, dislocation movement is only due to glide. At high temperature, screw dislocations may change their slip plane (cross slip) and dislocation may even move out of their slip plane by absorbing or releasing vacancies (dislocation climb) [START_REF] Geslin | Contribution à la modélisation champ de phase des dislocations[END_REF]. These two thermally activated processes are not considered in this study.

The description of dislocation in text books usually starts by considering straight dislocations. When b is perpendicular to the dislocation line vector l, the dislocation is called an edge dislocation. When b and l are parallel the dislocation is of screw character and in all other cases the dislocation has a mixed character. The glide of the dislocation implies a series of ruptures of atomic bonds followed by reattachment with the next atomic column. Due to the crystal periodicity, a critical resolved stress, called Peierls stress, must be overcome in order to make the dislocation glide. In this work, this Peierls stress is considered negligible, as it is the case in FCC crystals due to the dissociation of the dislocation core. In addition, the applied stress will be moderate, so that the velocity of the dislocation can be assumed proportional to the local resolved stress • defined by

• = b b • ‡ • n (I.1)
where n is the normal to the glide plane.

We now consider the case of an assembly of dislocations in a crystal. Because each dislocation generates strain fields in the crystal, a given dislocation will receive from the other dislocations a stress field that may lead to the movement of the dislocation. The elastic fields generated by a dislocation are decaying very slowly (as the inverse of the distance to the dislocation) and are also usually anisotropic. In addition, dislocations interact at short range (in particular by the formation of junctions), and the dislocation density may strongly evolve in time (source of dislocations).

All this makes it clear that the understanding of the collective behavior of dislocations and the emergence of complex patterns is a very challenging problem. An example of such a dislocation pattern in a monocrystal of copper is shown in figure I. [START_REF] Acharya | A model of crystal plasticity based on the theory of continuously distributed dislocations[END_REF]. The obtained pattern depends strongly on the time evolution and orientation of the applied load [START_REF] Takeuchi | Work Hardening of Copper Single Crystals with Multiple Glide Orientations[END_REF][START_REF] Kocks | Physics and phenomenology of strain hardening: the FCC case[END_REF][START_REF] Staker | The dislocation cell size and dislocation density in copper deformed at temperatures between 25 and 700 ¶[END_REF]. Scale invariant (fractal) patterns are also reported [START_REF] Zaiser | Fractal analysis of deformation-induced dislocation patterns[END_REF][START_REF] Zaiser | The flow stress of fractal dislocation arrangements[END_REF]. Revealed by typical transmission micrographs. Respectively after [START_REF] Kocks | Physics and phenomenology of strain hardening: the FCC case[END_REF][START_REF] Glazov | Low energy dislocation structures produced by cyclic deformation[END_REF].

The plasticity of heterogeneous crystals or polycrystals may be even more complex because the heterogeneities of the material introduce new length scales in the problem. This point is illustrated by the tensile test of a polycrystal of aluminium presented in figure I.5. After deformation, the orientation map of the polycrystal shows that the heterogeneity of the plastic field appears at the scale of the grain structure but also at a scale much lower than the grain size. After [START_REF] Zhao | Investigation of threedimensional aspects of grain-scale plastic surface deformation of an aluminum oligocrystal[END_REF].

From these examples, it appears clearly that the macroscopic mechanical behavior of a material cannot be understood and predicted without a deep analysis of all the intermediate length scale (called mesoscopic scales) between the micro and macro scales. The main difficulty is probably to derive relevant models to study the mechanical behavior at mesoscale. These models should in principle describe the transport of dislocations but also the consequences of the short and long range interactions between dislocations, the creation and annihilation of dislocations, etc. Such a complete physically based mesoscopic model is still far from being available in the scientific community.

I.3 Crystal plasticity models

The investigation of crystal plasticity goes through different types of tools specific for each scale, ranging from atomic scale to the macroscopic one: atomistic simulations, Discrete Dislocation Dynamic (DDD), Field Dislocation Mechanics (FDM) or Continuum Dislocation Dynamic (CDD) models, based on dislocation densities, and finally, conventional plasticity theories. Atomistic simulations follow the behavior of each atom of the lattice and are useful for the description of dislocation core interacting with different defects (interfaces, cracks, etc.). This is the most rigorous description of a dislocation but also the heaviest one and only a few dislocations can be simultaneously simulated. However, atomistic analysis can be used to calibrate DDD models. Even though less rigorous, it is currently the best way to understand the collective behavior of dislocations. The drawback is that simulations are limited in time and also in space (few thousands of dislocations) making it impossible to simulate complex microstructures. In fact, this tool is limited to simple crystals or, at least, to situations where there are only a few simple grain boundaries or interfaces, to analyse the impact of these defects on the dislocation dynamics [START_REF] Daveau | Interaction dislocations -joints de grains en déformation plastique monotone : étude expérimentale et modélisations numériques[END_REF]. It is worth mentioning that this description at the level of discrete dislocations can also be reached by a continuous Phase Field formulation based on the so-called "loopons" [START_REF] Rodney | Dislocations in crystals[END_REF].

The next scale, the mesoscopic scale, is the scale of our investigations in this thesis. Industrial development of polycrystalline materials are often approached by phenomenological crystal plasticity models based on macroscopic concepts and laws (Norton, Von Mises, etc.), see for example [START_REF] Gaubert | Modélisation des effets de l'évolution microstructurale sur le comportement mécanique du superalliage monocristallin AM1[END_REF]. Consequently, these models are built without taking into account the origin of plasticity, the dislocations. These models are sometimes improved with a gradient term of plastic strain or stress [START_REF] Forest | Stress gradient continuum theory[END_REF] in order to incorpore the size effect [START_REF] Cottura | Modélisation du champ de phase du couplage entre évolution microstucturale et propriétés mécaniques[END_REF] observed in the plastic behavior when plasticity develops in regions which are below a few microns in size. A better way is to approach the mesoscopic plasticity by handling dislocation densities.

A lot of works have already been done on dislocation density models. First studies were done a long time ago before the help of numerical resources. One of the first contribution was the Taylor law (1934). Taylor proposed that tangled dislocations are able to move only if the applied stress exceeds a critical value, called • c , which is inversely proportional to the distance between dislocation and thus proportional to Ô fl:

• c = -µb Ô fl (I.2)
where fl is the dislocation density, -is a phenomenological coefficient, µ is the shear modulus and b is the norm of the Burgers vector. This law was further generalized to the case of multiple glide systems [START_REF] Franciosi | Latent hardening in copper and aluminium single crystals[END_REF] and was then used in a phenomenological model developed to analyse each stage of the strain hardening of a monocrystal [START_REF] Devincre | Physical analyses of crystal plasticity by DD simulations[END_REF][START_REF] Kubin | Toward a physical model for strain hardening in ccc crystals[END_REF]. The next important contribution was the Kocks-Mecking model which proposes to define two types of dislocation density: the stored density and the mobile density [START_REF] Kocks | Physics and phenomenology of strain hardening: the FCC case[END_REF][START_REF] Kubin | Modeling dislocation storage rates and mean free paths in face-centered cubic crystals[END_REF][START_REF] Kubin | Dynamic recovery and its orientation dependence in facecentered cubic crystals[END_REF]. This model is built on the concept of the mean free path of a dislocation to reproduce dislocation storage. It integrates also the annihilation phenomenon between dislocations which leads to the restoration of the crystal. The Kocks-Mecking model has been extended to the case of crystal plasticity [START_REF] Teodosiu | Finite element simulation of the large elastoplastic deformation of multicrystals[END_REF] and has been used, for example, to study the mechanical behavior of polycrystal [START_REF] Teodosiu | Finite element simulation of the large elastoplastic deformation of multicrystals[END_REF] or Nickel based superalloy [START_REF] Cottura | Modélisation du champ de phase du couplage entre évolution microstucturale et propriétés mécaniques[END_REF][START_REF] Vattré | Strength of single crystal superalloys: from dislocation mechanisms to continuum micromechanics[END_REF], etc. Other formulations have been proposed to describe plasticity at the mesoscale such as FDM or CDD models. These formulations will be reviewed in the next chapter.

The important point is that all these models, describing plasticity using dislocation densities, require a phenomenological constitutive law to describe the crystallographic slip and the slip resistance. When changing the material under study or the type of loading, it is necessary to recalibrate these laws and sometimes to even change the law itself.

Therefore, it appears that it is important to work towards the definition of more physically based plasticity model, and the most relevant way to do that is to directly derive a mesoscopic model by coarse-graining the dynamics of individual dislocations. This route was opened in 1997 by I. Groma who proposed to use the tools of statistical physics to derive a dislocation density based model from the coarse-graining of the 2D dynamics of parallel edge dislocations.

Even if the geometry of the plastic model is very restrictive, this approach was the first one to be able to derive transport equation for dislocations densities. This approaches has also the advantage of pointing out the importance of the correlations of the dislocation positions in this coarse-graining procedure [START_REF] Groma | Origin of gradient terms in plasticity at different length scales[END_REF][START_REF] Groma | Dynamics of coarse grained dislocation densities from an effective free energy[END_REF]. Most importantly, the correlation term can be split into two contributions which are formally similar to the usual friction and backstress terms [START_REF] Groma | Spatial correlations and higher-order gradient terms in a continuum description of dislocation dynamics[END_REF], phenomenologically introduced in the classical plasticity models. Several papers have been devoted to the understanding and the computation of these correlation terms [START_REF] Zaiser | Statistical dynamics of dislocation systems: The influence of dislocation-dislocation correlations[END_REF]. The model was used to investigate dislocation patterning [START_REF] Groma | Investigation of dislocation pattern formation in a two-dimensional selfconsistent field approximation[END_REF][START_REF] Groma | Dislocation Patterning: From Micro-to Mesoscale Description[END_REF], and extended to the case of multiple slip systems [START_REF] Yefimov | Multiple slip in a strain-gradient plasticity model motivated by a statistical-mechanics description of dislocations[END_REF].

This manuscript reinvestigates this formalism and presents my contribution to this model. The first chapter details the coarse-graining procedure and aims at clarifying each step of the derivation, starting from the dynamics of discrete dislocations up to transport equations of dislocation density. In the second chapter, dislocation dynamics is used to estimate the correlation terms that emerges from the coarse-graining procedure. In the third chapter, the numerical implementation of the dislocation density model is discussed and an improved numerical scheme is proposed. Then the model is used to address the question of the spontaneous emergence of a dislocation density pattern when starting from an initial homogeneous density. The fourth chapter proposes a new FFT based scheme for numerically solving mechanical equilibrium. It is shown that this scheme, which could be used to improve the numerical solution of the dislocation density model, is able to provide accurate and non oscillatory mechanical field even in the case of very inhomogeneous anisotropic elasticity. Finally, the important results of the work are summarized and several perspectives are proposed.

Chapter 1 DENSITY BASED CRYSTAL PLASTICITY: A DERIVATION OF COARSE-GRAINED TRANSPORT EQUATIONS

Résumé du chapitre

Ce chapitre explique la mise en oeuvre d'une méthode de changement d'échelle permettant de passer de la dynamique des dislocations discrètes (DDD) à une théorie en densité de dislocations dans le cas d'une assemblée de dislocations coins parallèles se déplaçant suivant un seul système de glissement. Cette méthode qui reprend la dérivation initialement proposée par I. Groma [START_REF] Groma | Link between the individual and continuum approaches of the description of the collective behavior of dislocations[END_REF][START_REF] Groma | Spatial correlations and higher-order gradient terms in a continuum description of dislocation dynamics[END_REF], aboutit d'une part à des équations de transport sur les densités de dislocations et d'autre part aux différents termes contribuant à l'évolution de ces densités. Dans le cadre de la situation simple considérée, cette méthode permet de s'affranchir de la phénoménologie des modèles de plasticité cristalline utilisés en mécanique des matériaux. Après un état de l'art rappelant brièvement les travaux sur lesquels nous nous sommes appuyés, nous présenterons en détail la procédure de changement d'échelle. Nous commencerons ainsi par définir rigoureusement les densités de dislocations à l'aide de moyennes spatiale, temporelle et d'ensemble. Puis nous discuterons les différentes contributions à la force motrice d'évolution des densités de dislocations, à savoir la contrainte de champ moyen et les contraintes induites par les corrélations entre dislocations individuelles. Nous serons alors à même de proposer une signification physique claire à la contrainte dite 'backstress'. Dans un second temps, nous discuterons la brisure de symétrie émergeant dans les équations en densité de dislocations suite au changement d'échelle.

Introduction

Plasticity of crystalline solids involves the notion of dislocations. However, even today, conventional plasticity theories use mesoscopic variables and evolution equations that do not involve dislocations. This paradoxical situation is due to the enormous length and time scales that separate the description of plasticity at the level of individual dislocations and the macroscopic scale of engineering materials. This huge space and time separation renders the hope to use a discrete dislocation based approach out of reach for treating engineering problems. It could be argued that conventional or phenomenological plasticity theories are justified because, at the macroscopic scale, engineering materials always display some sort of disorder that gives to any macroscopic property or measure an inevitable averaging character. Hence, at the macroscale, plastic strain may be seen as resulting mesoscale from a space and time average over a huge number of individual dislocation glide events.

Nevertheless, conventional plasticity theories rely on strong approximations and on phenomenological laws that must be calibrated for each material, or, for each specific applications. Therefore, it is desirable to make a link between the micro and macro scales and to develop a mesoscopic plasticity theory that relies on a sound physical basis, i.e. that at least incorporates dislocation glide. The development of such a mesoscale theory is also crucial to better understand and simulate the materials behavior at length scales where the elastic interaction between dislocations becomes of the order of the interaction between dislocations and obstacles, such as precipitates in a matrix, small grains in a polycrystal or interfaces in nano-materials. At these scales, dislocations display collective phenomena that result in patterning and complex dynamic regimes. In these situations, plasticity cannot be described by a simple averaged plastic strain that obeys local time-dependent equations. Size-dependent effects and, most importantly, transport become fondamental. Conventional theories of plasticity are no longer valid and are unable to account for the complexity of the plastic activity because they lack the relevant internal length scale and do not incorporate transport.

These considerations motivate the development of continuum models in which dislocations are represented by continuous densities and in which the dynamics has conserved the transport character of the underlying dislocation glide.

Continuum dislocation representations often start from the Nye [START_REF] Nye | Some geometrical relations in dislocated crystals[END_REF] and Kröner [START_REF] Kröner | Allgemeine Kontinuumstheorie der Versetzungen[END_REF] representation of dislocations. This is the case of the Field Dislocation Model (FDM) proposed by Acharya [START_REF] Acharya | A model of crystal plasticity based on the theory of continuously distributed dislocations[END_REF][START_REF] Acharya | Driving forces and boundary conditions in continuum dislocation mechanics[END_REF] and developed subsequently by various authors [START_REF] Roy | Finite element approximation of field dislocation mechanics[END_REF][START_REF] Fressengeas | An elasto-plastic theory of dislocation and disclination fields[END_REF][START_REF] Taupin | Grain boundary modeling using elasto-plastic theory of dislocation and disclination fields[END_REF][START_REF] Fressengeas | A field theory of distortion incompatibility for coupled fracture and plasticity[END_REF]. The basic equations have been in fact known as early as the 60's [66,83] (see also volume 7, Theory of Elasticity, in the series Theoretical physics, by Landau-Lifshitz [START_REF] Landau | Theory of Elasticity[END_REF] and the chapter Crystal dislocations and the theory of elasticity by A. Kosevich in "Dislocation in Solids" [START_REF] Kosevich | C r y s t a l dislocations and the theory of the elasticity[END_REF]. The basic ingredient of the FDM is the dislocation density tensor -= ≠ curl -p ,where-p is the plastic distorsion tensor. When envisaged at the smallest scale, the tensor -represents all the dislocations and there is no need to introduce the concept of "geometrically necessary" or "statistically stored" dislocations (GND and SSD, respectively). The model is then exact, regardless of the atomic nature of the dislocations and provided that we accept that the dislocation velocity is simply proportional to the local resolved shear stress. However, being continuous by nature, the implementation of the model requires the use of a computational grid with a grid step significantly smaller than the Burgers vector length. This drastically limits the spatial length scale that can be investigated. Therefore, in order to reach a convenient macro scale, a change of scale must be performed to bridge the gap between the singular density tensor introduced above and a continuous one defined at an intermediate scale. There is of course no unique way to select this so-called "mesoscale". Obviously, the mesoscale must be larger than the average distance between dislocations and smaller than the characteristic length scale we want to investigate (average grain size in polycrystals, average distance between interfaces in multiphase alloys, etc.). The underlying averaging or "coarse-graining" procedure has of course been already mentioned in the context of the FDM [START_REF] Acharya | Size effects and idealized dislocation microstructure at small scales: Predictions of a Phenomenological model of Mesoscopic Field Dislocation Mechanics: Part I[END_REF][START_REF] Fressengeas | Computational Methods for Microstructure-Property Relationship[END_REF].

The crucial point is that the application of the coarse-graining procedure to the FDM equations leads to transport equations for the averaged one-body GND density in which the plastic strain rate inevitably depends on the correlations between the lower scale GND and velocity fields. This closure problem is often resolved by using a phenomenological velocity law borrowed from macroscopic plasticity models leading to the so-called Phenomenological Mesoscopic Field Dislocation Model (PMFDM) [START_REF] Acharya | Size effects and idealized dislocation microstructure at small scales: Predictions of a Phenomenological model of Mesoscopic Field Dislocation Mechanics: Part I[END_REF][START_REF] Roy | Size effects and idealized dislocation microstructure at small scales: Predictions of a Phenomenological model of Mesoscopic Field Dislocation Mechanics: Part II[END_REF]. The actual implementation of the mesoscale FDM thus suffers from the lack of a mathematically justified mesoscale plastic strain rate.

A more recent formulation of a Continuum Dislocation Dynamics (CDD) has been proposed by Hochrainer and its collaborators [START_REF] Hochrainer | Evolving Systems of Curved Dislocations: Mathematical Foundations of a Statistical Theory[END_REF][START_REF] Hochrainer | A three-dimensional continuum theory of dislocation systems: kinematics and mean-field formulation[END_REF]. It is based on a modified definition of the dislocation density tensor, in order to keep at mesoscale information concerning the geometry of the dislocations (in particular, line directions and curvatures). The necessity of using an averaging procedure to obtain a meaningful continuum model has also be pointed out in the context of the CDD formulation [START_REF] Sandfeld | Continuum modeling of dislocation plasticity: Theory, numerical implementation, and validation by discrete dislocation simulations[END_REF] (see also [START_REF] El-Azab | Statistical mechanics treatment of the evolution of dislocation distributions in single crystals[END_REF][START_REF] El-Azab | Statistical mechanics of dislocation systems[END_REF]), but a rigorous mathematical formulation of this coarse-graining procedure has not yet been proposed.

The first attempt to better treat the closure problem has been proposed by Groma [START_REF] Groma | Link between the microscopic and mesoscopic length-scale description of the collective behavior of dislocations[END_REF] and its collaborators [START_REF] Zaiser | Statistical dynamics of dislocation systems: The influence of dislocation-dislocation correlations[END_REF][START_REF] Groma | Spatial correlations and higher-order gradient terms in a continuum description of dislocation dynamics[END_REF]. This is the route that we follow below. A particular attention will be paid on the nature of the coarse-graining procedure and its consequences on the local stress fields that emerge from the time and space averaging process.

Density based plasticity theory: from the discrete to the continuum

Our aim here is to clarify the mathematics and physical aspects of the coarse-graining procedure that must be used to coarse-grain the dislocation dynamics from discrete to the continuum. Therefore, we consider the simplest situation, namely a 2D dislocation system with N edge dislocation lines parallel to the z-axis restricted to glide along the x-axis. The Burgers vector of dislocation i, i =1 to N , is noted s i b,w h e r es i is the sign of the dislocation i and b =(b, 0, 0). The numbers N + and N ≠ of negative and positive dislocations are considered equal and constant. We assume an overdamped motion: the glide velocity of the i th dislocation along the x-axis is simply proportional to the resolved Peach-Koehler force acting on the dislocation i,

dr i dt = Ms i b Q a N ÿ j" =i s j • ind (r i ≠ rj )+• ext R b , (1.1)
where M is the mobility coefficient equal to the inverse of the dislocation drag coefficient, • ext the external stress resolved in the slip system and • ind (r) the shear stress at position r generated by a positive dislocation located at the origin:

• ind (x, y)= µb 2fi(1 ≠ ‹) x(x 2 ≠ y 2 ) (x 2 + y 2 ) 2 (1.2)
where µ is the shear modulus and ‹ the Poisson coefficient. The first step is to define discrete dislocation densities:

fl + dis (r, t, {r 0 k })= q N i=1 " s i ,+1 "(r ≠ ri (t, {r 0 k })) fl ≠ dis (r, t, {r 0 k })= q N i=1 " s i ,≠1 "(r ≠ ri (t, {r 0 k })) (1.3)
where {r 0 k } refers to the initial positions of the N dislocations, " s,t is the Kronecker symbol and "(r) the 2D Dirac function. The notation ri (t, {r 0 k }) means that the trajectory of dislocation i depends on the initial dislocation positions {r 0 k }. By multiplying equation (1.1) by the Dirac function "(r ≠ ri (t, {r 0 k })) and taking its derivative with respect to r, we get the following transport equation for the discrete densities:

≠ ˆt fl s dis (r)=sM b • ˆr Y _ ] _ [ ⁄ r Õ " =r • ind (r ≠ r Õ ) ÿ s Õ =±1 s Õ fl s Õ dis (r Õ )fl s dis (r)dr Õ + • ext fl s dis (r) Z _ _ \ (1.4)
where, to simplify the notation, we write fl s dis (r) for fl s dis (r, t, {r 0 k }). Obviously, these transport equations link the time-dependence of the one-body densities to the product of two one-body densities, which is a direct consequence of the pairwise dislocation interactions. At this stage, the dislocation densities fl s dis (r) are highly singular. The next step is to introduce a coarse-graining procedure.

Coarse-graining procedure

We introduce now a coarse-graining procedure commonly used in statistical physics. We first define a space and time convolution window w(r, t) that we use to coarse-grain microscopic fields to mesoscopic ones:

f meso (r, t)= ⁄⁄ w(r ÕÕ ,t ÕÕ ) f micro (r + r ÕÕ ,t+ t ÕÕ ) dr ÕÕ dt ÕÕ .
(1.5)

The weighting function w(r, t) is non-dimensional and normalized. For simplicity, and without loss of generality, we choose w(r, t) to be separable:

w(r, t)=w L (r) w T (L) (t) (1.6)
where the functions w L (r) and w T (L) (t) are separately normalized:

⁄ w L (r) dr =1and ⁄ w T (L) (t) dt =1. (1.7)
The spatial linear dimension L of w L (r) should be of the order of the spatial resolution of the continuous model we seek and, obviously, significantly larger than the average distance between dislocations. The temporal width T (L) of the time window w T (L) (t) should, in all generality, depend on L (we comment on that point in section 1.2.4) and, for convenience, we choose w T (L) (t) to be non-zero only for t 6 0:

w T (L) (t) " =0 if t 6 0. (1.8)
Mesoscopic density fields may be defined through equation (1.5), but this is not enough to get a consistent continuous transport theory. First, we expect that the time evolution of the mesoscopic dislocation densities will be given by first-order transport (i.e. hyperbolic) equations. These equations must be supplemented by initial conditions at t =0which, of course, must be defined at mesoscale. In other words, the coarse-graining procedure should be such that, when applied to equation (1.4) and its initial condition given by the dislocation positions {r 0 k } at t =0, we end up with a set of mesoscopic transport equations supplemented by continuous initial conditions that do not depend on any specific initial set {r 0 k }.T h e r e f o r e ,i ffl s (r, t= 0), s = ±1, are given initial continuous densities, we must introduce a N -body probability density distribution P (r 0 1 ,...,r 0 N ) on the (discrete) initial positions {r 0 k } which is linked to the mesoscopic densities fl s (r, t= 0) in a way that we discuss below. The distribution P (r 0 1 ,...,r 0 N ) introduces a statistical ensemble on the initial discrete dislocation positions: P (r 0 1 ,...,r 0 N )dr 0 1 ...dr 0 N is the probability to have an initial dislocation configuration with dislocation 1 in a small volume dr 0 1 around position r 0 1 , dislocation 2 in a small volume dr 0 2 around position r 0 2 ,e t c . Now, the overall coarse-graining procedure is defined as the conjugate action of the space-time convolution window w(r, t) and the ensemble average defined by the probability density P (r 0 1 ,...,r 0 N ). The mesoscopic field X meso (r, t) associated with the discrete field X dis (r, t, {r 0 k }) is therefore defined by:

X meso (r, t)= N Ÿ k=1 ⁄ dr 0 k P (r 0 1 ,...,r 0 N ) ⁄ dr ÕÕ ⁄ dt ÕÕ w(r ÕÕ ,t ÕÕ )X dis (r + r ÕÕ ,t+ t ÕÕ , {r 0 k }). (1.9)
We refer to this coarse-graining procedure by the following short-hand notation:

X meso (r, t)=ÈÈX dis (r, t)ÍÍ P (1.10) where the double brakets refers to the space and time convolution and the lower index P to the ensemble average. The mesoscopic one-body and two-body densities are therefore defined by:

fl s (r, t)=ÈÈfl s dis (r, t, {r 0 k })ÍÍ P (1.11) and fl ss Õ (r, r Õ ,t)=ÈÈfl s dis (r, t, {r 0 k })fl s Õ dis (r Õ ,t,{r 0 k })ÍÍ P . (1.12)
We note that the two-body densities defined in equation (1.12) are continuous function of r and r Õ . This would not be the case if the coarse-graining procedure was limited to a space and time convolution. This is the second reason why we need to consider also an average over a statistical ensemble.

We can now precise the link, mentioned above, between the probability density P (r 0 1 ,...,r 0 N ), that defines the statistical ensemble, and the continuous dislocation densities fl s (r, t) that will be used as initial conditions for the mesoscopic kinetic equations. We consider that any discrete initial condition {r 0 k } on the N dislocation positions is extended to t<0: i =1to N and t 6 0: ri (t, {r 0 k })={r 0 i }.

(1.13)

Then, using the definition of the discrete densities (equation (1.3)) and the definition of the coarse-grained ones (equation (1.11)), we get:

fl s (r, t = 0) = N Ÿ k=1 ⁄ dr 0 k P (r 0 1 ,...,r 0 N ) ⁄ dr ÕÕ ⁄ dt ÕÕ w(r ÕÕ ,t ÕÕ ) N ÿ i=1 " s i ,s "(r + r ÕÕ ≠ ri (t ÕÕ , {r 0 k })).
(1.14) Using equations (1.6), (1.7) and (1.13), we obtain:

fl s (r, t = 0) = N ÿ i=1 " s i ,s N Ÿ k=1 ⁄ dr 0 k P (r 0 1 ,...,r 0 N ) w L (r 0 i ≠ r). (1.15)
Without loss of generality, we may restrict the stochastic variables r 0 i , i =1 to N ,t obes t a t i s t i c a l l y independent and to follow the same distribution function f (r).T h u s ,t h ed e n s i t yP (r 0 1 ,...,r 0 N ) may be factorized as:

P (r 0 1 ,...,r 0 N )=f (r 0 1 )f (r 0 2 ) ...f(r 0 N ) (1.16)
where the density f (r) is normalized:

⁄ f (r)dr =1.
(1.17)

Using equations (1.16) and (1.17), equations (1.15) become

fl s (r, t = 0) = N ÿ i=1 " s i ,s ⁄ w L (r 0 i ≠ r)f (r 0 i )dr 0 i (1.18)
and then fl s (r,

t = 0) = N s ⁄ w L (r 0 ≠ r)f (r 0 )dr 0 (1.19)
where N s is the number of dislocations of the sign s.U pt ot h ec o e fficient N s , the initial condition fl s (r, t = 0) is simply equal to the convolution of f (r), the distribution of initial discrete dislocation positions, with the convolution window w L (r). For given fl s and w L , equation (1.19) defines a unique function f (r). Thus, for prescribed initial mesoscopic dislocation densities fl s (r, t= 0) and a given spatial convolution window w L (r), the coarse-graining procedure introduced in equation (1.9) is completely and uniquely defined.

Coarse-grained kinetic equations

By a direct application of the coarse-graining procedure defined in equation (1.9) to equation (1.4), we get the following mesoscopic equations:

≠ ˆt fl s (r, t)=sM b • ˆr Y _ ] _ [ ⁄ r Õ " =r • ind (r ≠ r Õ ) ÿ s Õ s Õ fl ss Õ (r, r Õ ,t)dr Õ + • ext fl s (r, t) Z _ _ \ (1.20)
where the mesoscopic one-body and two-body densities fl s (r, t) and fl ss Õ (r, r Õ ,t) have been defined in equation (1.11) and (1.12). At this stage, no approximation has been introduced. Equations (1.20) are exact and contain the same information and complexity as equation (1.4) and, therefore, as equation (1.1). However, the time evolution of one-body densities fl s (r, t) is linked to the two-body dislocation densities fl ss Õ (r, r Õ ,t).I t is straightforward to realize that the time evolution of these two-body densities are themselves linked to the three-body densities, and so forth. Obviously, we are faced by the classical problem of closure that we meet in statistical physics when we try to replace a set of discrete degrees of freedom by a set of continuous densities.

As explained above, the next step is to solve the closure problem. This of course requires the introduction of some approximations. One way to do that is to analyse and possibly approximate the two-body correlations, defined by:

d ss Õ (r, r Õ ,t)= fl ss Õ (r, r Õ ,t) fl s (r, t) fl s Õ (r Õ ,t) ≠ 1 (1.21)
Using equation (1.21), the kinetic equation (1.20) becomes: .22) where the local stresses • s sc (r, t) and • s corr (r, t) are defined by:

≠ ˆt fl s (r, t)=sM b • ˆr C fl s (r, t) {• sc (r, t)+• s corr (r, t)+• ext } D . ( 1
• sc (r, t)= ÿ s Õ s Õ ⁄ r Õ " =r • ind (r ≠ r Õ ) fl s Õ (r Õ ,t) dr Õ (1.23) and • s corr (r, t)= ÿ s Õ s Õ ⁄ r Õ " =r • ind (r ≠ r Õ ) d ss Õ (r, r Õ ,t) fl s Õ (r Õ ,t) dr Õ . (1.24)

Mean field stress

≠ ˆt fl s (r, t)=sM b • ˆr C fl s (r, t) {• sc (r, t)+• ext } D . (1.25)
The local stress exerted on the dislocations of sign s does not depend on s and is simply the sum of the external stress • ext and the stress • sc (r, t) generated by all the one-body densities and defined in equation (1.23):

• sc (r, t)= ⁄ • ind (r ≠ r Õ ) ÿ s Õ s Õ fl s Õ (r Õ ,t)dr Õ = ⁄ • ind (r ≠ r Õ )Ÿ(r Õ ,t)dr Õ (1.26)
where we introduced the polar or GND (Geometrically Necessary Dislocation) density:

Ÿ(r, t)= ÿ s Õ s Õ fl s Õ (r, t). (1.27)
As • sc (r, t) does not incorporate any correlation effects, it may be called a mean field stress or, as it closes the theory, a self-consistent stress [START_REF] Groma | Link between the microscopic and mesoscopic length-scale description of the collective behavior of dislocations[END_REF].

Correlation-induced local stresses

We want now to go beyond the mean field approximation and incorporate the correlations. In other words, the correlation stress • s corr (r, t) (equation (1.24)) will be now taken into account. These correlations should be approximated in order to close the theory. We need first to discuss the time and spatial variations of the correlation functions d ss Õ (r, r Õ ,t).I t has already been observed [START_REF] Zaiser | Statistical dynamics of dislocation systems: The influence of dislocation-dislocation correlations[END_REF][START_REF] Groma | Spatial correlations and higher-order gradient terms in a continuum description of dislocation dynamics[END_REF] that the correlation length of d ss Õ (r, r Õ ,t) is finite and of the order of a few average dislocation spacings. Consequently, if the width of the convolution window is sufficiently larger than the mean dislocation spacing, the correlations d ss Õ (r, r Õ ,t), for a fixed point r and as a function of r Õ , decrease to zero before the one-body densities fl s Õ (r Õ ) vary significantly. Therefore, within the domain around point r where they are non-zero d ss Õ (r, r Õ ,t) may be considered as a function of (r≠r Õ ) and of the local one-body densities fl s (r, t):

d ss Õ (r, r Õ ,t) ƒ d ss Õ (r ≠ r Õ , {fl s },t) (1.28)
where the notation {fl s } refers to {fl s ,s= ±}. Now, we comment on the time dependence of the correlations. We recall that the coarse-graining procedure introduced above (see equations (1.6) and (1.9)) involves a time convolution. A width T (L) for the time window must be selected. Due to their short-range nature in space, we argue that, provided L is large enough, the time needed by the correlations d ss Õ (r, r Õ ,t) to reach a stationary state is much smaller than the characteristic time of the time evolution of the one-body densities. Consequently, if the width T (L) of the time convolution is large enough 1 , the coarse-grained correlations reach a stationary state which is dependent on the local one-body densities only: the explicit time dependence in d ss Õ (r, r Õ ,t) disappears and shows up only implicitly through the time dependence of the one-body densities fl s (r, t). In short, equation (1.28) becomes:

d ss Õ (r, r Õ ,t) ƒ d ss Õ (r ≠ r Õ , {fl s }) (1.29)
Now, using again the short-range nature of the correlations discussed above, we note that fl s Õ (r Õ ,t) in equation (1.24) may be expanded to 1 st -order around r.T h e l o c a l s t r e s s d e fi n e d i n equation (1.24) is then split into two terms:

• s corr (r, t)=≠• s b (r, t) ≠ • s f (r, t) (1.30) with • s f (r, t)=≠ ÿ s Õ s Õ fl s Õ (r, t) ⁄ r Õ " =r • ind (r ≠ r Õ ) d ss Õ (r ≠ r Õ , {fl s }) dr Õ (1.31)
and

• s b (r, t)=≠ ÿ s Õ s Õ ˆfls Õ (r, t) ˆr ⁄ r Õ " =r (r Õ ≠ r) • ind (r ≠ r Õ ) d ss Õ (r ≠ r Õ , {fl s }) dr Õ . (1.32)
At this stage, the coarse-grained kinetic equations given in equations (1.22) read: We first need to identify the symmetry properties of the correlation functions d ss Õ .A c c o r d i n gt o their very definition (equation (1.3) and (1.12)), we obviously have:

≠ ˆt fl s (r, t)=sM b ˆr Ë fl s (r, t) Ó • ext + • sc (r, t) ≠ • s f (r, t) ≠ • s b (r, t) ÔÈ (1.33)
d s Õ s (r Õ ≠ r, t)=d ss Õ (r ≠ r Õ ,t) (1.34)
Next, using the discrete kinetic equation (1.1) and its symmetry properties, we obtain the following properties:

d ss Õ (r ≠ r Õ ,t; • ext )=d ss Õ (r ≠ r Õ ,t; ≠• ext ) (1.35) d ss Õ (x ≠ x Õ ,y≠ y Õ ,t; • ext )=d ss Õ (x Õ ≠ x, y ≠ y Õ ,t; ≠• ext ) (1.36)
where s = ≠s and where the dependence of the correlations on the external stress • ext has been explicitly pointed out. Now, in order to get a physical insight into the local stresses • s f and • s b ,w ed i s c u s st h e i rd ependences on the external stress. They inherite this stress-dependency through the correlations d ss Õ , which of course depend on the local stress. Within the spirit of the present coarse-graining procedure, which inevitably leads to a hierarchy of independent and successive many-body densities, we consider that the stress dependence of the k-body densities is due to the stress generated by the correlations up to order (k≠1). Therefore, the stress dependence of the correlations d ss Õ is due to the sum of the external stress and the mean-field stress • sc (r). We note • (r) this sum: • (r)=• ext +• sc (r).U s i n g the symmetry properties given in equation (1.35) and (1.36), it is straightforward to show that the local stresses • s f and • s b defined in equations (1.31) and (1.32) display the following properties:

• s f (r, t, ≠• (r)) = ≠• s f (r, t, • (r)) (1.37) and • s b (r, t, ≠• (r)) = • s b (r, t, • (r)) (1.38)
where the local stress dependency has been explicitly added. These properties clarify the physical meaning of the local stresses • s f and • s b . • s f (r) changes its sign with the local stress • (r) and, as we will see below, is positive when • (r) is positive, whereas • s b is invariant with respect to a change of sign of • (r). As a consequence, • s f , which always opposes the local stress • ext + • sc (see equation (1.33)), plays the rôle of a friction term whereas • s b ,w h i c hb r e a k st h es y m m e t r yw i t hr e s p e c tt oa reversal of the local stress • (r), may generate a Bauschinger effect and a translation of the elastic domain. Therefore, • s b (r) plays the rôle of a backstress.

Broken symmetry in the kinetics of the coarse-grained signed dislocation densities

It is important to notice that, according to equation (1.33), the local stress fields experienced respectively by the positive and negative dislocation densities are different: the correlation-induced stress components • s f and • s b depend on the sign s. In other words, the symmetry that exists at the discrete scale (positive and negative discrete dislocations at the same point r have opposite velocities) is broken at mesoscale: the velocities of positive and negative dislocation densities are not simply of opposite sign. This broken symmetry is the direct consequence of a mesoscale description and its associated coarse-graining procedure: the averaging process required to build a continuous description generates kinetic equations for one-body densities that inevitably incorporate two-body correlations which, in all generality, break the lower-scale symmetry.

In order to be more specific, we analyse explicitly the friction stresses • + f and • ≠ f experienced by the positive and negative dislocation densities, respectively. According to equations (1.31), we have:

• + f (r, t)=≠fl + (r, t) ⁄ • ind (r ≠ r Õ ) d ++ (r ≠ r Õ , {fl s (t)}) dr Õ + fl ≠ (r, t) ⁄ • ind (r ≠ r Õ ) d +≠ (r ≠ r Õ , {fl s (t)}) dr Õ (1.39)
and

• ≠ f (r, t)=≠fl + (r, t) ⁄ • ind (r ≠ r Õ ) d ≠+ (r ≠ r Õ , {fl s (t)}) dr Õ + fl ≠ (r, t) ⁄ • ind (r ≠ r Õ ) d ≠≠ (r ≠ r Õ , {fl s (t)}) dr Õ . (1.40)
Using the symmetry property given in equation (1.34), it is easy to show that the terms that depend on d ++ and d ≠≠ are equal to zero. Therefore, the previous equations reduce to:

• + f (r, t)=fl ≠ (r, t) ⁄ • ind (r ≠ r Õ ) d +≠ (r ≠ r Õ , {fl s (t)}) dr Õ (1.41) and • ≠ f (r, t)=≠fl + (r, t) ⁄ • ind (r ≠ r Õ ) d ≠+ (r ≠ r Õ , {fl s (t)}) dr Õ . (1.42)
Again, using the symmetry properties of equation (1.34), it is easy to show that the integrals in equations (1.41) and (1.42) differ only by their sign. Thus, we have:

• + f (r, t)=fl ≠ (r, t)A(r, t) (1.43) and • ≠ f (r, t)=fl + (r, t)A(r, t) (1.44) with A(r, t)= ⁄ • ind (r ≠ r Õ ) d +≠ (r ≠ r Õ ) dr Õ . (1.45)
Thus, when the signed densities fl + (r, t) and fl ≠ (r, t) are different, which is the generic situation, the friction stresses • + f and • ≠ f are different, which is sufficient to break the symmetry between the velocities of the positive and negative dislocation densities. To better understand this broken symmetry in physical terms, we note that fl ≠ (r, t)d +≠ (r ≠r Õ , {fl s (t)}) may be interpreted as the excess (with respect to the uncorrelated state) of negative dislocations in the surrounding of a positive dislocation that sits at point r. Equation (1.41) tells us that this excess of negative dislocations at r is at the origin of the friction stress • + f experienced by a positive dislocation. There is of course no reason for this excess of negative dislocations around a positive dislocation to be exactly the opposite of the excess of positive dislocations around a negative one. Therefore, the friction stresses • + f and • ≠ f ought to be different. 2 Now, to better visualize this broken symmetry in the signed kinetic equations, we introduce the half sums and half differences of the friction and back stresses:

• f (r)=(• + f (r)+• ≠ f (r))/2 •f (r)=(• + f (r) ≠ • ≠ f (r))/2 • b (r)=(• + b (r)+• ≠ b (r))/2 •b (r)=(• + b (r) ≠ • ≠ b (r))/2 (1.46)
Using equation (1.31) and (1.32), we see that these stresses are linked to the correlations d ss Õ (r ≠ r Õ ) as follows:

• f (r)= 1 2 fl(r) ⁄ • ind (r ≠ r Õ ) d +≠ (r ≠ r Õ ) dr Õ , (1.47) • b (r)=≠ 1 4 ˆfl ˆr ⁄ (r Õ ≠ r) • ind (r ≠ r Õ ) Ó d ++ (r Õ ≠ r) ≠ d ≠≠ (r Õ ≠ r) Ô dr Õ ≠ 1 4 ˆŸ ˆr ⁄ (r Õ ≠ r) • ind (r ≠ r Õ ) Ó d ++ (r Õ ≠ r)+d ≠≠ (r Õ ≠ r)+d ≠+ (r Õ ≠ r)+d +≠ (r Õ ≠ r) Ô dr Õ , •f (r)=≠ 1 2 Ÿ(r) ⁄ • ind (r ≠ r Õ ) d +≠ (r ≠ r Õ )dr Õ , (1.48) •b (r)=≠ 1 4 ˆŸ ˆr ⁄ (r Õ ≠ r) • ind (r ≠ r Õ ) Ó d ++ (r Õ ≠ r) ≠ d ≠≠ (r Õ ≠ r) Ô dr Õ ≠ 1 4 ˆfl ˆr ⁄ (r Õ ≠ r) • ind (r ≠ r Õ ) Ó d ++ (r Õ ≠ r)+d ≠≠ (r Õ ≠ r) ≠ d +≠ (r Õ ≠ r) ≠ d ≠+ (r Õ ≠ r) Ô dr Õ
where Ÿ(r, t) is the GND density defined in equation (1.27) and fl(r, t) the total dislocation density:

fl(r, t)= ÿ s fl s (r, t). (1.49)
By definition, • f (r) and • b (r) are the components of the friction and back stresses experienced by a dislocation independently of its sign, whereas •f (r) and •b (r) are their sign-dependent counterparts.

2. In fact, this broken symmetry could already be pointed out earlier when we wrote the coarse-grained kinetic equations in the form of equation (1.22), where the dislocation-induced stress was split into the mean field stress •sc(r) and the correlation-induced stress • s corr (r, t) defined in equation (1.24). Using the fact that the stress function • ind (r) is odd (see equation (1.2)), it is easy to realize that • s corr (r, t) would be independent of s if and only if the correlations d ss Õ are such that the product fl s Õ (r Õ ,t)d ss Õ (r, r Õ ,t) is equal to the opposite of fl sÕ (r ÕÕ ,t)d ss Õ (r, r ÕÕ ,t),w h e r er ÕÕ and r Õ are symmetric points with respect to r. There is of course absolutely no reason for this to be fulfilled, even if, due to the short range nature of the correlations, point r ÕÕ and r Õ may be restricted to be very close to each other.

Using these stresses, equations (1.33) become:

≠ ˆfl+ (r, t) ˆt = M b ˆr Ë fl + (r) {• ext + • sc (r) ≠ • f (r) ≠ • b (r) ≠ •f (r) ≠ •b (r)} È (1.50)
and

≠ ˆfl≠ (r, t) ˆt = ≠M b ˆr # fl ≠ (r) {• ext + • sc (r) ≠ • f (r) ≠ • b (r)+• f (r)+• b (r)} $ .
(1.51)

Similar equations have already been proposed [START_REF] Zaiser | Statistical dynamics of dislocation systems: The influence of dislocation-dislocation correlations[END_REF][START_REF] Groma | Spatial correlations and higher-order gradient terms in a continuum description of dislocation dynamics[END_REF] 3 , but without the "symmetry-breaking" stresses •f (r) and •b (r) and with a sign-independent backstress • b (r) limited to the term that depends on the gradient of the polar (GND) density Ÿ(r),i . e . t ot h e2 nd term in the right hand side of equation (1.48).

Conclusion

In this chapter, we have clarified the mathematical procedure needed to coarse-grain the dislocation dynamics from the discrete to the continuum. In particular, we have emphasized that the coarse-graining procedure requires a space and time convolution, supplemented by an average on a statistical ensemble. We also argued that, if the width L of the spatial correlation and the width T (L) of the associated time convolution are both large enough, the mesoscopic two-body correlations may be considered locally invariant by translation and stationary at the scale of the characteristic evolution time of the one-body densities. In other words, we may write d ss Õ (r, r Õ ,t) ƒ d ss Õ (r≠r Õ , {fl s (r, t)}). We have explained that the coarse-graining procedure generates correlation-induced stresses • s f and • s b that have specific physical interpretations. The stress • s f always opposes the local stress • (sum of the applied stress and the mean-field stress) and is anti-symmetric with respect to • ; therefore, • s f is a friction stress. The stress • s b is invariant with respect to a change of sign of • . Consequently, it generates a translation of the elastic domain: it is a backstress.

We have also shown that the sign-independent backstress • b ,w h i c hi su s u a l l yl i m i t e dt oat e r m that depends on the gradient of the GND density, contains also a term that depends on the gradient of the total density.

Finally, we argued that the friction and back stresses • s f and • s b , which are sign-dependent, break the symmetry of the kinetic equations: positive and negative dislocation densities do not experience the same local stress and, therefore, display velocities which are not strictly opposite.

The next chapter is devoted to a numerical analysis of the friction and back stresses. In particular, we investigate the rôle of the width L of the convolution window used to define the coarse-graining procedure and its impact on the correlation-induced stresses.

Chapter 2 NUMERICAL ESTIMATION OF THE CORRELATION-INDUCED LOCAL STRESSES WITH 2D DDD SIMULATIONS

Résumé du chapitre

Ce chapitre se consacre à l'estimation des termes de correlation introduits dans le chapitre 1 qui sont nécessaires pour les équations de transport des densités de dislocations. Ces termes, le terme de friction et le terme de 'backstress', dépendent de cartes de corrélations qui peuvent être évaluées dans une situation stationnaire. Par conséquent, il est possible de quantifier ces corrélations. Ces cartes portent l'information sur la probabilité de positionnement relatif entre deux dislocations discrètes. Ainsi, ces cartes sont estimées dans ce chapitre grâce à des statistiques sur un grand nombre de simulations DDD. Une partie préliminaire du chapitre détaille l'implémentation des équations cinétiques des dislocations discrètes en 2D, abordant plusieurs problématiques numériques comme le traitement des dipôles, l'impact d'un système pseudo-périodique ou encore l'analyse des temps caractéristiques des simulations. Il s'est avéré important de s'affranchir de ces artéfacts afin de mesurer avec précision ces cartes de corrélations et leur dépendance à la contrainte appliquée et à la taille de la simulation.

Une deuxième partie illustre les différentes cartes obtenues et propose une interprétation. Enfin, dans une dernière partie, ces cartes sont utilisées pour quantifier les termes de corrélation dans le cadre d'un modèle à deux dimensions.

Introduction

As described in chapter 1, the transport equations of a mesoscale dislocation density theory contain correlation induced stresses, in particular the friction term and the backstress term. These terms depend on the correlation functions d ss Õ , which must be computed through a coarse-graining procedure which has been developed in chapter 1. As explained in section 1.2.2, if the width L of the spatial convolution window is large enough, d ss Õ (r, r Õ ,t) may be considered as a function of (r≠r Õ ) and of the local densities fl s (r, t). The aim of this chapter is to recourse to 2D Discrete Dislocation Dynamics (DDD) to compute numerically these correlation maps and, then, to analyse the flow and back stresses as a function of the external stress.

The computation of these maps requires an ensemble average that is obtained by averaging the result of numerous discrete dislocation simulations. Due to the strong short-range and long-range interactions between dislocations, the numerical treatment of dislocation dipoles and of boundary conditions is usually subtle in Discrete Dislocations Dynamic (DDD) simulations. In order to determine rigorously these correlation maps and their dependences, such as on the applied stress, it is crucial to minimize and control the numerical artefacts. The aim of this chapter is to emphasize important difficulties and results already present in the 2D geometry considered in chapter 1. We consider only straight edge dislocations of one slip-system which involves only two types s of dislocation, called + and ≠. From now onwards, the axis ę1 of the basis is aligned with the slip direction of all dislocations. We characterize a positive dislocation with a Burgers vector toward +ę 1 and a negative dislocation for a Burgers vector toward ≠ę 1 . Moreover, we assume that the geometrically necessary dislocations (GND) are negligible compared to the total number of dislocations. The first part of this chapter details the implementation of equations for DDD in 2D. A second part illustrates the obtained correlation maps and proposes interpretations. Finally, using these maps, an analysis of all dependences of correlation terms is performed.

Preliminary requirement: 2D simulations of discrete dislocation dynamics

In this preliminary part we discuss the numerical treatment of the discrete dislocation dynamics used in 2D for the estimation of correlation maps. One of the main difficulties is to reproduce correctly the behavior of close dislocations, such as dipole of dislocations. In fact, these dipoles can artificially oscillate due to the time discretization. In addition, the time discretization may help a dislocation assembly to overcome a high but narrow energy barrier. In most of discrete dislocation studies, these artefacts have no important consequences but in the present investigation they strongly affect the correlation maps. The first section of this part addresses the numerical discretization of equations of DDD. The second section puts on ground the conditions of simulations in a pseudoperiodic system and proposes some tests to calibrate this time discretization. In a third section, we analyze the collective behavior of an assembly of dislocations.

Assumptions for numerical treatment

The kinetics of a distribution of dislocations is governed by the equation (1.1). We write this equation for our simplified case:

ˆx s i (t) ˆt = Mbs • Q a N ÿ j" =i s j • • + ind 1 x s i (t) ≠ x s Õ j (t),y s i (t) ≠ y s Õ j (t) 2 + • ext R b (2.1)
where • ext is the external stress resolved in the slip system and • + ind (x, y) the shear stress at position (x, y) generated by a positive dislocation located at the origin which, within isotropic elasticity, reads as:

• + ind (x, y)= µb 2fi(1 ≠ ‹) • x(x 2 ≠ y 2 ) (x 2 + y 2 ) 2 (2.2)
where µ is the shear modulus and ‹ the Poisson's ratio. Along any direction, this function behaves as 1/r where r is the distance between two dislocations. This induces very strong stresses between close dislocations and a weak stress even for far dislocations. The first step is to settle the simulations in order to represent correctly the mesoscopic quantities defined by the coarse-graining process. In fact, the statistic has to be performed on many simulations with exactly the same densities fl s .A good way to guaranty this density is to use a periodic distribution with a period corresponding to L. This periodicity could bring numerical artefacts, discussed at the end of the chapter. Now, we link the density of a distribution to the number of dislocations included in one period. In 2D, this period has the size L◊L and this link is fl s = N s /L 2 where N s is the number of dislocations of the type s. This density is controlled by two parameters, N s and L. In order to handle only one parameter, we proceed to an adimension of equation (2.2). All distances are normalized by the characteristic length of the problem which is the average distance between dislocations. Thus, all distances x are adimensionalized by x = x Ô fl. For the sake of clarity, all times t are also adimensionalized by t = tM µb 2 fl/(2fi(1≠‹)). Finally, using a first order scheme for the time discretization, the equation

(2.2) becomes xs i ( t + ∆ t) ≠ xs i ( t) ∆ t = s Q a N ÿ j" =i s j • (x s i ≠ xs Õ j )((x s i ≠ xs Õ j ) 2 ≠ (ỹ s i ≠ ỹs Õ j ) 2 ) ((x s i ≠ xs Õ j ) 2 +(ỹ s i ≠ ỹs Õ j ) 2 ) 2 + 2fi(1 ≠ ‹)• µb Ô fl R b (2.3)
where the stress induced is written in terms of xs i and xs j . The adimensional equation has now two relevant dependences. From then on, we mainly discuss with adimensional quantities: x, t, • . As detailed, this equation is implemented with an explicit discretized time differential with the adimensional time step ∆ t. On the contrary, we keep the exact spatial position of each dislocation to describe short-range correlations with high accuracy. pattern contains small walls of identical type of dislocations and also local coupling between opposite dislocations. In case of two isolated opposite dislocations, it is named a dipole.

Dipole behavior: time discretization and pseudo-periodicity

An elementary system of dislocations is the dipole. Due to its simplicity and to its importance in dislocation pattern, it is very useful for the investigation of the calibration of the time discretization. Thus, we now consider two close dislocations of opposite Burgers vector inside the main box. The dipole size is chosen much smaller than the box size L to investigate the case of an isolated dipole. This dipole size is defined with the difference between the y-coordinates of the two dislocations, called "y. Due to the dipolar effect, the stress generated by a dislocation dipole decays as 1/r 2 where r is the distance to the dipole. Consequently, if the dipole size is much smaller than the box length the periodic copies of the dipole will create a very weak stress inside the main box. In this case, the dipole can be assimilated to an isolated. The equilibrium position of an isolated dipole under stress is obtained when the applied resolved stress is equal to the resolved stress generated by the other dislocation:

2fi(1 ≠ ‹) µ "y b • = x(1 ≠ x 2 ) (x 2 + 1) 2 "y. (2.4)
The graphic solution of this equation is given in figure 2.2. For an applied stress below the critical value

• crit = b "y µ 2fi(1≠‹) 1 4
, the equilibrium distance along ę1 between the opposite dislocations is x eq .F o ra na p p l i e ds t r e s sj u s tb e l o w• crit , the equilibrium distance along x is x min =( Ô 2 ≠ 1)"y corresponding to the angle -= fi/8. The maximum time step that can be obtained by imposing that during a time step ∆ t, the displacement of the dislocation ∆x is much lower than the width of the stress peak around x min /"y in figure 2.2. In adimensional units, this width is around "y Ô fl.U s i n g

• crit as a maximum value for the resolved stress acting on the dislocation, we obtain the condition

∆ t 6 Cste • " ỹ2 (2.5)
where Cste is a numerical constant. After testing many different configurations, we found that an acceptable time step is ∆ t ƒ " ỹ2 . For an applied stress above • crit ,t h e r ei sn om o r ee q u i l i b r i u m position. The numerical test of this behavior is presented in figure 2.3 where the size of a dislocation dipole under stress in our DDD simulations with N rep =2 and "y/L =3.9 • 10 ≠2 is compared to the result of an isolated dipole. After relaxation, we measure the relative position of the blue dislocation through the -angle. The grey dislocation illustrates the equilibrium position in case of no external stress which gives an angle of fi/4 radian. When adding a positive external stress, the stationary position of the blue dislocation is shifted to the left. Figure 2.3 shows the evolution of the -angle as a function of the external stress. The comparison between the results of our periodic simulation and the theory of an isolated dipole shows a very good agreement and the small difference is due to the stress created by the copies of the dipole. After this static investigation, we study the dynamics of a dipole to calibrate the time discretization. The dipole is a good candidate because the formation of a dipole seems to be the fastest phenomenon inside a collective behavior of dislocations. In order to confirm this aspect, we will also carry out some tests on collective dislocations, discussed after. The present test is performed for an external stress taken just below the previous threshold delimiting the equilibrium states. The goal is to find the largest ∆ t which reproduces correctly the stop of the blue dislocation coming from the right of The right side illustrate the system. The left side compare the approximation of a dipole in our pseudo-periodic system compared to the theory. The dipole size is "y =3.9

• 10 ≠2 • L.
the red one. In fact, if we take a too large ∆ t, the dislocation will go over this equilibrium position and will continue gliding toward the left. Figure 2.4 shows the strain rate 'p of a dipole function of the time and for simulations with different ∆ t. If the simulation approaches the equilibrium state, the strain rate will tend to zero. For the largest time step, ∆ t = 10 ≠3 , the negative dislocation goes straight on without any influence of the red dislocation. When we decrease the time discretization to 10 ≠4 or 10 ≠5 , we obtain a wrong behavior with oscillations of the blue dislocation around the red dislocation, alternatively on its right and its left sides. Finally, for ∆ t = 10 ≠6 , the negative dislocation stops on its equilibrium position. After testing many different sizes of dipole, we found that an acceptable time step is ∆ t = Cste • " ỹ2 where the constant is close to 1.

d p d t t ∆ t 0 5.0e -4 1.0e -3 1.5e -3 2.0e -3 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 1e -3 1e -4 1e -5 1e -6
Another time discretization method was tested using an adaptative time step. The time step was reduced if the estimation of the energy increment was positive. We found that in a system containing more than ten dislocations, this method was unable to avoid oscillations of dislocation positions. In addition, the problem of using a too high time step discussed above was not solved by adding the energy criterion. The behavior of the dipoles inside the microstructure were therefore not always correct. As conclusion, we have not pursued further this route and a constant time step has beed used in the following simulations. The next investigation on dipoles is the evaluation of the impact due to the pseudo-periodicity. In fact, the passage of a dislocation through an interface delimiting the main box is traduced by a disappearance of a dislocation at one extremity of the pseudo-periodic system and a creation of another dislocation at the opposite extremity. That creates a stress jump received by all dislocations inside the main box and alters the energy of the system. Moreover, this phenomenon creates artificial correlations at the range of L, seen on correlation maps shown after. The present test, illustrated by figure 2.5, quantifies this artificial stress jump for two dislocations inside the main box with opposite Burgers vector. N rep is the number of layers of boxes surrounding the main box. This between the two main dislocations in the main box in order to emphasize the stresses created by the copied dislocations on the main dislocations. We still consider a blue dislocation as fixed and a red dislocation slipping along ę1 with the coordinate x. The right hand side of the figure is a sketch of this test with the active area delimited by the cut-off. Dislocations inside the circle create a stress contribution on the fixed dislocation whereas dislocations outside do not contribute. The left hand side of the figure shows the evolution of the total stress received by the fixed dislocation while the red ones move. The test is performed for two different cut-off, 1 and 10 times the box size and the number of layer N rep is taken bigger than the cut-off.F o rt h es m a l l e s tc u t -o ff,w efi n dau n i q u e stress jump which is of the order of the jump without cut-off revealed on the previous test. When the cut-off increases, this jump is subdivided to small jumps. In fact, each jump is equivalent to the stress created by one dislocation positioned at the distance r co . Consequently, the jump amplitude behaves as the inverse of r co . However, even small jumps alter significantly the energetic landscape of the system. A good solution is to add a smooth transition rather than a sharp cut-off.I t i s traduced by multiplying

• s Õ ind ( r ) with 0.5[1≠tan ≠1 ((| r |≠r co )/A s )]
where A s is the parameter of the size of the smooth. Thus, a dislocation, going out of the delimited area, will slowly disappear. The green curve illustrates the same simulation of the blue one with r co = 10 but with a smoothed cut-off. Consequently, we obtain a continuum stress received by each dislocation. We observe that this additional smooth allows only few layers to have a correct energetic landscape. In our simulations we have used the typical value of A s =5• 10 ≠2 . This value has been selected by considering its influence on the value of • f as it is discussed below.

Collective behavior

Now that the behavior of a dipole is correctly reproduced, we analyse the collective behavior of an assembly of dislocations. First of all, we test whether the selection of the time step based on the behavior of a dipole is also relevant for an assembly of dislocations. Secondly, we characterize the two stages of simulations: the transient stage and the stationary stage. From this section, we restrain the study to the case N + = N ≠ . We analyse an history of a distribution through the strain rate of plasticity as a function of time. Figure 2.7 shows the plastic strain rate evolution for different time steps ∆ t. This is an evolution of an initial random distribution containing N = 100 dislocations during its relaxation with no external load. As expected, we converge toward the same history when ∆ t tends to zero. In fact, if we con- sider the history with ∆ t = 10 ≠7 as a perfect referential history, the relative error of the history with ∆ t = 10 ≠6 compared to this reference is smaller than 1%. However, for larger time discretizations, this error is non-negligible.

d p d t t 3e -4 -3e -4 -6e -4 1e -2 2e -2 3e -2 4e -2 0 0 ∆ t
In an assembly of dislocations, we expect that the movement of dislocations will be correctly reproduced by equation (2.3) if the criterion ∆t Ø " ỹ2 is fulfilled for all dipoles in the assembly. In other words, the time step is limited by the square of the size along ę2 of the smallest dipole " ỹmin .I n order to confirm this point, we measure the smallest " ỹ inside different distributions. Simultaneously, we measure also the correct ∆ t giving a negligible relative error below than 1%. We conclude that the law ∆ t ƒ " ỹ2 , established in section 2.2.2, is reliable for all collective behaviors. Now, we analyse the collective behavior as a function of the applied stress • Prior to the investigation of the correlation maps, we first analyse the overall collective behavior in term of the plastic strain rate. The analysis relies on an average on a statistical ensemble of disordered initial configurations. Each distribution of this ensemble is randomly initialized but with a restriction on the relative distance " ỹ between two dislocations. More precisely, we impose the minimal value " ỹmin = 10 ≠2 and we select the time step ∆ t = " ỹ2 min . This restriction prevents the occurence of distribution with very small dipoles which would require a too high numerical cost. After this initialization, a quick relaxation of each distribution is performed with no external load.

We now discuss our simulation results together with the results obtained by [START_REF] Miguel | Dislocation Jamming and Andrade Creep[END_REF]. The comparison can only be qualitative because the 2D DDD model in [START_REF] Miguel | Dislocation Jamming and Andrade Creep[END_REF] slightly differs from our model, in particular because annihilation of dislocation can occur. Following [START_REF] Miguel | Dislocation Jamming and Andrade Creep[END_REF] we present in figure 2.8 the evolution of log( 'p ) as function of log(t). We observe qualitatively the same behavior in each investigation. For each value of • ,therea ret w o specific stages: the transient stage and the stationary stage. The stationary stage corresponds either to a static state, if the stress is low enough, or to a steady state, associated to a permanent plastic flow, if the stress is large enough. During the transient stage, the strain rate decreases as function of time and follows the Andrade law [START_REF] Cottrell | Andrade creep[END_REF][START_REF] Cottrell | Strain hardening in Andrade creep[END_REF][START_REF] Cottrell | Logarithmic and Andrade creep[END_REF][START_REF] Cottrell | A microscopic theory of Andrade creep[END_REF] which has been experimentally observed [START_REF] Trouton | On the Viscous Flow in Metals, and Allied Phenomena[END_REF]. This law predicts that the decrease of the strain rate follows a power law with a coefficient ≠2/3. This law is perfectly reproduced by both investigations. A power law fit is illustrated on each side by black lines. The second state is the stationary state which is characterized by a time invariance of the strain rate. We observe two distinctive stationary stages separated by a stress threshold • and estimated around 0.05 for the first study and 0.3 for our study. Below this threshold, the simulations tend toward a frozen state associated with a strain that decreases to zero. Above this threshold, the simulations tend toward a steady state with a non-zero strain-rate constant over time. This threshold defines the activation of the plastic flow which should correspond to the Taylor threshold (I.2). This threshold is well characterized on the left study but can only be estimated between 0.2 and 0.4 in our study. We also observe a difference on the amplitude of the strain rate between the two investigations also due to the incorporated physics. In fact, the in left investigation, dipoles smaller than a characteristic length are annihilated. This avoids the oscillations discussed in section 2.2.2. This restricts the simulations to low enough densities a physics present only inside high local dislocation density. In opposition, in the right investigation, we keep all dipoles. Both investigations are relevant but correspond to different ranges of dislocation density. Now we focus on the duration of the transient stage. Figure 2.8 makes it clear that this duration strongly depends on the applied stress, and diverges when the applied stress comes close to the critical stress • and . As a consequence, it is numerically very costly to get information on the stationary regime close to • and .

In the next section, we investigation the correlation maps as a function of the applied stress • .F o r values of • that lead to a steady state, we simulate the flow up to reduced time t = 10.P r a c t ically, this guaranties that the stationary regime will be significantly established. Figure 2.9 shows that it is indeed the case for applied stresses above • =1.68. In fact, the characteristic time of the left curve (• =1.68) is of the order of t =4. For the right curve (• =3.37), it is more difficult to define precisely a characteristic time but it is definitely even smaller than in the previous situation. The last point of our current simulations is the influence of the numerical parameter N ,r e p r esenting the number of dislocations inside the main box of our pseudo-periodic systems. The right side of figure 2.8 precise this parameter for each curve. The curves obtained with • =0.084 are performed for two distinctive numbers: N = 50 and N = 200. We observe that this number has a significant influence on the results This point will be addresses when investigating the friction term in section 2.3.3. We first analyse the dependence of the correlation maps.

Correlation terms analysis using the coarse-graining procedure

In this part we proceed to the estimation of correlation terms. First, we measure correlated position between dislocations on a large quantity of realizations in the framework of the coarsegraining procedure developed. Then, we define the correlation maps. Finally, we use these maps to estimate the friction and the backstress terms. We consider here a situation with fl + = fl ≠ ,i . e . Ÿ = fl + ≠fl ≠ =0.

Implementation of the coarse-graining procedure

As explained in chapter 1, the basic ingredient of our density-based theory is a coarse-graining procedure. This procedure induces a spatial convolution defined by the width L, a time convolu-tion defined by the width T (L) and, finally, an average over a statistical ensemble of initial conditions.

The spatial convolution window is simply a constant window function of size L,t h el i n e a rs i z e of the DDD simulation box. As discussed in section 1.2.1, L should be of the order of the spatial resolution of the continuous model we want to develop and, of course, sufficiently larger than 1/ Ô fl, the average distance between dislocations. Meanwhile, this garanties that L will always be significantly larger than the range of the correlations d ss Õ . As a consequence, the pertinent choice for the time window, that in all generality should depend on L,i st os e l e c tT (L) of the order of the average time needed by the dislocations to glide along the distance L. This guaranties that T (L) is long enough but still smaller than the characteristic time of the evolution of the one-body densities.

Finally, this space and time convolution is supplemented by a statistical average over an ensemble of random initial dislocation configurations, as explained in section 1.2.1.

Estimation of correlation maps

This section presents the procedure for the obtention of correlation maps. For a given set (s, s Õ ), we proceed as follows. For each initial random dislocation configuration, we let the system evolve for a time of the order of the time needed by dislocations to travel through the simulation box. We then select typically the configurations that belong to the last third of the total simulation time and this procedure is repeated for a set of initial configurations. Let N c the number of configurations selected at the end of this overall procedure. For each of these N c configurations, we identify all the dislocation pairs of type (s, s Õ ) and report the vector ų = r ≠˛Õ on a 2D frame. This leads to a discrete maps, noted Γ ss Õ , that consists in N c N s N s Õ points, where N s and N s Õ are the numbers of dislocations of sign s and s Õ ,r e s p e c t i v e l y .

From these discrete maps, we can obtain an estimation of fl ss Õ (r)dr (and thus of d ss Õ (r))b y counting the number of points inside the surface element dr at position r. The obtained maps are presented in figure 2.10 with an averaging resolution 256 ◊ 256. Note that the quality of the maps strongly depend on the chosen resolution. The d ss Õ (r) maps can than be used to compute the correlation terms. However, due to stiff gradients close to r =0, the accuracy of this approach is limited. A better way to proceed is to directly estimate the integrals in the definition of the correlation induced stresses using the discrete set of points Γ ss Õ without explicitly computing the correlation maps d ss Õ (r), i.e. for any function f (r),w eu s e :

⁄ d ss Õ (r)f (r)dr ƒ 1 N c N s N s Õ fl s fl s Õ ÿ roeΓ ss Õ f (r) ≠ ⁄ f (r)dr (2.6)
where fl s is the density of dislocations of sign s inside the L◊L box. Each column in figure 2.10 represents an ensemble of measures for a specific applied stress, from zero to • =0. orientation of the dipoles is close to fi/4 radians in the absence of applied stress (figure 2.10 top row, left). When a stress is applied, the d +≠ correlation maps are no longer symmetric with respect to the y-axis and the angle with the y-axis of the most probable dipole orientation is smaller than fi/4 radians. These results are in agreement with the correlation maps presented in [START_REF] Miguel | Dislocation Jamming and Andrade Creep[END_REF][START_REF] Ispánovity | The probability distribution of internal stresses in externally loaded 2D dislocation systems[END_REF]. The fact that we prevent two dislocations to be closer than " ỹmin along the y-axis leads to a narrow band along the x-axis with a negative value. We have verified that the influence of the value of " ỹmin on the correlation terms • f and • b is negligible if " ỹmin is lower than 10 ≠2 . A second artefact appears specifically on the two maps obtained with • =0 .15µb Ô fl. We observe diffuse and weak bands parallel to the y-axis at the boundaries of the computational domain. This artefact is minimized by increasing the cut-off and the with of the smoothing function introduced in section 2.2.2. In any case, we verified numerically that these diffuse bands have no impact on the values of the correlation terms. Now, we enumerate the symmetries that the correlation maps may display. -First, dislocations are restricted to glide along the x-axis and the dislocation-dislocation interaction • ind (x, y) is even in y. Therefore, inspection of the discrete kinetic equation (2.1) shows that the correlations d ss Õ are symmetric with respect to the x-axis:

d ss Õ (x, y)=d ss Õ (x, ≠y). (2.7)
-Also, as already mentioned (see equation (1.34)), we have:

d ss Õ (x, y)=d s Õ s (≠x, ≠y). (2.8)
-T h e r e f o r e ,w h e r es = s Õ , the maps are centro-symmetric:

d ss (x, y)=d ss (≠x, ≠y).

(2.9)

There are supplemental symmetries in specific cases. If N + and N ≠ are equal, which is the situation considered here, we have:

N + = N ≠ ∆ d ++ (x, y)=d ≠≠ (x, y). (2.10)
Also, in the absence of stress, and because • ind (x, y) is even in x, the correlations are symmetric with respect to the y axis:

• ext =0 ∆ d ss Õ (x, y)=d ss Õ (≠x, y). (2.11)
The results on figure 2.10 confirm that the correlations d ++ and d +≠ display small correlation lengths. The larger correlation length is obtained with • =0.15µb

Ô fl. For higher applied stresses, the correlation length is smaller. In any case, correlations are essentially limited to a close neighborhood of the origin. When we use a "y min small enough, we observe a divergence close to the origin. This behavior is consistent with the 1/r behavior predicted at short range in [START_REF] Miguel | Dislocation Jamming and Andrade Creep[END_REF].

Friction term

In this section we use the previous correlation maps to estimate the friction term given by (see equation (1.47)):

• f (r)= µb 2fi(1 ≠ ‹) fl(r, t) 1 2 ⁄ rÕ " =0 d +≠ ! x Õ ,y Õ , fl, •,L " x Õ (x Õ2 ≠ y Õ2 ) (x Õ2 + y Õ2 ) 2 dx Õ dy Õ (2.12)
where, on top of the obvious external stress and density dependencies, we explicitly point out an L dependence. In order to be consistent with the rescaled correlation maps, we adimensionalize also all distances by 1/ Ô fl. The friction term becomes:

• f (r)= µb 2fi(1 ≠ ‹) Ò fl(r, t) ⁄ rÕ " =0 1 2 d +≠ A xÕ , ỹÕ , 2fi(1 ≠ ‹)• µb Ô fl ,L Ô fl B xÕ (x Õ2 ≠ ỹÕ2 ) (x Õ2 +ỹ Õ2 ) 2 dx Õ dỹ Õ (2.13)
which may be written as:

2fi(1 ≠ ‹)• f (r) µb Ô fl = f A 2fi(1 ≠ ‹)• (r) µb Ô fl ,L Ô fl B (2.14)
where f is the integral estimated numerically by using the previous correlation maps.

Generally speaking, we may expect that the coarse-graining length L will show up in the coarsegrained quantities that result from the coarse-graining procedure. The important point is that we are dealing here with a situation where many length scales may emerge from the complex spatial and dynamical coupling that governs the dislocation dynamics. It is indeed well known that, most often, dislocations self-organized themselves into complex patterns that display length scales much larger than the average dislocation spacing, such as dislocation walls in cyclic loading [73] or even seemingly fractal structures [START_REF] Zaiser | Fractal analysis of deformation-induced dislocation patterns[END_REF] with no characteristic length scale [82]. In such situations, when many different large length scales are physically present, an averaging procedure at a given intermediate length scale will generate a continuous theory which is scale dependent. In the present context, it means that the correlation-induced stresses generated by coarse-graining may definitely display an L-dependence.

Therefore, in order to investigate this important feature, we consider below different values for L. In fact, as the only pertinent quantity is L Ô fl,w ea n a l y s ed i fferent values for

Ô N = L Ô fl,w h e r e
N is the total number of dislocations.

The computations are performed for different applied stresses. The limits of this stress dependence are analytically predictable. In fact, when • tends to zero, there is an axial symmetry on d +≠ with respect to y-axis, as explained in the previous section. Consequently, the function f integrates an odd function, hence equals to zero. Secondly, when • tends to the infinity, the individual interactions,

• s Õ ind , in the DDD equation (2.
3) is negligible compared to the stress • . Thus, the movement of each dislocation depends only on this applied stress. Consequently, two dislocations with opposite movements cannot be correlated contrary to two dislocations of same type. That gives the nullity of the map d +≠ which necessarily implies the nullity of the function f . Figure 2.11 illustrates the numerical results between these two limits and for three different values for the parameter

Ô N = L Ô fl.
These computations are performed with simulations with an artificial smooth at a radius r co =2, introduced in section 2.2.2. This smooth allows the use of a small pseudo-periodic systems. In fact, we observe that between r co =1 and r co =2 the relative error on • f values is smaller than one percent. It is important to note that this smooth is the main parameter of the band artefact seen on figure 2.10 on the second column. In fact, this artefact is strongly limited by this smooth. We take A s equals to 5 • 10 ≠2 , representing the characteristic width of the smooth, which is sufficient to have no impact on correlation term values. The results presented in figure 2.11 are consistent with the predicted behavior when the applied stress tends to zero and to infinity. Indeed, for • =0, we obtain • f =0 and for high stresses, • f decreases to zero. In fact, the friction stress displays two different regimes. The first one is for small applied stresses (• up to approximately 1.5) where the friction term is approximatively linear with a slope close to 1. Therefore, the friction term opposes almost totally the applied stress. This is associated to a quasi-static state where there is no effective flow of dislocations. For higher applied stresses, we observe a permanent flow of dislocation. These behaviors are in agreement with the observations on DDD simulations. Now, we comment of the dependence of the friction stress • f with the parameter L Ô fl. The physical origin of the increase of the number of dipoles with L (at constant dislocation density) is that the coarse-graining procedure involve a time convolution with a width T (L) of the order of the travelling time over the length L. Therefore, the probability that a given dislocation meets another dislocation of opposite sign during the time T (L) increases with L. In brief, the longer L, the higher the number of dipole that have the time to form. However, we note that this physical phenomena may be here disturb by the use of periodic boundary conditions because a dislocation may travel through the simulation box more than once. This undesirable effect may be avoided with a careful numerical monitoring of T (L), which has not been done here. Therefore, the L-dependence observed in figure 2.11, even if it has a true physical origin, is not perfectly quantitative.

2π(1 -ν)τ f µb √ ρ L √ ρ = 10 L √ ρ 14.1 L √ ρ = 20 2π(1 -ν)τ µb √ ρ 0 0 2 4

Backstress term

In this section we proceed to an estimation of the backstress term using correlation maps. In the same way as for the friction term, we adimension all distances by 1/ Ô fl.

We limit the analysis to zero external stress and, again, to a situation with fl + = fl ≠ = fl/2.T h e n , equation (1.35) leads to

• s b (r)= µb 2fi(1 ≠ ‹) 1 fl ˆfl(r, t) ˆx ÿ s Õ ⁄ 1 2 d ss Õ A xÕ , ỹÕ , 2fi(1 ≠ ‹)• µb Ô fl B xÕ2 (x Õ2 ≠ ỹÕ2 ) (x Õ2 +ỹ Õ2 ) 2 dx Õ dỹ Õ . (2.15)
rewritten as

• + b (r)= µb 2fi(1 ≠ ‹) 1 fl ˆfl(r, t) ˆx (C ++ + C +≠ ) . ( 2.16)
and

• ≠ b (r)=≠• + b (r) (2.17)
where the symmetry properties given in equations (1.34), (1.35) and (1.36) have been used. The quantities C ++ and C +≠ are given by

C ++ = 1 2 ⁄ d ++ (x Õ , ỹÕ , • =0,L Ô fl) xÕ2 (x Õ2 ≠ ỹÕ2 ) (x Õ2 +ỹ Õ2 ) 2 dx Õ dỹ Õ (2.18)
and

C +≠ = ≠ 1 2 ⁄ d +≠ (x Õ , ỹÕ , • =0,L Ô fl) xÕ2 (x Õ2 ≠ ỹÕ2 ) (x Õ2 +ỹ Õ2 ) 2 dx Õ dỹ Õ (2.19)
where the symmetry of d ss Õ with respect to the x-axis has been taken into account. For L Ô fl = 20, the numerical estimation gives C ++ ƒ 1.09 and C +≠ ƒ≠0.23 with a relative error smaller than 10 ≠3 .

In view of its definition (see section 1.2.4), the physical meaning of this backstress term is distinctively different from the backstress term proposed in other dislocation density formulations [START_REF] Evers | Non-local crystal plasticity model with intrinsic SSD and GND effects[END_REF][START_REF] Bayley | A comparison of dislocation induced back stress formulations in strain gradient crystal plasticity[END_REF][START_REF] Geers | Second-order crystal plasticity: internal stress effects and cyclic loading[END_REF] which do not consider correlations. In these studies, the "backstress" term can be viewed as originating from a decomposition of the self-consistent stress given by equation (1.23) by splitting the integral into two parts:

• sc (r)= ⁄ r Õ " =r Ÿ(r Õ ,t)• + ind (r ≠ r Õ )dr Õ = ⁄ r Õ " =r |r Õ i ≠r i |6R Ÿ(r Õ ,t)• + ind (r ≠ r Õ )dr Õ + ⁄ r Õ " =r |r Õ i ≠r i |>R Ÿ(r Õ ,t)• + ind (r ≠ r Õ )dr Õ (2.20)
where, in principle, R should be small enough for the gradient of Ÿ to be considered constant at the scale R.E x p a n d i n gŸ(r Õ ,t) in the first integral around the point r, we obtain for the first integral the following estimation:

• local ƒ µb 2fi(1 ≠ ‹) ˆŸ(r, t) ˆx ⁄ r ÕÕ " =0 |r ÕÕ i |6R x 2 (x 2 ≠ y 2 ) (x 2 + y 2 ) 2 dxdy = µb 2fi(1 ≠ ‹) 1 fl ˆŸ(r, t) ˆx • I. (2.21)
Now, we proceed to the same adimensionalization in order to bring out an integral, noted I,w h i c h depends only of the length R = R Ô fl. This integral can be computed analytically and increases monotonously with respect to R.T h eb a c k s t r e s st e r m• local is equivalent to the backstress term proposed by [START_REF] Bayley | A comparison of dislocation induced back stress formulations in strain gradient crystal plasticity[END_REF]. In order to appreciate the amplitude of our correlation-based backstress term • s b , we compare it to this backstress term • local . This comparison is meaningful because of the similar mathematical expressions of both terms. The minimal amplitude for • local is obtained with the minimal usable value of R. The correlations on figure 2.10 impose a lower bound of the order of 10. Then, for R = 10, we obtain I ƒ 85.84. In comparison, the equivalent constant in • s b is equal to 0.86. It is almost hundred times smaller than the coefficient of classical backstress term.

Conclusion

This chapter detailed the numerical procedure, used to apply the coarse-graining procedure presented in chapter 1 to DDD simulations. The first part was dedicated to the identification of a precise and stable numerical scheme adapted to the quantities that control the correlation maps. In fact, we opt for a pseudo-periodic system in order to handle the appropriate dislocation density. That imposes to manage the artefacts linked to this system. Moreover, this part proposes a rigorous investigation of the calibration of the time discretization of DDD equations in order to reproduce correctly the behavior of dipoles. The second part detailed the statistic measures on DDD which is used to build the correlation maps. Then, we used these maps to estimate the correlation terms. We observe in particular on L-dependence of the coarse-grained friction stress. Generally speaking, this length-scale dependence is not surprising, regarding the frequently observed patterns that dislocation dynamics often generate. These patterns generally exhibit characteristic length scales much larger than the average distance between dislocations. Therefore, a coarse-graining procedure based on a length scale L smaller than these configurational length scales will inevitable lead to correlation-induced stresses that are L-dependent. In the present oversimplified situation, where parallel dislocations are limited to a single glide system, the L-dependence has been linked to the dynamical formation of short-range dipoles associated to the spatial and time scales of the coarsegraining procedure.

Chapter 3 SIMULATIONS OF DISLOCATION DENSITY PATTERNING WITH ONE SLIP SYSTEM IN 2D Résumé du chapitre

Ce chapitre propose en première partie un schema d'integration numérique pour des équations de transport de densité du type (1.22), formulées en 2D. Ce nouveau schéma est inspiré du schéma Lax-Wendroff Haute-Résolution ainsi que de construction de schémas conservatifs. De plus, l'équation d'Orowan doit être discrétisée de manière similaire afin d'avoir une cohérence entre le champs plastique calculé et son origine physique, les GND. Ce traitement numérique permet de limiter fortement des artefacts tels que la diffusion numérique, les instabilités numériques ou encore les problèmes de conservation des champs. En deuxième partie, une simulation 2D utilisant le modèle formulé illustre l'évolution de densités de dislocations initialement homogènes sollicitées par un cisaillement du matériau. Il est observé une emergence d'un auto-arrangement des dislocations, ce qui confirme les instabilités prédites sur le comportement du modèle.

Introduction

This chapter treats the numerical implementation of transport equations such as the ones developped in chapter 1. This numerical aspect is not much detailed in the literature on dislocation densities. In the first part, we propose and compare different numerical schemes for these transport equations. The behavior of these schemes is independent on the stress components incorporated into the velocity fields of the dislocation densities. Consequently, for the sake of simplicity, we adopt the simplified case of chapter 2: we consider only one slip system with edge dislocations in an isotropic material. Thus, there are only two types of dislocation density, noted fl + and fl ≠ , depending on the direction of the Burgers vector. Moreover, we restrain the model with a GND density negligible compared to the total dislocation density. Consideration of the GND density is useful because there is a direct link through the Kröner equation, between the GND density and the plastic strain, which in turn can be used to compute the stress field generated by the dislocations. We will fix some inconsistencies related to this link. After an illustration of artefacts attached to classical numerical schemes, we set up a new scheme inspired by the Lax-Wendroff scheme used in fluid mechanics to avoid these artefacts. The second part of this chapter proposes a 2D simulation with our implemented model complemented by a stability analysis of the initial homogeneous state.

Transport equations of density of edge dislocations

We consider the 2D plasticity model developed in chapter 1. For the sake of completeness, we report here the transport equation (1.50) and (1.51):

≠ ˆfl+ (r, t) ˆt = M b ˆr Ë fl + (r) {• (r) ≠ • f (r) ≠ • b (r) ≠ •f (r) ≠ •b (r)} È (3.1) ≠ ˆfl≠ (r, t) ˆt = ≠M b ˆr # fl ≠ (r) {• (r) ≠ • f (r) ≠ • b (r)+• f (r)+• b (r)} $ . ( 3.2)
where

• (r)=• ext + • sc (r)
is the sum of the external stress and the self consistant stress, i.e. the stress generated by the one-body dislocation densities: The stress tensor ‡ sc (r), which is linear with respect to the one-body dislocation density, can be in fact computed through the usual mechanical equilibrium law, in which the source term is the spatial derivative of the plastic strain generated by the dislocation glide. More precisely, if ' P kl (r) is the plastic strain, ‡(r) is the solution of the following mechanical equilibrium problem:

• sc (r)= ⁄ r Õ " =r Ÿ(r Õ ,t)• + ind (r ≠ r Õ )dr Õ . ( 3
ˆ ‡ij (r) ˆrj =0 (3.4) with Hooke's law ‡ ij (r)=⁄ ijkl (' kl (r) ≠ ' P kl (r)) (3.5)
supplemented by the boundary condition

⁄ ijkl 'kl = ‡ ext ij (3.6)
where 'kl is the spatial average of the total stress tensor ' kl (r).T h ep l a s t i cs t r a i n' P (r) that appears in equation (3.5) is linked to the dislocation glide through Orowan's law:

ˆ'p ij ˆt = ÿ s p s ij ˆ"s ˆt with ˆ"s ˆt = fl s (r, t)b s v s (r) (3.7)
where p s ij is the orientation tensor of the glide system s. Finally, with the system of equations coupling Orowan equation (3.7) and the mechanical equilibrium problem given in equations (3.4)-(3.6), we compute easily the mesoscopic local stress • needed by the transport equations (3.1)-(3.2). From a numerical point of view, the resolution of the system of equations (3.4)-(3.6) is easier than the computation of the integrals that appear in the expression (3.3) of the self-consistent stresses. The resolution of the system with the mechanical equilibrium needs specific solvers such as fixed-point FFT-based algorithms which are presented in chapter 4 with new enhanced solvers.

In brief, the kinetic equations that we must integrate consist in the transport equations (3.1) and (3.2) and the Orowan equation (3.7):

Y _ _ _ _ ] _ _ _ _ [ Transport equation: ˆt fl + (r, t)+ ˆx # fl + (r, t)v + (r) $ =0 ˆt fl ≠ (r, t)+ ˆx [fl ≠ (r, t)v ≠ (r)] = 0 Orowan equation: ˆt ' p (r, t)=fl + (r, t)bv + (r) ≠ fl ≠ (r, t)bv ≠ (r) (3.8)
where the notation of the component of the strain tensor is ' p = ' p 12 © ' p xy and the velocities are given by:

v + (r)=Mb •{• ≠ • f ≠ • b ≠ •f ≠ •b } v ≠ (r)=≠Mb •{• ≠ • f ≠ • b +• f +• b } (3.9)
where M = M + = M ≠ is the mobility coefficient and each stress contribution corresponds to a resolved shear stress. Now we discuss each of the stress components that appear in equations (3.9). First, we propose a fit of the friction term • f , based on the estimations performed in section 2.3.3. • f could be naively fitted on the curves displayed in Figure 2.11, but it appears that it is more appropriate to fit (• ≠ • f ) rather than • f . Figure 3.1 illustrates this fit. As already proposed in [START_REF] Miguel | Dislocation Jamming and Andrade Creep[END_REF], we use a power law of the form: We note that the exponentis almost constant, whereas the prefactor B seems to display a variation with L Ô fl. We have therefore decided to perform a new fit of the three curves for which the exponent is fixed to the value 1.71 (Figure 3.1). We observe that, for each of the three situations analysed, the power law fit matches correctly the numerical points, except the first one, which corresponds to the smaller analysed stress value. Therefore, it seems that the dislocation dynamics display two distinct regimes, eventhough the present analysis is not precise enough to identify with precision a "critical" stress that would separate this two regimes. More intensive computations would be necessary to explore this point. Nevertheless, we decided to consider that the power law fit are valid even in the limit • ae 0, which implies that the expression (• ≠ • f ) is always finite when • is finite, which seems a reasonable hypothesis1 . 3), respectively. We restrain the analysis to Ÿ π fl.I nt h a tc a s e sa tz e r oa p p l i e ds t r e s s ,d ++ (r Õ ≠ r) and d ≠≠ (r Õ ≠ r) can be assumed equal. Because d ++ (r Õ ≠r) and d ≠≠ (r Õ ≠r) are practically insensitive to the applied stress [START_REF] Zaiser | Statistical dynamics of dislocation systems: The influence of dislocation-dislocation correlations[END_REF], • b (r) given in equation (1.48) can be approximated by

• ≠ • f (r)=sign(• ) • B A 2fi(1 ≠ ‹) µb Ô fl B -≠1 •|• | -. ( 3
• b (r) ƒ≠ 1 2 ˆŸ ˆr ⁄ (r Õ ≠ r)• ind (r ≠ r Õ ){d ++ (r Õ ≠ r)+d +≠ (r Õ ≠ r)}dr Õ (3.11)
where the correlation maps d ss are computed at zero applied stress. Using dimensionless quantities and symmetry properties given in equations (1.34-1.36), we obtain

• b (r) ƒ (C ++ ≠ C +≠ ) µb 2fi(1 ≠ ‹) 1 fl ˆŸ ˆr (3.12)
where C ++ and C +≠ are dimensionless quantities defined be equations (2.18) and (2.19) respectively.

Using the same analysis and the same approximations for the •b term, we get

•b (r) ƒ (C ++ + C +≠ ) µb 2fi(1 ≠ ‹) 1 fl ˆfl ˆr . (3.13)
Now, we discuss, again within the approximation Ÿ π fl, the influence of the friction stress component •f , whose definition is given in equation (1.3). According to equations (3.9), it has to be compared to • ≠ • f where • = • ext + • sc . Using the power law analysis presented in equations (3.10) and the fact that according to equations (1.47) and (1.3), we have

•f = ≠ Ÿ fl • f ,w ec a ns h o w that, in the limit • ∫ Ÿ fl µb Ô fl 2fi(1≠‹) , •f is negligible with respect to • ≠ • f .
Finally, neglecting •f and using the total and GND densities fl and Ÿ, the transport equations (3.9) can be rewritten:

ˆfl(r, t) ˆt + Mb ˆx [{• ≠ • f ≠ • b } Ÿ(r, t) ≠ •b fl(r, t)] = 0 ˆŸ(r, t) ˆt + Mb ˆx [{• ≠ • f ≠ • b } fl(r, t) ≠ •b Ÿ(r, t)] = 0 (3.14)
In the following, we investigate the numerical treatment of these equations, we proceed to the stability analysis and then to numerical simulations.

Numerical treatment of transport equations

This part investigates the numerical treatment of the equations (3.14). Firstly, the discretization of the transport equations raises specific issues. Some of these are regularly addressed in others contexts, such as numerical diffusion problems in the transport equations of fluid mechanics. Secondly, we need to discretize the Orowan equation. However, this equation is linked to the transport equations. In fact, the strain field, generated by Orowan equation, has to be consistent with the GND density involved in the transport equations. This link is given by the Kröner equation [START_REF] Kröner | Allgemeine Kontinuumstheorie der Versetzungen[END_REF]. Consequently, we will present the particularity of the discretization of the Orowan equation.

Requirement of an advanced scheme discretization

This first section treats only the transport equation of Ÿ but exactly the same treatment would apply to any dislocation density transport equations. The first step is to discretize the time. Implicit discretization appears difficult because the velocity is dependent of dislocation densities. Thus, we choose an explicit discretization leading to

Ÿ(r, t + ∆t) ≠ Ÿ(r, t) ∆t + ˆfl(r, t)v(r, t) ˆx =0. (3.15)
where ∆t is the time step. The second step is the spatial discretization of the divergence operator. Figure 3.2 helps us to describe each possible scheme on a simple grid. For the sake of simplicity, we present here only the direction ę1 . ∆x is the grid step and each discrete field is defined on this grid with the subscript i. For the moment, we consider the velocity v as homogeneous in space.

After managing scheme artefacts, we will discuss some problems for an inhomogeneous velocity. The simplest discretisation scheme of a transport equation is the Euler scheme with a central difference:

-E u l e rs c h e m e : ˆfl(r, t)v ˆx ---- where | i informs about the position where the discrete divergence is computed and | E means the use of the Euler scheme. However, as predicted by stability analysis, a transport equation discretized with an explicit time scheme and a centered space scheme is unconditionally instable [START_REF] Buffat | Modélisation numérique des écoulements compressibles[END_REF]. An alternative scheme used for transport equations is the Up-Wind scheme. This scheme treats the divergence with an off-centered difference depending on the sign of the velocity:

E i = v • fl i+1 ≠ fl i≠1 2∆x (3.16)
-U p -W i n ds c h e m e : ˆfl(r, t)v ˆx - - - - UW i = v • Y _ ] _ [ fl i ≠fl i≠1 ∆x if v>0 fl i+1 ≠fl i ∆x if v<0 (3.17)
where | UW means the use of the Up-Wind scheme. This scheme is stable if we take a time step smaller than the critical value ∆t crit = ∆x/v [START_REF] Buffat | Modélisation numérique des écoulements compressibles[END_REF]. The main drawback of this alternative scheme, well known in fluid mechanics, is the numerical diffusion: strong gradients of a field of density field will decrease. Consequently, the field will become smoother over time. Figure 3.3 illustrates this problem with the time evolution of a square shaped dislocation density in constant and homogeneous velocity field v. The expected evolution is the uniform translation to the left of the square shaped at velocity v.T h ed i ffusion of the blue curves in figure 3.3 shows that the UW numerical scheme correctly reproduces the translation of the dislocation density but the shape of the density peak becomes smoother and broader. This numerical diffusion increases when ∆t decreases.

To reduce the numerical diffusion we turn towards an improved scheme [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF] proposed by Lax-Wendroff,a b b r e v i a t e dLW .T h i s m e t h o d e s t i m a t e s t h e d i ffusion of an Euler scheme. Then, it proposes to add an anti-diffusion term to the discretized equation based on a central difference. However, this central difference coupled with our explicit temporal discretization is still unstable even though oscillations are attenuated by this anti-diffusion term. Moreover, the densities can evolves to negative values when the local gradients are strong. Recent investigations [START_REF] Ma | High-Resolution Simulation of Multidimensional Crystal Growth Industrial & Engineering Chemistry Research[END_REF] propose an e ws c h e m em i x i n gt h eUW scheme and the LW scheme, called Lax-Wendroff High-Resolution, abbreviated by LW -HR. This scheme is an evolution of the equation (3.15) with an anti-diffusion term

Ÿ i (t + ∆t) ≠ Ÿ i (t) ∆t + ˆfl(r, t)v ˆx - - - - UW i + T ad (∆x, ∆t)=0 (3.18)
where the divergence is treated similarly to the UW scheme and T ad is the anti-diffusion term. This new term is also expressed differently depending on the sign of the velocity:

T ad (∆x, ∆t)=v 3 1 ≠ |v|∆t ∆x 4 • Y _ _ ] _ _ [ (fl i+1 ≠fl i )" i ≠(fl i ≠fl i≠1 )" i≠1 2∆x if v>0 (fl i ≠fl i≠1 )" i≠1 ≠(fl i+1 ≠fl i )" i 2∆x if v<0 (3.19)
This term looks like the second order off-centered discretization [START_REF] Ma | High-Resolution Simulation of Multidimensional Crystal Growth Industrial & Engineering Chemistry Research[END_REF] but with a ponderation with a function ", called the limiter. This function is essential because the magnitude of the anti-diffusion term is too large if the shorter wavelengths (compatible with the grid) show up, leading to undesirable oscillations. This function depends on the local gradient through the intermediate function ◊: snaps at t 1 and t 2 =2• t 1 in blue for the UW scheme and in green for the LW -HR scheme. The simulations are performed with a time step equals to 2.6 • 10 ≠2 ∆t crit . As mentioned before, the blue one decreases strongly with a numerical diffusion which is worst when ∆t decreases. Note that in the special case of an homogeneous velocity, this numerical diffusion can be totally removed by ∆t = ∆t crit . But this critical time discretization is unreachable for a realistic case. The green curve looses only the smallest wavelengths but keeps the amplitude for any time discretization. In fact, the LW -HR scheme generates only a small regularisation effect but no diffusion, in opposition to the UW scheme.

" i = |◊ i | + ◊ i 1+|◊ i | with ◊ i = fl i ≠ fl i≠1 fl i+1 ≠ fl i if fl i+1 " = fl i else if " i =1. ( 3
The case of homogeneous velocity is now under control but additional issues appear for the inhomogeneous situations. Before testing an inhomogeneous velocity, we need to precise the index of the discrete field v involved in each product fl • v in the discrete equation. Indeed, equation (3.18) involves several products summarized by fl j • v where j can take the values i or i ± 1.T w o intuitive possibilities appears for the choice of the index of the velocity: v i or v j . For any choices, the scheme does not guaranty the conservation of the density. However, the conservation law is a characteristic of the transport equation (1.22). Figure 3.4 illustrates a simple case where the density is not conserved. This test is similar to the previous test with an initial square wave of density but with a sinusoidal imposed velocity shown in dotted line. The blue curve illustrates the density at the two times t 1 and t 2 , obtained with LW -HR scheme with the choice v i for each product detailed before. We see a strong increase of the size of the wave at the time t 1 which is in contradiction to the conservation law. In order to fix this issue, we need to turn towards an other type of discretization schemes based on the concept of density flux. In a discrete space, this flux characterizes the flow between two neighbors nodes. Consequently, this flux, called f , is defined between each nodes of the main grid. The choice is to keep only one index for all grids involving half-integers for the index of the flux. This staggered grid allows a building of schemes which guaranty the conservation law of the density. These schemes need to verify this structure

Ÿ i (x, y, t + ∆t) ≠ Ÿ i (x, y, t) ∆t + f fl i+ 1 2 ≠ f fl i≠ 1 2 ∆x =0 (3.21)
where f fl i+ 1 2 is the flux of the density fl from fl i+1 to fl i . Thus, we proposes here a new scheme inspired by the UW scheme to avoid oscillations and including anti-diffusion terms inspired by LW -HR.T h e definition of the flux is

f fl i+ 1 2 = v i+ 1 2 • Y _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ [ C fl i + 1 2 " i A 1 ≠ - -v i+ 1 2 - -∆t ∆x B (fl i+1 ≠ fl i ) D if v i+ 1 2 > 0 C fl i+1 + 1 2 " i A 1 ≠ - -v i+ 1 2 - -∆t ∆x B (fl i ≠ fl i+1 ) D if v i+ 1 2 < 0 (3.22)
As wished, this new scheme is analytically equivalent to the LW -HR scheme in case of homogeneous velocity. Figure 3.4 illustrates also the evolution of the square wave with this new scheme in green. This evolution is perfectly conservative. Note that, this scheme imposes the discrete velocity field in the same position of the flux. Therefore, the driving force terms (3.9) need interpolated values of 1/fl on the staggered grid. Moreover, the gradient of fl inside • b is simply the central difference between two close nodes. We perform systematic numerical tests to attempt to predict the stability of our scheme. It seems to have the same stability range than the UW scheme. To conclude the validation of this new numerical scheme, figure 3.6 shows the impact of ∆t. We observe exactly the

Requirement of specific treatment of mechanical fields

This section treats the discretization of the Orowan equation. As explained in the previous section, the local stress • , computed from the mechanical equilibrium (3.4), is defined on the staggered grid illustrated by the blue circles on figure 3.5. Moreover, this mechanical equilibrium needs, in input, the plastic strain ' p . Consequently, the Orowan equation, which gives this plastic strain, has to be computed on this staggered grid. In this way, the simplest discretization is

' p i+ 1 2 (t + ∆t)=' p i+ 1 2 (t)+ 1 2 • ∆t • b • fl i+ 1 2 • v i+ 1 2 (3.23)
where Ÿ i+ 1 2 is the average of Ÿ i and Ÿ i+1 . The behavior of the model (3.8) with this discretization scheme leads to numerical divergences. The problem comes from the inconsistency of two related fields, namely the GND density Ÿ and the plastic strain ' p . In fact, the Kröner equation [START_REF] Kröner | Allgemeine Kontinuumstheorie der Versetzungen[END_REF] relates the plastic strain to its physical origin, the GND. This law is reformulated for our simple case by

-=curl- in our 2D case ≠≠≠≠≠≠≠≠ae -13 = bŸ = ≠2 ˆ'p ˆx (3.24)
where -is the field of Nye tensor and -is the distortion tensor of the lattice. Both tensors are second order tensors. In order to guaranty this consistency, the Orowan equation has to be discretized exactly in the same manner than transport equations. Consequently, the Orowan equation computes the flow of Ÿ with the expression (3.22) and incorporates, curiously, the anti-diffusion term. Then, the discretization is

' p i+ 1 2 (t + ∆t)=' p i+ 1 2 (t)+ 1 2 • ∆t • b • f fl i+ 1 2 . (3.25)
Thus, thanks to this new discretization, we well obtain the Kröner link (3.24). In fact, by multiplying equation (3.21) with the Burgers vector and using this discretization (3.25), we obtain

b Ÿ i (t + ∆t) ≠ Ÿ i (t) ∆t = ≠b f fl i+ 1 2 ≠ f fl i≠ 1 2 ∆x = ≠ 2 ∆t AA ' 1+ 1 2 ≠ ' 1≠ 1 2 ∆x B t+∆t ≠ A ' 1+ 1 2 ≠ ' 1≠ 1 2 ∆x B t B
(3.26) where (•) t means the field at the time t. We well recognize an explicit time derivative of the discrete Kröner equation:

bŸ i (t)=≠2 A ' 1+ 1 2 (t) ≠ ' 1≠ 1 2 (t) ∆x B . (3.27)
It is important to note that the initial fields of each simulation have to validate this discrete equation.

A convenient way to proceed is to define an initial strain field and use equation (3.27) to compute the initial Ÿ field. numerical test is an evolution of a positive dislocation attracted by a negative dislocation. Each dislocation are represented respectively by an initial bump of density of fl + and fl ≠ .T h e s eb u m p s are localized in a 2D simulation at the index (466, 31) for the positive dislocation and (558, 31) for the negative dislocation in a box of 1024◊64 nodes. Due to this close positioning, we consider these dislocations as two infinite isolated walls of dislocations along ę2 even though there is boundary conditions. The behavior of the red discrete dislocation, predicted by equation (2.3), should be attracted by the blue discrete dislocation. As well as this red dislocation is defined by a smoothed bump, this field will be spread during its shifting. In fact, the closest side of the bump to the blue dislocation is more attracted by the other side. Figure 3.7 illustrates only the behavior of the bump of fl + but there is a similar behavior of the bump of fl ≠ . The initial bump of fl + is plotted in black and the shifted field in green is obtained from its transport equation. The two other fields, the blue dots and the red dots, are indirectly obtained from the plastic strain field with the Orowan equation followed by the Kröner equation. The red field uses the simple discretization (3.23) of the Orowan equation and the blue field uses the consistent discretization (3.25). We observe that the naive discretization of Orowan equation gives a field distinctively different than the green expected field. This artefact deteriorates strongly simulations during the time and the simulations tend to diverge.

In opposition, the blue field, obtained by the consistent discretization, gives exactly the same field than the green one.

Finally, all artefacts linked to the first system of equations (3.8) are fixed. However, others artefacts will come from the coupling with the second equation, the mechanical equilibrium (3.4). In fact, the FFT algorithms used to compute the mechanical equilibrium treat not correctly small wavelengths in input, in particular on the plastic strain field. Thus, if the plasticity model predicts small wavelengths on the plastic strain field, it will create oscillations on the output of the mechanical equilibrium which is the local stress field • . These oscillations coupled with the dynamic of transport equations will tend to numerical divergences. Chapter 4 proposes two new solvers for the second system to avoid any oscillations.

Spontaneous emergence of heterogeneities

In this part, we aim at investigating (i) whether dislocation patterning, as observed in stage I and II of hardening of monocrystal [START_REF] Kocks | Physics and phenomenology of strain hardening: the FCC case[END_REF], could emerge spontaneously from the simple physical picture given by our 2D model (see part 3.2); and (ii) if so, whether the correlations are decisive or not. For that purpose, and in the spirit of what has already been done for previous models [START_REF] Groma | Link between the individual and continuum approaches of the description of the collective behavior of dislocations[END_REF][START_REF] Groma | Investigation of dislocation pattern formation in a two-dimensional selfconsistent field approximation[END_REF][START_REF] Groma | Dislocation Patterning: From Micro-to Mesoscale Description[END_REF], we perform in this section the stability analysis of an infinite system with an homogeneous dislocation density submitted to an applied load.

Stability analysis of an homogeneous state

The stability analysis relies on the following steps: first, small perturbations are applied to the initial homogeneous fields. Then, the PDEs are linearized with respect to this small perturbations around the initial state. Considering the Fourier transforms of the perturbations, the coupled PDEs are transformed into a system of coupled ODEs that can be diagonalized to get the dispersion relations, i.e. the growth rates of the disturbance versus their wavelength. First, let us recall the basic equations of the model, that in the limits Ÿ π fl and

• ext ∫ Ÿ fl µb Ô fl 2fi(1≠‹)
,c a nb ew r i t t e na s( s e e section 3.2):

Y _ _ ] _ _ [ ≠ ˆfl(r, t) ˆt = Mb ˆx [{• ext + • sc ≠ • f ≠ • b } Ÿ(r, t) ≠ •b fl(r, t)] ≠ ˆŸ(r, t) ˆt = Mb ˆx [{• ext + • sc ≠ • f ≠ • b } fl(r, t) ≠ •b Ÿ(r, t)] (3.28)
where

•b = µb 2fi(1 ≠ ‹) 1 fl (C ++ + C +≠ ) ˆfl ˆx , (3.29) • b = µb 2fi(1 ≠ ‹) 1 fl (C ++ ≠ C +≠ ) ˆŸ ˆx , (3.30) • ext ≠ • sc ≠ • f = B A 2fi(1 ≠ ‹) µb Ô fl B -≠1 È• ext + • sc Í -, (3.31)
where Èf Í -= sign(f ) |f | -and with

• sc = ⁄ • ind (r ≠ r Õ )Ÿ(r Õ )dr Õ . (3.32)
From now on and without loss of generality, we consider a strictly positive external stress • ext .A sj u s t recalled, the analysis is performed on a model restricted to Ÿ π fl. Moreover, the self-consistent stress depends only on Ÿ. Consequently, at the initial time t 0 ,w eh a v e• sc π • ext .T h u s ,• ext + • sc (r) Ø 0 and we can safely write È• Í = • . The initial state is assumed to be homogeneous with no GND, i.e. fl(r, t 0 )=fl 0 and Ÿ(r, t 0 )=0, and a small noise is added to both fields. Next, we make appear the following quantities: a characteristic frequency

Γ µ = Mb 2 µfl 0 2fi(1 ≠ ‹) , ( 3.33)
which gives an adimensional time t = t Γ µ , (3.34) the adimensional wave vector qx = q x Ô fl 0 , (3.35) and the adimensional length x = x Ô fl 0 , (3.36) the adimensional stress

•ext = 2fi(1 ≠ ‹)• ext µb Ô fl 0 , ( 3.37)
and By applying small perturbations "fl(r)=fl(r)≠fl 0 and "Ÿ(r)=Ÿ(r), equations (3.28) can be expanded up to first order in these perturbations to give:

•ind (x, y)= x(x 2 ≠ ỹ2 ) (x 2 +ỹ 2 ) 2 , ( 3
Y _ _ ] _ _ [ ≠ ˆ"f(r, t) ˆ"t = B • - ext ˆ"Ÿ ˆx ≠ (C ++ + C +≠ ) ˆ2" f ˆx 2 ≠ ˆ"Ÿ(r, t) ˆ"t = B-• -≠1 ext ˆ• sc ˆx + B (3 ≠ -) 2 (• ext ) -ˆ"f ˆx ≠ (C ++ ≠ C +≠ ) ˆ2"Ÿ ˆx 2 (3.40)
It is worth noting that the self-consistent stress (3.32) depends on "Ÿ. This system of linear PDEs is Fourier transformed (appendix A) to obtain a system of ODEs:

Y _ _ _ ] _ _ _ [ ≠ ˆf ˆt = B • - ext iq x Ÿ ≠ (C ++ + C +≠ )(≠q 2 x ) f ≠ ˆŸ ˆt = B-• -≠1 ext iq x A ≠4ifi qx q2 y q4 B Ÿ ≠ (C ++ ≠ C +≠ )(≠q 2 x ) Ÿ + B (3 ≠ -) 2 • - ext (iq x ) f (3.41)
where we have used the adimensional Fourier transform of the individual stress

•ind (q)=≠4ifi qx q2 y q4 (3. 42
)
where i is the imaginary number, q =(q 1 , q2 , q3 ) is the adimensional wave vector and f stands for the Fourier transform of f . In matrix form, equation (3.41) reads:

ˆt A " f(q, t) " Ÿ(q, t) B = D(q) • A " f(q, t) " Ÿ(q, t) B (3.43)
where D is a 2-dimensional matrix defined by

D(q)=Mb Q a ≠(C ++ + C +≠ )q 2 x ≠iB • - ext qx ≠iB (3≠-) 2 • - ext qx ≠ 3 -B • -≠1 ext 4fi q2 x q2 y q4 +(C ++ ≠ C +≠ )q 2 x 4 R b . (3.44)
As usual, solving (3.43) amounts to find the eigenvalues of D. Its characteristic polynomial reads:

P = ⁄ 2 ≠ ⁄ • A + E (3.45)
where

A = ≠ A (C ++ + C +≠ )q 2 x + A -B • -≠1 ext 4fi q2 x q2 y q4 +(C ++ ≠ C +≠ )q 2 x BB (3.46)
and

E =(C ++ + C +≠ )q 2 x A -B • -≠1 ext 4fi q2 x q2 y q4 +(C ++ ≠ C +≠ )q 2 x B + B • - ext qx B (3 ≠ -) 2 • - ext qx . (3.47)
Finally, we obtain the eigenvalues of D

⁄ ± = A ± Ô A 2 ≠ 4E 2 . (3.48)
The real part of ⁄ ± is always negative, stating that the homogeneous state is always stable. For qx =0, ⁄ ± is equal to zero implying that the perturbations with wave vector parallel to the y axis will neither grow nor decrease.

We first consider the mean field solution i.e. correlation terms are neglected. To recover this case, we set the fitting parameters to -=1 and B =1 and suppress all the terms that depend on C ++ or C +≠ , which amounts to set C ++ = C +≠ =0. Equation (3.48) then becomes:

2⁄ ± Γ µ = ≠ sin 2 (2") 4 ± î ı Ù A sin 2 (2") 4 B 2 ≠ 1 (2fi) 2 • 2 ext q 2 x (3.49)
as predicted by [START_REF] Groma | Link between the individual and continuum approaches of the description of the collective behavior of dislocations[END_REF][START_REF] Groma | Investigation of dislocation pattern formation in a two-dimensional selfconsistent field approximation[END_REF] and where " is the angle between the x-axis and the wave vector q:

1 4 sin 2 (2")= q2 x q2 y q4 . (3.50)
We note that the real parts of ⁄ ± are always strictly negative, except when q is parallel to either the x or the y axis, in which case they are equal to zero, meaning that the corresponding perturbations are marginally stable [START_REF] Groma | Link between the individual and continuum approaches of the description of the collective behavior of dislocations[END_REF][START_REF] Groma | Investigation of dislocation pattern formation in a two-dimensional selfconsistent field approximation[END_REF]. Now, we consider -=1. negative, except when q x =0, in which case they are equal to zero. This means that perturbations with a wave vector parallel to the y axis are marginally stable.

Simulation of dislocation pattern

In this section, we present numerical simulations of the model discussed in the previous stability analysis. We also compare the results to the ones obtained with a mean field model to highlight the importance of the correlation-induced stresses. These numerical results should be considered as a preliminary study, as only a few situations have been investigated.

We have chosen a homogeneous initial state with densities fl 0 =2e 14 m ≠2 and Ÿ =0 and we have added to these fields small perturbations of the order of 10 12 m ≠2 .Ap e r i o d i c10 ≠4 ◊10 ≠4 m 2 box is discretized into 256 2 nodes such that Ô fl 0 ∆x =3.9 to be consistent with the linear stability analysis in figure 3.8. We apply a stress •ext =5.34, and we choose ∆t = ∆t crit /2 to ensure the stability of the numerical scheme.

Snapshots of fl and Ÿ are displayed in top row and bottom row of figure 3.9, respectively. At t = 347 we observe the emergence of a pattern on the GND density field. This pattern consists in small walls perpendicular to the glide direction, alternatively consisting of positive and negative dislocations. The anisotropy of the pattern is clearly evidenced in the Fourier map shown in figure 3.10. We indeed observe that the amplitude of the Fourier modes of the GND density field are confined to a small domain centred on a q-vector with a finite qx and zero qy components. We note that the emergence of this pattern contradicts the stability analysis presented above, which predicts that fluctuations with non-zero qx should collapse. Hence, non-linearities are probably already at play.

Conclusion

In a first part, we have proposed to apply an accurate numerical scheme developed to handle convection equations in fluid dynamics to the transport equations of dislocation densities. This scheme mixes ingredients from the Lax-Wendroff High Resolution scheme and from a conservative scheme. Moreover, we have shown that a particular attention must be paid to remain consistent with the transport equations when discretizing the Orowan equation to get the plastic strain field from the dislocation densities. It must be stressed that these schemes are worth to be applied to any plasticity model based on dislocation densities transport equations.

In a second part, we performed the linear stability analysis of our model, both in a mean-field approximation and in a complete model that incorporates the correlation-induced stresses (in the limit Ÿ π fl). In the mean field situation, we found that fluctuations with wave vectors parallel or perpendicular to the glide direction are marginally stable. When correlation-induced stresses are taken into account, only fluctuations with q-vectors perpendicular to the glide direction are marginally stable.

Finally, we have presented a preliminary numerical simulation of the model. We observed that the simulated microstructures do not follow the predictions of the stability analysis. Interestingly, when correlation-induced stresses are taken into account, we found that a pattern consisting of alternating positive and negative dislocation walls perpendicular to the glide direction emerges, whereas no heterogeneities show up in the mean field version of the model.

Chapter 4 NUMERICAL ESTIMATION OF MECHANICAL FIELDS FOR INHOMOGENEOUS ANISOTROPIC MEDIA

Résumé du chapitre

Ce chapitre traite de la résolution numérique de l'équilibre mécanique dans un matériau hétérogène, afin d'obtenir en tout point les champs de déformation et de contrainte. La connaissance de ces champs est utile pour de nombreuses applications et en particulier à la mise en oeuvre numérique des modèles de plasticité en densité de dislocations présentés dans les chapitres précédents. Les algorithmes que nous avons considérés dans ce travail sont basés sur une méthode de point fixe utilisant des Transformées de Fourier Rapides (FFT). Ces méthodes sont des alternatives aux Méthodes Elements Finis (FEM). Les méthodes FFT sont construites sur des grilles régulières rendant aisé un grand nombre de couplages. En particulier, ces grilles sont aussi adaptées à la discrétisation des equations de transport de densité de dislocations vue au chapitre 3. Quatre modèles sont présentés ici. Le premier, notés (1), correspond à des travaux pionniers basés sur la publication initiale de Moulinec et Suquet [START_REF] Moulinec | A fft-based numerical method for computing the mechanical properties of composites from images of their microstructures[END_REF] et nous étudions également une variante de ce modèle utilisant une méthode de différences finies pour l'évaluation des opérateurs différentiels (modèle (2)). Ensuite, deux nouveaux modèles, notés (3) et (4), sont proposés sur des grilles décalées. Ce sont des méthodes de différences finies permettant de résoudre les artefacts numériques des modèles (1) et (2) conduisant à des oscillations sur les champs mécaniques. Quand les constantes élastiques sont isotropes (ou anisotrope alignée sur la grille), les modèles (3) et (4) sont identiques et permettent l'obtention de champs mécaniques non oscillants même pour de très fortes inhomogénéités (milieux poreux). Les modèles (3) et (4) différent par leur traitement des anisotropies élastiques: le modèle (3) utilise des interprétations entre sites premiers voisins alors que le modèle (4) se base sur une construction variationelle pour laquelle l'énergie élastique est décomposée sur les subvoxels définis par les grilles décalées. En utilisant des configurations modèles, nous montrons que seul le modèle (4) permet d'obtenir des champs mécaniques non oscillants quelque soit l'anisotropie et l'inhomogénéité.

Introduction

As mentioned before in chapter 3, the numerical solution of the proposed dislocation density model requires the knowledge, at each time step, of the stress field in the material. This field can be obtained by numerically solving the mechanical equilibrium. An efficient numerical solver is then needed. This chapter is devoted to the proposition of a new solver. Note also that an efficient solver could also be of interest in other situations, such as for the homogenization of the mechanical properties of an heterogeneous material.

A classical way of solving mechanical equilibrium is to use a Finite Element Method (FEM). However, for the case of an evolving microstructure, the numerical mesh has to be permanently re-evaluated in order to match to the evolutive microstructure. This remeshing is usually very time consuming. Another way has been proposed by Moulinec and Suquet [START_REF] Moulinec | A fft-based numerical method for computing the mechanical properties of composites from images of their microstructures[END_REF] using a scheme based on Fast Fourier Transforms (FFT). This method discretizes the media on a regular grid rather than a mesh and uses the efficiency of FFT. Moreover, this regular grid can be easily coupled with a wide range of information such as experimental images. The FFT-based solver is also commonly used in phase field models describing elastically inhomogeneous microstructures [START_REF] Boussinot | Phase-field simulations with inhomogeneous elasticity: Comparison with an atomic-scale method and application to superalloys[END_REF][START_REF] Cottura | Role of elastic inhomogeneity in the development of cuboidal microstructures in ni-based superalloys[END_REF], heterogeneous plasticity [START_REF] Geers | Coupled glide-climb diffusion-enhanced crystal plasticity[END_REF] or viscoplasticity [START_REF] Idiart | Infinite-contrast periodic composites with strongly nonlinear behavior : Effective-medium theory versus full-field simulations[END_REF] or crack description [START_REF] Herrmann | Linear and elastic-plastic fracture mechanics revisited by use of fourier transforms -theory and application[END_REF]. Finally, these FDM algorithms are built with the fixed point method which gives very simple algorithm. The drawback of this method is on the quality of the description of interfaces only localized on regular nodes. After the seminal paper of Moulinec and Suquet, several improvement have been proposed for the FFT-based solvers. First, an enhanced algorithm, proposed by [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF], accelerates the convergence of the iterative fixed point. Then, [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF] proposes an other building with an augmented Lagrangian which allows infinite inhomogeneities of elastic stiffness such as for porous materials [START_REF] Monchiet | Numerical homogenization of nonlinear composites with a polarizationbased ffti t e r a t i v es c h e m e[END_REF]. Then, a variational framework has been proposed in [START_REF] Brisard | Fft-based methods for mechanics of composites: A general variational framework[END_REF] to handle infinite elastic contrast. However, all these methods are often only discussed in the context of the homogenization of mechanical properties and almost all applications are performed using isotropic elasticity. Because we are interested here in alloy where the microstructure may be simultaneously elastically inhomogeneous and anisotropic, we have to discuss the accuracy of these FFT-based methods for the determination of the full mechanical fields in such a complex alloys. A first step towards this goal has been very recently made by F. Willot who showed that classical FFT schemes may lead to oscillating fields [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF], and a new spatial discretization grid was proposed to reduce the numerical oscillations. The aim of the present chapter is to propose an alternative and very efficient FFT scheme able to obtain accurate and non oscillating mechanical fields. This chapter investigates several FFT-based schemes. We distinguish two groups of discretization grids. Section 4.2.1 presents two classical methods, (1) and (2), using a single regular grid based on works of [START_REF] Moulinec | A fft-based numerical method for computing the mechanical properties of composites from images of their microstructures[END_REF]. Section 4.2.2 proposes two alternative methods, (3) and (4), using staggered grids inspired, among others, by geophysical works [START_REF] Bale | Staggered grids for 3d pseudospectral modelling in anisotropic elastic media[END_REF][START_REF] Saenger | Modeling the propagation of elastic waves using a modified finite-difference grid[END_REF][START_REF] Zeng | A staggered-grid finite-difference method with perfectly matched layers for poroelastic wave equations[END_REF] or coherency loss [START_REF] Geslin | Investigation of coherency loss by prismatic punching with a nonlinear elastic model[END_REF]. Methods (1), ( 2) and ( 3) are implemented with a fixed point. Method ( 4) is implemented with a direct solver but the equivalent fixed point algorithm in Fourier space is also expressed. In order to emphasize the importance of a variational scheme (method (4)), part 4.3 and part 4.4 present tests for the comparison of all methods and expose the benefits of staggered grids and especially the necessity of a physically consistent algorithm.

Numerical methods: classical and alternative approaches

In this part, we present several Fourier-based numerical schemes for solving inhomogeneous and anisotropic elasticity problems. As we consider linear elasticity, the problem to be solved can be expressed as

‡ ij (r)=C ijkl (r)(' kl (r) ≠ ' ú kl (r)) ; rj ‡ ij (r)=0 (4.1)
where r =(r 1 ,r 2 ,r 3 ) describes the space in three dimensions, C ijkl (r) is the elastic stiffness tensor and ‡ ij (r) the stress tensor. In equation (4.1), the summation over repeated indices is assumed, ' ú kl (r) is a local eigenstrain tensor which may originate from a phase transformation or from plasticity, and ' ij (r) is the total strain tensor which derives from the total displacement field ų =(u 1 ,u 2 ,u 3 ):

' ij (r)= 1 2 A ˆui (r) ˆrj + ˆuj (r) ˆri B . ( 4.2)
As we consider a system with periodic boundary conditions, the usual boundary condition is replaced by a condition on the average value of the strain or stress fields. In this paper, we hold the system at a constant average strain È' ij Í but the extension to the constant average stress condition is straightforward.

All the following schemes use an homogeneous reference material with a stiffness tensor C 0 ijkl , and the local stiffness heterogeneities are

"C ijkl (r)=C ijkl (r) ≠ C 0 ijkl (4.3)
where C 0 ijkl influences the convergence of the algorithm and a good choice appears to be a linear average of matrix stiffness C m ijkl and precipitate stiffness C p ijkl [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. In order to benefit from the efficiency of the FFT algorithms, we use regular grids and the schemes are written in Fourier space. For simplicity, we consider a grid spacing d equal in each direction. The simulation box size in each direction is then L i = N i d where N i is an integer, assumed even. The definition of the Fourier transform ' f (q) of a periodic function f (r) is recalled in appendix A.

Single grid

In the first two methods, we use a single regular grid for all mechanical fields.

Method (1): Continuum equation of equilibrium in Fourier space

As proposed by [START_REF] Moulinec | A fft-based numerical method for computing the mechanical properties of composites from images of their microstructures[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF], the differential problem is transformed into an integral problem using the Lippmann-Schwinger equation. In Fourier space, we obtain

' ' ij (q)=≠Γ ijkl ' • kl (q) (4.4)
where the hat sign denotes the Fourier transform and the polarisation tensor is

• ij (r)="C ijkl (r) ' kl (r) ≠ C ijkl (r) ' ú kl (r). (4.5)
The strain Green tensor is

Γ (1) ijkl (q)= 1 4 Ë1 q j G (1) ik + q i G (1) jk 2 q l + 1 q j G (1) il + q i G (1) jl 2 q k È (4.6)
where G ik stands for the displacement Green tensor defined by

G (1) il ≠1
(q)=C 0 ijkl q j q k . (

Because periodic boundaries are assumed, the Fourier space is discrete. To use FFT algorithms, the Fourier space is limited to the first Brillouin Zone ' B 0 defined as the set of vectors 2fi(-

1 /L 1 , -2 /L 2 , -3 /L 3 ) where -i are integers such as ≠N i /2 < -i AE N i /2.
The problem is then solved using the following fixed-point algorithm:

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ -Initialization: ' n=0 ij (r)=È' ij Í
-Iteration: 1: Calculation of the polarization tensor:

• ij (r)="C ijkl (r)' n kl (r) ≠ C ijkl (r)' ú kl (r) 2: ' • (q) © Fourier transform of • (r) 3: Actualization: ‰ ' n+1 ij (q)=≠Γ ijkl (q) ' • kl (q)( q " = 0) 4: ' n+1 (r) © Inverse Fourier transform of ‰ ' n+1 (q)
Steps 1 to 4 are repeated until convergence and we choose a practical stopping criterion on the variation of the deformation fields between two iterations:

e ' = ÈÎ' n+1 (r) ≠ ' n (r)ÎÍ Î' ref Î (4.8)
where ' ref is a reference deformation tensor, È•Í denotes the spatial average, and where to Î'Î = Ò q ij ' 2 ij is the norm of the tensor. In the following simulations performed at fixed average strain, we use ' ref = È'Í when this quantity is non zero; and ' ref = ' ú otherwise. When a typical value ‡ ref for the stress fields in the material is known, the more physically based criterion

e ‡ =  ÈÎdiv ‡Î 2 Í ÎC m • ' ref Î (4.9)
is also computed to compare the different schemes and discretizations.

Due to the non-symmetric choice of the Fourier domain ' B 0 ,thesymmetryΓ ijkl (≠q)= Γijkl (q) is not ensured at the boundary (ā stands for the complex conjugate of a). This property can be enforced by replacing Γ ijkl (q) by (Γ ijkl (q)+ Γijkl (≠q))/2 or by C 0 ijkl ≠1 at the boundary where Γijkl is the complex conjugate of Γ ijkl . Conforming to [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF], utilization of C 0 ijkl ≠1 decreases significantly oscillations in case of small eslastic inhomogeneities. In opposition, for strong inhomogeneities, this choice compromises the stability of the algorithm. Because this paper present results on very inhomogeneous microstructures, such as materials containing pores, we have used the first choice.

Method (2): Fixed point algorithm on discretised equation (4.10)

The second FFT based method is built on the discretized mechanical equilibrium in real space on a regular grid B 0 :

D (ii) j [C 0 ijkl 1 D (ii) k [u l ] ≠ ' ú kl 2 + "C ijkl (' kl ≠ ' ú kl ) ] =0 (4.10)
where

D (ii) i [ • ]
is the central difference along the direction i between two nodes distanced of 2d. Going to Fourier space, an expression similar to equation (4.4) is obtained for each vector q of ' B 0 but the symmetrised strain Green tensor is now given by

Γ (2) ijkl (q)= 1 4 C 3 G (2) ik q (2) j + G (2) jk q (2) i 4 • q (2) l + 3 G (2) il q (2) j + G (2) jl q (2) i 4 • q (2) k D (4.11)
where q

(2)

i = sin q i d d . The displacement Green tensor Gij is defined inside the first Brillouin zone ' B 0 by G (2) il ≠1 (q)=C 0 ijkl q (2) j q (2) k . (4.12)
At the boundary of the first Brillouin zone, G

ij is set to zero. Then, a fixed point scheme similar to the one presented in subsection 4.2.1.1 is used to obtain elastic equilibrium.

When considering diffuse interface approaches, such as phase field models [START_REF] Boussinot | Phase-field modeling of bimodal microstructures in nickelbased superalloys[END_REF][START_REF] Cottura | A phase field model incorporating strain gradient viscoplasticity: Application to rafting in Ni-base superalloys[END_REF][START_REF] Cottura | Phase field study of acicular growth: Role of elasticity in widmanstatten structure[END_REF][START_REF] Cottura | Role of elastic inhomogeneity in the development of cuboidal microstructures in ni-based superalloys[END_REF], only small q vectors are used (typically |q| < fi/(5d)), therefore q (2) i ƒ q i and the definitions of the strain Green tensor (4.6) and (4.11) lead to very similar results.

Staggered grids schemes

In this section, we propose alternative schemes based on staggered grids [START_REF] Geslin | Investigation of coherency loss by prismatic punching with a nonlinear elastic model[END_REF][START_REF] Bale | Staggered grids for 3d pseudospectral modelling in anisotropic elastic media[END_REF][START_REF] Saenger | Modeling the propagation of elastic waves using a modified finite-difference grid[END_REF][START_REF] Zeng | A staggered-grid finite-difference method with perfectly matched layers for poroelastic wave equations[END_REF] with a particular emphasis on the anisotropic elasticity cases. Method (3), easier to implement, is adapted for restrained anisotropic cases and method (4) is powerful for any anisotropy. The characteristics of the staggered grid is to define each component of the displacement vector (u 1 ,u 2 ,u 3 ) on different grids, as shown in figure 4.1. The aim of this choice is to define the strain components using a central difference scheme based on nearest neighbors nodes, in order to avoid the short wavelength oscillations observed with classical central difference schemes. Figure 4.1 also indicates the grids on which the strain components are defined. Due to the staggered grids of the displacement, the central differences defining the strain tensor are also on staggered grids. The diagonal terms ' ii are all on the circle-grid (in blue). The non diagonal components ' ij are defined on a grid which is shifted by a vector d 2 (ę i + ęj ) with respect to circle grid. For example, the ' 12 components are known on the cross-grid (in red) in figure 4.1. More generally, the components of all the rank-2 symmetric tensors A ij (e.g. ', ' ú , • ,...) a r ed e fi n e d on the same grids than the ones used for '. The mechanical equilibrium (4.1) and the definition of the polarization tensor contain terms of the type the sum of terms A kl which may not be on the same grid. In order to solve this problem, the A kl components need to be interpolated on the grid on which B ij is defined. In this work wa have two ways of performing these interpolations. The first one is based on a nearest neighbor scheme (Method (3)) and the second one is obtained from a variational principle. These two approaches are successively presented in the following subsections.

B ij = C ijkl A ij . Therefore the component B ij is v - ' v - ' v - ' v - ' v T c v T c v T c v T c j j j j d d
(α, β) (α +1, β) (α, β +1) (α +1, β +1) 12 12 12 12 (α -1, β) (α, β) (α -1, β +1) (α, β +1) u 1 (α -1, β) u 1 (α, β) u 1 (α -1, β +1) u 1 (α, β +1) u 2 (α, β) u 2 (α +1, β) u 2 (α, β +1) u 2 (α +1, β +1) - ' d A ii A 12 A 13 A 23

Method (3): Linear interpolation of stress components

Method (3) is a FFT based method built from the discrete mechanical equilibrium where gradient operators are computed using central differences between two nearest neighbor nodes. This central difference along i direction is written

D (i) i [ • ]: D (i) j [C 0 ijkl 1 D (i) k [u l ] ≠ ' ú kl 2 + "C ijkl (' kl ≠ ' ú kl ) ] =0. (4.13)
In Fourier space this operator reduces to multiplication by iq

(3) j ,w h e r eq (3)
j =sin(q j d/2)/(d/2). As expected, q(3) approaches the wave vector q as |q| tends to 0. As explained above, we need an interpolation scheme to estimate the values kl component of a rank-2 symmetric tensor on the grid where ij components are defined. As detailled in appendix C, this can be done by averaging the values of the kl components on the four sites that in a first nearest neighbor position from the considered site of the ij components. We therefore introduce the interpolation operator I (ij,kl) (A kl ) defined as Thus, we express the energy of our system by

I (ij,kl) (f )= ij\kl (i, i)( 2 , 3) (1, 3) (1, 2) (i, i) (2, 3) (1, 3) (1, 2) Q c c c c a fL 0 11 L 10 1 L 11 0 L 011 fL 110 L 101 L 101 L 1 10 fL 0 11 L 110 L 10 1 L 01 1 f R d d d d b (4.14) n n n n @ @ @ ◊ ' ú ij ◊ C ijij ◊ ' ú ij ◊ C ijij ◊ ' ú ij ◊ C ijij ◊ ' ú ij ◊ C ijij ¶ ' ú ii ¶ C iijj ¶ C iiij ¶ ' ú ii ¶ C iijj ¶ C iiij ¶ ' ú ii ¶ C iijj ¶ C iiij ¶ ' ú ii ¶ C iijj ¶ C iiij
F el = 3 d 2 4 D ÿ R ÿ s ÿ ijkl 1 2 C ijkl ( R, s)' el ij ( R, s)' el kl ( R, s) (4.21)
detailed in appendix D. The energy density is considered homogeneous in each subvoxels. The subscript D is the dimension of our space and 'el ij ( R, s) is linked to the displacement detailed in appendix B. The dependences R and s scan each subvoxels. The strain in a subvoxel 'el mn ( RÕ , s) is related to the strain ' el mn ( RÕ ) defined on the staggered grid (figure 4.1) by

'kl ( RÕ , s)=' kl ( RÕ ≠ ◊kl (s)) (4.22)
where ◊kl (s) is the shift vector defined by

Y _ ] _ [ ◊kl (s)= d 2 (ę k + ęl ) ≠ 2[(s • ęk ) • ęk +(s • ęl ) • ęl ] for k " = l ◊kk (s)=0 (4.23)
The last precision of equation (4.21) is the choice of the elastic coefficient C ijkl used for a subvoxel. In order to illustrate, we express the energy density of one subvoxel by splitting in three different contributions according to stiffness tensor symmetries:

f el = f iijj el +2f i" =j ijij el +4f i" =j iiij el with f ijkl el = 1 2 C ijkl ' el ij ' el kl . (4.24)
The choice of the stiffness coefficient for each contribution, summerized in figure 4.2, is given by:

-The first contribution handles only diagonal strain terms defined on the circle-grid. Our choice is to use the elastic stiffness of the same grid ( ¶ C ijij ).

-The second contribution handles only non-diagonal term defined on the cross-grid. Our choice is to use the elastic stiffness of the same grid ( ◊ C iijj ).

-The third contribution handles diagonal but also non-diagonal terms. Our choice is to use the elastic stiffness defined on the circle-grid ( ¶ C iiij ).

Signs ¶, ◊ designate respectively fields on the circle-grid and on the cross-grid.

In order to be consistent with the periodic conditions, we split the displacement field: ų = ų0 +"ų where ų0 = È'Í • r is the linear part, which is non-periodic and "ų is the fluctuation around this linear part. Thus, "ų can be consistent with the periodicity. After substracting ų0 and applying the second law of thermodynamics, we obtain the following linear system on the fluctuation "ų:

A ij • "U j = B i with "U j©{k,-,-} = "u k (-, -), B i©{l,-Õ ,-Õ } = ˆFel ˆul (-Õ , -Õ ) , A ij©{l,-Õ ,-Õ }{k,-,-} = ˆ2F el ˆul (-Õ , -Õ )ˆu k (-, -) , ( 4.25)
where U is a 1-dimensional matrix containing fluctuation with an adapted storage indexed by j (j scan all component k and all discrete spaceand -). B is also a 1-dimensional matrix containing the first derivative of the energy by the displacement and indexed by i (i is a similar index than j). Then, A is a 2-dimensional matrix, the Hessian matrix, containing the second derivative of the energy by the displacement indexed by i and j. The first and second derivation of the energy are given respectively by appendix D and appendix E. Note that, because of a translational invariance by ų, the system have an infinite solutions. In this way, to restrain to one solution, we constrain ų(-=0, -= 0) to zero. We solve the system with the MUMPS library which is a fast numerical solver for sparse system using Approximate Minimum Degree automatic quasi-dense row detection method (QAMD). The direct solving appears to be faster than the fixed point method when we need more than the order of hundred iterations. However, direct solvers can be easily massively parrallelized.

An alternative solving is to implement a fixed point algorithm similarly to the third first methods. Based on the same staggered grids, the algorithm is quite similar to the method (3) with two small modifications. Firstly, the Green tensor of the displacement Γ Secondly, the polarization tensor of a subvoxel is now

•ij ( R, s)=" Cijkl ( R, s)' kl ( R, s) ≠ Cijkl ( R, s)' ú kl ( R, s). (4.28)
For aligned case, the algorithm of method ( 3) and ( 4) are equivalent. Note the method (4) has been very recently (August 2015) and independently proposed in [START_REF] Schneider | Computational homogenization of elasticity on a staggered grid[END_REF]. However, no anisotropic simulations are presented in this reference.

Comparison of methods: Case of cubic symmetry aligned with the grid

This part treats the restricted 'aligned case'. As mentioned in subsection 4.2.2.1, in this case, there is no stress interpolation in method [START_REF] El-Azab | Statistical mechanics of dislocation systems[END_REF]. In this way, method (3) is similar to method (4). Thus, only methods (1), ( 2) and (3) will be compared. Firstly, we illustrate intrinsic artefacts with the simple case of an homogeneous material, section 4.3.1. In a second time, we observe a strong degradation in case of inhomogeneous material, section 4.3.2. These analysis are performed in two dimensions and should be equivalent in 3D. Simulations are achieved with a simple microstructure containing square precipitates in a periodic matrix. Square is favored because it stays perfectly described in a discrete space and moreover it owns a singularity at each corner. In a continuum theory, mechanical fields diverge at these singularities, which is an interesting problem for schemes behavior analysis. The width of the precipitate is equal to a quarter of the matrix period. The single grid attached to method (1) and ( 2) has the same position of the round-grid attached to method (3). Moreover, the interface of the precipitate is sharp and positioned on the cross-grid. That is why, a linear interpolation of matrix and precipitate stiffness tensor is taken for ◊ C ijkl at the interface. The appendix F gathers simulations of the situation with a wide range of parameters in order to have exhaustive comparisons. The following analysis are based on these comparisons.

Illustration of classical scheme artefact

This section considers homogeneous stiffness fields with Lame coefficient ⁄ = 100 GPa and µ = 50 GPa. The precipitate is affected by an eigenstrain ' ú 11 = ' ú 22 =0.01. Fixed point is always carried out up to the convergence (e ' 6 10 ≠10 ). For each method we proceed to different discretizations, called N-run where N characterizes the 2D-discretization: N 2 = {64 2 ; 1024 2 }. Some others discretizations are proposed by appendix F. We observe that method (1) and (3) converge toward the same solution and method (2) oscillates around this solution. This oscillation, coming from classical central difference, stays for any discretization. Note that, this oscillation looks like steps and is weaker in specific case of square precipitate described with even nodes.

The second artefact, localized around corners, is attached to the curvature of the interface. In order to quantify this artefact, figure 4.4 shows the error of the solution of 64-run compared to a perfect reference. The maps are centered around one corner and show 20◊20 nodes, as illustrated on the map of figure 4.3. We take the higher resolution (1024-run with method (1)) as reference and, in order to compare, we average it on a 64◊64 grid. The overall pattern of error looks like the same for all methods but with additional oscillation for method (1) and ( 2). This error analysis is also done for 256-run (256◊256 grid) and we observe the same pattern, rigorously the same for method (1) and (3) and quite similar for method (2) du to it non-local artefact. Consequently, this pattern depends only on the curvature of the discrete interface. In this way, we can define a characteristic length of this error. For the square case, oscillations spread up to six nodes.

The characteristic amplitude of an artefact is quantified by the standard deviation of this error on the entire space which is 3.5e ≠6 , 1.5e ≠4 and 8.4e ≠6 , respectively for method (1), ( 2) and (3). We notice that method (2) is strongly worse than method (1) and (3) and method (1) appears to be the best. We observe the same behavior for ' 12 analysis and for any discretization (16-run, 64-run, 256-run). Finally, these artefacts remain very soft for method (1) and (3) but will strongly increase for method (1) in case of inhomogeneity. Note that, recent works propose solutions to limit oscillations with new schemes tested for 'aligned case' [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF][START_REF] Gélébart | Filtering material properties to improve fft-based methods for numerical homogenization[END_REF][START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF].

To conclude the case of an homogeneous media, we propose the classical test of two dislocations inside the same slip system, gliding along ę1 and with opposite Burgers vectors. The Kröner equation (3.24) gives the equivalent strain dilatation ' ú xy defined between these two dislocations. This strain field is spread on 16 nodes along ę1 and one node along ę2 . Figure 4.5 illustrates the resulting fields of ‡ yy obtained, in the order, by method (1), ( 2) and ( 3). This test is interesting because we observe Method ( 1) Method ( 2) 2) and (3). The two dislocations are positioned at [START_REF] Daveau | Interaction dislocations -joints de grains en déformation plastique monotone : étude expérimentale et modélisations numériques[END_REF][START_REF] Fressengeas | An elasto-plastic theory of dislocation and disclination fields[END_REF] and [START_REF] Geers | Coupled glide-climb diffusion-enhanced crystal plasticity[END_REF][START_REF] Fressengeas | An elasto-plastic theory of dislocation and disclination fields[END_REF]. The profil, coming from the method (3) and illustrated by the dotted line on the third map, is compared to the analytical solution. a significant oscillations on method (1) and (2) which emphasize the benefit of method (3). The field ‡ yy created by an edge dislocation can be obtained analytically. The profil on figure 4.5 is compared to the mathematical solution. This profil goes through the negative dislocation as illustrated by the dotted line on the map of method (3). We observe a perfect matching excepted near the dislocation core and at the boundary of the box. Firstly, the discrete aspect cannot reproduce the divergence predicted approaching the core. Secondly, due to the symmetry of the periodicity the strain field is imposed to zero at the boundary whereas a non zero field in case of an infinite empty space.

2.03e -1 1.58e -1 1.13e -1 6.75e -2 2.25e -2 -2.25e -2 -6.75e -2 -1.13e -1 -1.58e -1 -2.

Influence of elastic constants inhomogeneity

We keep the previous microstructure configuration with a modification of Lame coefficient in order to have a quasi-porous microstructure: ⁄ m = 100 GPa and µ m = 50 GPa for the matrix and (⁄ p ,µ p ) = 10 ≠3 • (⁄ m ,µ m ) for the precipitate. Moreover, we impose a macroscopic strain È' 11 Í = È' 22 Í =0.01 and we have no eigenstrain. The oscillation of the corner artefact appears to have the same characteristic length as homogeneous case but with a sizable amplitude. In addition, the opposite case with hard precipitates inside a soft matrix is proposed in appendix F (contrast equals to 10 3 ). We observe significant oscillations in each case on all mechanical fields (' 11 , ' 12 , ‡ 11 , ‡ 12 ). The worst degradation appears always inside the concave part (precipitate), on strain field for quasi-porous case and on stress field for soft matrix case. The alternative method (3) seems to fix these oscillations for any contrast and any discretization.

Comparison of methods: General anisotropic case

This part treats general anisotropy holding new contributions (C i" =j iiij " =0 )w h i c he m p h a s i z e s importance of interpolations. Firstly, a comparison between stress interpolation (method (3)) and energy based interpolation (method (4)) is illustrated in section 4.4.1. Then, we illustrate a realistic issue in titanium alloys which is reachable by method (4) (section 4.4.2).

Benefit of a physical interpolation (method (4)) in case of strong anisotropic inhomogeneity

We still keep the previous microstructure configuration with square precipitate but with a cubic crystalline material. Moreover, the cubic symmetry and the grid symmetry form an angle of 18 ¶ which gives the following stiffness tensor for the matrix in the natural basis of the grid (Mandel notation in GPa): First point, method (3) presents a vertical blue line localized at the interface. This permanent local artefact is linked to the discretization d. Regarding the global response, we notice a worse accuracy for method (3) compared to method (4). Method (3) gives distinctively a different fields for 64-run, 256-run or 1024-run, especially inside the precipitate. That means a requirement of a very high discretization of the space to converge to the solution. In opposition, method (4) gives a good solution with 64-run which is qualitatively the same field than 1024-run. The standard deviation of the error (defined in section 4.3.1) for 64-run is about 3.6e ≠2 for method (3) and 3.4e ≠3 for method (4) which confirm the poor accuracy of method [START_REF] El-Azab | Statistical mechanics of dislocation systems[END_REF]. Method (4) appears to be efficient for any contrast, such as hard precipitate in a soft matrix.

C m = Q c c c c c c c c

Voronoï microstructure for titanium alloy containing a subsurface crack

Dwell fatigue on titanium alloy, imposed by engine cycles, owns internal cracks growth which are strongly dependent of the characteristic direction of the microstructure. In fact, the particular grain orientations due to the fatigue allows propagation of cracks through interfaces on several phases [START_REF] Islchak | A simple model to account for the role of microtexture on fatigue and dwell fatigue lifetimes of titanium alloys[END_REF]. In order to predict the wide path of cracks we need to know the mechanical full-fields attached to an heterogeneous anisotropic microstructure on several phases. Finite element method, the classical way to analyse such a problem, describes very well the vacuum of a cracks inside a microstructure. Powerful in two dimensions, the meshing of a large number of phases in 3D involves costly simulations. Method [START_REF] Acharya | A model of crystal plasticity based on the theory of continuously distributed dislocations[END_REF] proposes a new way to obtain a high accuracy mechanical field which is very competitive. The right side of figure 4.8 shows the Von Mises stress inside the microstructure. The first half shows the Voronoï microstructure and the second half shows a monocrystal. Heterogeneity of the microstructure alters significantly the response, especially around the crack tip which will lead the crack propagation.

Conclusion

In the previous chapter, we have shown that the numerical situation of the dislocation density model requires the use of an efficient solver able to provide the stress field in an arbitrary microstructure. In this chapter, we have tested four FFT-based schemes for solving mechanical equilibrium in inhomogeneous and anisotropic microstructures. We have analyzed the full mechanical fields (strain and stress) and a particular attention has been paid to the numerical oscillation that may arise. Firstly, we consider methods based on a single grid coming from the initial work [START_REF] Moulinec | A fft-based numerical method for computing the mechanical properties of composites from images of their microstructures[END_REF]: Method ( 1) is based on a continuum equilibrium equation in Fourier space and Method (2) is based on Finite Difference Method for computing differential operators. These methods lead to significant numerical oscillations whose amplitude increases with the elastic inhomogeneity. Then, we propose two new numerical schemes built on staggered grids. In method (3) the equilibrium condition div ‡ =0 is discretized on the staggered grids. In case of isotropic elasticity (or simple symmetries aligned with the grid), we obtain two independent conditions on the two different grids. In that situation, we show that the numerical solution leads to accurate non oscillating mechanical fields. In the case of an arbitrary symmetry, interpolations of quantities from one grid to the other are required. using the most simple interpolation scheme using the nearest neighbor sites, method (3) has been implemented and tested. We obtain non oscillating fields but in highly inhomogeneous and anisotropic cases, the accuracy of the field may be moderate. Finally, we propose another method where discretization is performed on the same staggered grids. In this approach, the space is divided into subvoxels and a discrete elastic energy is written. Minimization of this discrete energy leads to a new scheme, formally very similar to the one of the method (3). Using model configurations, we show that this scheme leads to both non oscillating and accurate mechanical fields for any inhomogeneous anisotropic media such as quasi-porous material or hard precipitate in soft matrix. Method (1), ( 2) and ( 3) are implemented with a fixed point algorithm and method (4) with a direct solver. To conclude, this new scheme is very promising for new issues, for example for dwell fatigue with crack propagation in complex 3D microstructures.

OVERALL CONCLUSION

Plasticity of crystalline solids is mainly due to the movement dislocations. However, even today, conventional plasticity theories use mesoscopic variables and evolution equations that do not involve dislocations. This strongly limits the validity range of these theories. In this work, we have followed the route opened by I. Groma in order to derive a physically based plasticity model using dislocation densities. This approach, consists in explicitly performing the coarse graining procedure of the kinetics of discrete dislocations.

In a simple 2D geometry of parallel edge dislocations, this approach leads, at mesoscale, to transport equations of dislocation densities in which all ingredients can be related to the underlying physics of individual dislocations.

In the first chapter, we have clarified the mathematical procedure needed to coarse grain the dislocation dynamics from the discrete to the continuum. We have in particular emphasized that the coarse graining procedure requires a space and time convolution, supplemented by an average on a statistical ensemble. The coarse graining procedure generates correlation-induced stresses • s f and • s b that depend on the signe s of the dislocations and that have specific physical interpretations:

• s f is a friction stress because it is anti-symmetric and always opposes the local stress • ; • s b is a backstress because it is invariant with respect to a change of sign of • and thus generates a translation of the elastic domain.

We have also shown that the friction and back stresses • s f and • s b , which are sign dependent, break the sign reversal symmetry present at the discrete scale: at mesoscale, dislocations of opposite sign do not display opposite velocities.

We have also shown that the sign-independent component of • s b ,w h i c hi su s u a l l yl i m i t e dt oa term that depends on the gradient of the GND density, contains also a term that depends on the gradient of the total density.

The second chapter has been devoted to the computation of the correlation induced stress terms • s f and • s b using a dynamics of discrete parallel edge dislocations. Numerical difficulties arises from the wide range of stresses generated by a dislocation as well as from their long-range character. An efficient numerical scheme has been devised, tested and implemented. Correlation maps of the dislocation positions have been extracted from these simulations, and these maps have them been used to compute the correlation induced stress terms • s f and • s b . We have observed a dependence of these terms to the size L of the coarse graining procedure. Generally speaking, an L dependence may not be surprising because dislocation patterns are often observed on characteristic sizes much larger than the average dislocation spacing. In the present situation, where we have a single glide system with parallel glide dislocations, the L dependence has been linked to the dynamical formation of short range dipole associated to the spatial and time scales of the coarse-graining procedure. However, this L dependence remains to be more clarified and quantitatively analyzed.

Chapter 3 addresses the question of the numerical solution of the dislocation density kinetic model. It is shown that the numerical implementation requires a well controlled scheme. We have used a multigrain spatial discretization where fluxes are computed half ways from the grid points on which dislocation densities are defined. Transport equations are solved using a scheme inspired by Lax-Wendroff High Resolution scheme in order to remove the numerical diffusion observed in more simple schemes. Finally, we have also shown that the numerical scheme used to solve Orowan equation must be consistent with the one used for the transport equations of dislocation densities, in particular to ensure that Kröner equation is always fulfilled.

We have then selected a simple dislocation density model in which a friction term and a backstress have been derived from the coarse graining procedure. We first performed a linear stability analysis of our model, both in a mean-field approximation and in a complete model that incorporates the correlation-induced stresses (in the limit Ÿ π fl). In the mean field situation, we found that fluctuations with wave vectors parallel or perpendicular to the glide direction are marginally stable. When correlation-induced stresses are taken into account, only fluctuations with q-vectors perpendicular to the glide direction are marginally stable.

Finally, we have presented a preliminary numerical simulation of the model. We observed that the simulated microstructures do not follow the predictions of the stability analysis. Interestingly, when correlation-induced stresses are taken into account, we found that a pattern consisting of alternating positive and negative dislocation walls perpendicular to the glide direction emerges, whereas no heterogeneities show up in the mean field version of the model.

The fourth chapter proposes a new FFT-based scheme for numerically solving mechanical equilibrium in an heterogeneous material. This scheme used a multigrain spatial discretization (staggered grids) and leads to the minimization of a well defined discrete elastic energy. The results show that the new scheme is the only one to provide, in all test configurations, accurate and non oscillatory mechanical fields even in the case if very inhomogeneous and anisotropic elasticity. The accuracy of the new scheme is tested on model configurations and the results are compared to the one obtained with elastic solvers from the literature. This scheme uses a multigrain spatial discretization (staggered grids) and is based on the minimization of a discrete elastic energy. The accuracy of the new scheme is tested on model configurations and the results are compared to the ones obtained with other FFT-based solvers from the literature. The results show that the new scheme is the only one to provide, for all considered configurations, accurate and non oscillatory mechanical fields even in the case of very inhomogeneous and anisotropic elasticity. In addition, this scheme has the advantage ti make and well defined assumptions for the elastic behavior of interfaces.

Several perspectives can be put forward for this work. First, the new FFT scheme proposed in this work to solve the mechanical equilibrium could be compared to the recent approach proposed by [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] based on a rotated staggered grid. Preliminary investigation using isotropic elasticty have shown that the scheme proposed by W. Willot still leads to slightly oscillated field as opposed to the new scheme [START_REF] Schneider | Computational homogenization of elasticity on a staggered grid[END_REF]. Second, in the case of high elastic inhomogeneity, the convergence of the scheme could be improved for example using the accelerated scheme proposed in [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF], or going to more advanced approaches [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF]. The efficiency of the algorithm could also be improved using a Fixed Point Method on the dis-placement field variables, as for example in [START_REF] Boussinot | Phase-field simulations with inhomogeneous elasticity: Comparison with an atomic-scale method and application to superalloys[END_REF]. This will save computer memory but also computer time as less FFT will be necessary for each iteration.

The formalism of a physical crystal plasticity model developed in this thesis is a promising way but stays currently far from its accomplishment. A realistic 3D case involves much more physics, such as curved dislocation lines, climbing, cross-slip, etc. not considered for the moment. The overall building of such a 3D model should be performed with the two same stages as our simplified 2D case: firstly, a mathematical coarse-graining procedure (in the spirit of the one developed in chapter 1) and secondly, the estimation the correlation-induced local terms that inevitably will emerge (in the spirit of chapter 2). But each stage will be of course necessarily adapted to the more complex 3D situation. We attempt here to foresee the obstacles of each stage.

Here, we propose a reflexion for the first stage. In a 3D space, dislocations are curved, as illustrated in figure I.1. Consequently, the coarse-graining formalism needs to incorporate this curved aspect. Some research teams propose already different ways for the development of a curved dislocation formulation [START_REF] Hochrainer | A three-dimensional continuum theory of dislocation systems: kinematics and mean-field formulation[END_REF][START_REF] Sandfeld | Numerical implementation of a 3D continuum theory of dislocation dynamics and application to micro-bending[END_REF][START_REF] Zaiser | Some steps towards a continuum representation of 3D dislocation systems[END_REF][START_REF] Zaiser | Modelling size effects using 3D densitybased dislocation dynamics[END_REF][START_REF] Hochrainer | Continuum dislocation dynamics: towards a physical theory of crystal plasticity[END_REF][START_REF] Xia | Computational modelling of mesoscale dislocation patterning and plastic deformation of single crystals[END_REF]. Here, we suggest some important points for these developments. In an infinite single crystal, a dislocation is necessarily produced by a Frank-Read source. The left hand side of figure 4.9 illustrates this mechanism (1-1). Two defects fix locally the dislocation which will turn around when a shear stress is applied. That leads to a loop, illustrated by the symbol (2). Afterward, this loop will meet others dislocations (forest mechanism), illustrated by the right hand side of figure 4.9. At this moment, the initial perfect loop (2) tends to a loop which, locally, display ondulations. If the forest obstacle may be considered as randomly distributed, we could consider that the characteristic length scale of these ondulations is approximately time independent. In fact, the curvature of a dislocation carries two distinctive physics which should emerge differently at the mesoscopic scale. The overall curvature (or average curvature) should represent an increase of the dislocation density. And then, the small additional deviations (the local ondulations) should act only on a 3D friction term equivalent to our • s f . In other words, some relevant assumptions could be performed after splitting the curvature in two parts, an average curvature and its complementary fluctuations. The numerical treatment of a 3D dislocation density model should be quite similar to the 2D case (see chapter 3).

Here, we propose a reflexion on the second stage which would be devoted to the estimation of the driving force of the 3D model. The statistical procedure performed in 2D in chapter 2 is difficult and needs very strong assumptions such as one slip system, edge dislocation, Ÿ π fl,e t c . I nt h e general 3D case, such statistics will appear considerably more complex as well as they will require exorbitant simulation costs. Moreover, the correlation terms will have a high number of dependencies. For example, the friction term • s f in the simple case of two slip-systems of edge dislocations gives three integrals of correlation maps and each integral will have five dependencies. Consequently, the investigations will be on 15 dependencies. Such studies, based on a systematic statistical method could be unattainable (without mentioning the numerical aspect). We suggest that a direct study of the strain rate (3.7), linked to the driving force, could be a good complementary way to intuit the relevant ingredients of this driving force, more precisely, the velocity field for each slip system. First hand, this study could be performed on 3D DDD simulations including more physics of the plasticity such as junctions, cross-slip, etc. Second hand, the support of some specific experimental results, exhibiting the microscopic behavior, could be a good test to confirm the intuited driving force. For example, figure 4.10 illustrates the the tensile tests of monocrystals of copper which shows the strong dependence of the lattice orientation on the plastic strain. q kl C ijkl A kl is the sum of terms that are not known on the same grid. In order to express all components of B ij on the same grids as the other symmetric tensors, a local interpolation is required. In this method, we use the most simple interpolation based on the use of the four nearest neighbors. Because the position of these four points depends on the components ij and kl of the tensors, we introduce the interpolation operator: I (ij,kl) whose action on a given discrete field f is gathered in the following table

I (ij,kl) (f )= ij\kl (i, i)( 2 , 3) (1, 3) (1, 2) (i, i) (2, 3) (1, 3) (1, 2) Q c c c c a fL 0 11 L 10 1 L 11 0 L 011 fL 110 L 101 L 101 L 1 10 fL 0 11 L 110 L 10 1 L 01 1 f R d d d d b (C.2)
The interpolation scheme is

L s 1 s 2 0 (f )= f (-, -, ")+f (-+ s 1 , -, ")+f (-, -+ s 2 , ")+f (-+ s 1 , -+ s 2 , ") 4 . (C.3)
where (-, -, ") are the grid indexes and s i = ±1. The interpolation is similar by index rotation for L 0s 2 s 3 and L s 1 0s 3 . Using the above interpolation scheme, the polarisation term in equation (C.1) becomes

• ij = I (ij,kl) (C 0 ijkl ' ú kl ≠ "C ijkl (' kl ≠ ' ú kl )) (C.4)
where the summation over indices in subscript position is implicit.

89

This appendix proposes systematic comparisons of numerical results for different parameters in case of a square precipitate in a matrix. This precipitate is aligned to the grid of a periodic system and its width is equal to a quater of the period L. The boundary conditions are held by the average strain È'Í. We compare on each line the influence of the discretization on the maps ' xx , ' xy and ‡ xx , ‡ xy in GPa. These comparisons are performed for each method, for different elastic inhomogeneities and in case of isotropy and anisotropy. The elastic inhomogeneity is expressed by the contrast ratio r defined by C p = r • C m where C p and C m are respectively the elastic tensor of the precipitate and the matrix. There is three distinctive cases, the soft matrix with r = 10 3 ,t h eh a r dm a t r i x with r = 10 ≠3 and the homogeneous media with r =1. The two inhomogeneous cases are held at È' xx Í = È' yy Í =0.01. The homogeneous case owns a precipitate dilatation rather than an applied strain in order to obtain an heterogeneous response. It is traduced by an eigenstrain ' ú xx = ' ú yy =0.01 inside the precipitate and È'Í =0. For the isotropic case we take ⁄ m = 100 GPa and µ m = 50 GPa, and for the anisotropic case we take . The error is defined on the entire space by the difference between the simulation and a perfect reference. For the three weakest discretizations, we take the 1024 2 simulations as reference and average it on grids corresponding to the smaller discretizations, respectively 16 2 , 64 2 and 256 2 . Then, we propose three information:

C m = Q c c c c c c c c
-Min: the minimal value of this error in the entire map.

-Max: the maximal value of this error in the entire map.

-SD: the standard deviation of all values of this error. Note that each equivalent simulation performed with the four methods, we impose the same color scale. Moreover, this color scale is taken in order to illustrate both information: the global field and the artefacts. 1.78e -3 1.78e -3 1.78e -3 4.16e -1

-1.78e -3 -5.33e -3 -8.89e -3 -1.24e -2 -1.60e -2 ‡ xx 4.07e -1 2.39e -1 7.03e -2 -9.80e -2 -2.66e -1 -4.35e -1 -6.03e -1 -7.71e -1 -9.40e -1 -1.11e 0 ‡ xy 1.62e 0 1.26e 0 9.01e -1 5.41e -1 1.80e -1 -1.80e -1 -5.41e -1 -9.01e
-1.78e -3 -5.33e -3 -8.89e -3 -1.24e -2 -1.60e -2 ‡ xx 4.07e -1 2.39e -1 7.03e -2 -9.80e -2 -2.66e -1 -4.35e -1 -6.03e -1 -7.71e -1 -9.40e -1 -1.11e 0 ‡ xy 1.62e 0 1.26e 0 9.01e -1 5.41e -1 1.80e -1 -1.80e -1 -5.41e -1 -9.01e
-1.78e -3 -5.33e -3 -8.89e -3 -1.24e -2 -1.60e -2 ‡ xx 4.07e -1 2.39e -1 7.03e -2 -9.80e -2 -2.66e -1 -4.35e -1 -6.03e -1 -7.71e -1 -9.40e -1 -1.11e 0 ‡ xy 1.62e 0 1.26e 0 9.01e -1 5.41e -1 1.80e -1 -1.80e -1 -5.41e -1 -9.01e
3.44e -1

2.71e -1

1.99e -1

1.27e -1 4.16e -1

5.47e

3.44e -1

2.71e -1

1.99e -1

1.27e -1 4.16e -1

5.47e

3.44e -1

2.71e -1

1.99e -1

1.27e -1 4.16e -1

5.47e

3.44e -1

2.71e -1

1.99e -1

1.27e -1 2.07e -2

5.47e

1.72e -2

1.37e -2

1.02e -2

6.72e -3 2.07e -2

3.22e

1.72e -2

1.37e -2

1.02e -2

6.72e -3 2.07e -2

3.22e

1.72e -2

1.37e -2

1.02e -2

6.72e -3 2.07e -2

3.22e

1.72e -2

1.37e -2

1.02e -2

6.72e -3 -4.87e -4.87e Finally, we have developed an efficient numerical scheme to integrate the transport equations as well as a multigrid spectral scheme to solve elastic equilibrium associated to an arbitrary eigenstrain in an elastically heterogeneous and anisotropic medium.

3.22e

Figure I. 1 :

 1 Figure I.1: Illustration of crystal defects in a crystal: On the left, there is an interface (dim. 2 )between two phases with different orientations. On the top, there is a vacancy (dim. 0 ). On the right, there is a portion of dislocation loop (dim. 1 ) enclosing an area sheared by an elementary translation b, the Burgers vector. l is the dislocation line vector. Drawing inspired by[START_REF] Read | Dislocations in Crystals[END_REF].

Figure I. 3 Figure I. 3 :

 33 Figure I.3: Translation from the left to the right of an edge dislocation due to a shear load of the cubic lattice symbolized by two black pointers.

1 µm 5 µm

 15 Figure I.4: Dislocation-pattern of dislocations (black regions) in a monocrystal of copper at stage II of the strain hardening on the left and cycled at a strain amplitude in the plateau on the right. Revealed by typical transmission micrographs. Respectively after [64, 71].

Figure I. 5 :

 5 Figure I.5: Tensile test response of a polycrystal of aluminium: Above, The lattice orientation of each phase of the microstructure. Below, the inter-phases pattern of plasticity. The fluctuations of dislocation density are obtained by the DIC technique. After [112].

Figure 2 . 2 e 1 - 6 for b toward + e 1 for b toward -e 1 AdimensionalFigure 2 . 1 :

 22161121 Figure 2.1: Example of a pseudo-periodic simulation with imposed shear stress • =1,with50 positive dislocations in red and 50 negative dislocations in blue.

2π( 1 -Figure 2 . 2 :

 122 Figure 2.2: Equilibrium of stress contributions of the negative dislocation (equation (2.4)). The right hand side of the equation is represented by the red curve and the left hand side by the black line. The equilibrium position is represented by the crossing at x eq .

6 Figure 2 . 3 :

 623 Figure 2.3: Evolution of the stationary position of a dipole of dislocations with opposite Burgers vector function of the external stress.The right side illustrate the system. The left side compare the approximation of a dipole in our pseudo-periodic system compared to the theory. The dipole size is "y =3.9 • 10 ≠2 • L.

Figure 2 . 4 :

 24 Figure 2.4: Evolution of the strain rate function of the time for different time discretization.

Figure 2 . 5 :Figure 2 . 6 :

 2526 Figure 2.5: Stress jump due to the pseudo-periodicity. The right side shows the system for different N rep and the left side shows the value of the stress jump depending on N rep .parameter is illustrated by the right hand side of figure2.5. The plot on the left hand side shows the stress jump received by a fixed dislocation when the other one goes through the interface. We observe that the jump cannot be attenuated by adding more and more layers of boxes. This behavior comes from the fact that an additional layer involves n extra dislocations where each contribution is proportional to the inverse of the distance from the main box (see equation (2.2)). The contribution of each dislocation being proportional to 1/n,t h ep r o d u c tn • 1/n gives a jump equivalent to one dislocation positioned in the first neighbor box. A solution could be to copy more boxes toward x than toward y but we need a lot of copies to have a negligible jump and it will considerably affect the simulation cost. The best solution appears to put an artificial cut-off r co on the range of the stress induced by a dislocation • s ind (r).T h i sc u t -o ff is given in number of box size L.B e l o w this cut-off the stress is unchanged and above the stress is zero. The test illustrated by figure 2.6 still involves two dislocations in the main box with opposite Burgers vector and positioned with the largest possible spacing between each dislocation: "y =0.5L. This spacing minimizes the interactions

Figure 2 . 7 :

 27 Figure 2.7: Relaxation of a random initial distribution with 100 dislocations and without external load: Strain rate of the distribution function of the time and depending on different discretization ∆ t. " ỹmin is equal to 10 ≠2 .

Figure 2 . 8 :

 28 Figure 2.8: The averaged strain-rate of a high quantity of distributions during the relaxation for different applied stresses. The left investigation with 400 dislocations is performed by [77] with reviewed scales. The right side is a new investigation with a pseudo-periodic system describe previously involving different number of dislocations: N = 50, N = 100 or N = 200. Our averaged strain-rate is obtained with 10 3 distributions.

Figure 2 . 9 :

 29 Figure 2.9: The averaged strain-rate of 10 3 distributions during the relaxation for two high applied stresses • =1.68 and • =3.37 with N = 100.

Figure 2 . 10 :

 210 Figure 2.10: Correlation maps of size L◊L, averaged on a 256◊256 grid and for different applied stresses. Simulations are performed with N + = N ≠ = 200, r co =2 and ∆ t = 10 ≠4 . Top row: the function d +≠ (r). Bottom row: the function d ++ (r).

Figure 2 . 11 :

 211 Figure 2.11: Numerical results of the friction term • f depending on the applied stress and the period.

Fig- ure 2 .Figure 2 . 12 :

 2212 Figure 2.12: Zoom of the correlation maps d +≠ (top row) and d ≠≠ (bottom row). The adimensional size of the zoom is 5◊5. Simulations are performed for different size of main box of periodic-system. The left size is performed with Lx = 10 and Ly = 20 and the right side is performed with Lx = Ly = 20.

. 3 )

 3 In fact, • ext and • sc (r),a n dt h e r e f o r e• (r), should be considered as the shear components of the stress tensors ‡ ext , ‡ sc (r) and ‡(r)= ‡ ext + ‡ sc (r),r e s p e c t i v e l y .

. 10 )

 10 where the fitting coefficients B anddepend on the coarse-graining parameter L Ô fl.T h er e s u l t sf o r the three different values of L Ô fl that have been considered in Figure2.11 are presented in Table3.1.

Figure 3 . 1 :

 31 Figure 3.1: Symbols: Numerical estimation of • ≠ • f , coming from section 2.3.2 for three different values of L Ô fl. Lines: corresponding fit using a power law with an exponent -=1.71.

Figure 3 . 2 :

 32 Figure 3.2: Illustration of the grid along x. The subscript of fl is similar for every densities.

. 20 )Figure 3 . 3 2 ρFigure 3 . 3 :

 2033233 Figure 3.3 illustrates also the evolution of the square wave with the LW -HR scheme. There is two

Figure 3 . 4 :

 34 Figure 3.4: Movement of a square wave of dislocation density shifting on the left with a nonconstant velocity inside a 1D periodic grid with 256 nodes. Both simulations are performed with ∆t =2.6 • 10 ≠1 ∆t crit .T h es n a p sa tt 1 and at t 2 are taken respectively at the iteration 440 and 880.

Figure 3 .

 3 5 illustrates in blue this new intermediate grid in 1D.

Figure 3 . 5 :

 35 Figure 3.5: Grid along x with the localization of the flux f between 2 nodes of the main grid.

Figure 3 . 6 :

 36 Figure 3.6: Movement of a square wave of dislocation density shifting on the left with a non-constant velocity inside a 1D periodic grid with 256 nodes. The snaps at t 1 and at t 2 are taken respectively at the iteration 440 and 880 for the blue curve and at the iteration 4400 and 8800 for the green curve.

Figure 3 .κFigure 3 . 7 :

 337 Figure 3.7: Simulation of a positive dislocation attracted by a negative dislocation inside a box of 1024◊64 nodes. Both dislocations are defined respectively by a bump of density of fl + and fl ≠ .T h e side of fl + is represented here (equivalent behavior for fl ≠ ): the initial state in black, the shifted state in green. The dot fields represent Ÿ computed by the discrete Kröner equation: the red dots use a naive discretization and the blue dots use a consistent discretization.

3 •Figure 3 . 8 :

 338 Figure 3.8: Ÿ(⁄ ≠ ) (left) and Ÿ(⁄ + ) (right), vs the wavelength q of the perturbations for an applied stress • ext =5.34.

Figure 3 . 9 :

 39 Figure 3.9: Initial configuration (left) and dislocations density fields at t = 347 (right) obtained with the model including the correlation based local stresses. The top row displays fl≠fl 0 and the bottom row Ÿ.T h ea p p l i e ds t r e s si s• ext =5.34.

Figure 3 .

 3 Figure 3.10: Left: snapshot of GND density field Ÿ at t = 347 (figure 3.9). Right: Amplitude map of the corresponding Fourier transform.

Figure 3 . 11 :

 311 Figure 3.11: Initial configuration (left) and dislocations density fields at t = 347 (right) obtained with the mean field model. The top row displays fl≠fl 0 and the bottom row Ÿ.T h ea p p l i e ds t r e s si s •ext =5.34. The color scale is chosen to allow a direct comparison with Figure 3.9.

u 1 u 2 u 3 dFigure 4 . 1 :

 2341 Figure 4.1: Staggered grids for mechanical fields in 2D on the left and in 3D on the right.

Figure 4 . 2 :

 42 Figure 4.2: One voxel containing four sub-voxels with homogeneous elastic energy with the illustrated eigenstrain and stiffness tensors decomposition.

  l q(3) m e iq• Pil e ≠iq• Pmn ' J (il,mn) C 0 ilmn(4.26) where' J (il,mn) = 1 2 Dÿ s e iq•( ◊il (s)≠ ◊mn(s)) .(4.27)

Figure 4 .

 4 Figure 4.3 shows a comparison between the three methods for the highest discretization (1024run) in order to confine local artefact and stand out, far from corners, non-local artefacts. On the right side, there is a profil of ' 11 along the black line illustrated by the map on the left side.

Figure 4 . 3 :

 43 Figure 4.3: On the left: Map of ' 11 (Method (1)) with the localisation of the profil (black line) and the localisation of figure 4.4 (black square). On the right: Profil of ' 11 through the interface for each method.

Figure 4 . 4 :

 44 Figure 4.4: Difference between ' (64) 11 (64-run) and ' (1024) 11 (1024-run averaged on a 64 ◊ 64 grid) for each method. Only one corner of the precipitate is exhibited, localized by the square in figure 4.3.

Figure 4 . 5 :

 45 Figure 4.5: Fields of 2fi(1 ≠ ‹) ‡ yy /µ computed in a box with 64◊64 nodes. The three maps are, in the order, obtained by method (1), (2) and (3). The two dislocations are positioned at[START_REF] Daveau | Interaction dislocations -joints de grains en déformation plastique monotone : étude expérimentale et modélisations numériques[END_REF][START_REF] Fressengeas | An elasto-plastic theory of dislocation and disclination fields[END_REF] and[START_REF] Geers | Coupled glide-climb diffusion-enhanced crystal plasticity[END_REF][START_REF] Fressengeas | An elasto-plastic theory of dislocation and disclination fields[END_REF]. The profil, coming from the method (3) and illustrated by the dotted line on the third map, is compared to the analytical solution.

Figure 4 .Figure 4 . 6 :

 446 Figure 4.6: Maps of ' 11 (64-run) for each method, containing a quarter of the periodic system with the middle of the precipitate matching with the corner of the figure.

(4. 29 )

 29 As required, we hold new contributions C m 1112 = 19 GPa. We still keep the same contrast (C p ijkl /C m ijkl = 10 ≠3 )a n dt h es a m ed r i v i n gf o r c eÈ' 11 Í = È' 22 Í =0.01. Method (1) and (2) keep the same artefact, already illustrated with the 'aligned case'. As for method (3) and (4), figure4.7 shows the resulting field of ' 11 . Appendix F proposes also this test with the opposite contrast for all methods.

Figure 4 . 7 :

 47 Figure 4.7: Maps of ' 11 , containing a quarter of the periodic system with the middle of the precipitate matching with the corner of the figure. Above: 64-run. Below: 1024-run. Left: method (3). Right: method (4).

  Figure 4.8 illustrates a titanium alloy response of a one-directional load along the vertical axis È' 22 Í =0.01 where we consider only --phases assumed as a Voronoï microstructure. The elastic stiffness is given by the following tensor expressed in the natural basis of the grid (Mandel notation in GPa): of figure 4.8 shows the angle representing the phase orientation in the natural basis.

Figure 4 . 8 :

 48 Figure 4.8: Example of a Voronoï microstructure of titanium alloy with a crack. On the left: the map of the angle (in gradian) of the cubic symmetry. On the right: the map of Von Mises stress in GPa (method (4)) for a polycrystalline and a monocrystalline.

1 1 2 Figure 4 . 9 :

 249 Figure 4.9: On the left: Illustration of a perfect Frank-Read source inside a shear load simulation, after[START_REF] Fivel | Frank-Read multiplication mechanism simulations[END_REF]. On the right: Illustration inside the plane of a curved dislocation of its movement through a forest of dislocations, after[START_REF] Foreman | The bowing of a dislocation segment[END_REF].

Figure 4 .

 4 Figure 4.10: Load-elongation curves of copper single crystals with various multiple-glide orientations at 22 ¶ C. After [95].

99

 99

F- 4 - 8 .67e - 4 - 2 .54e - 3 - 4 .22e - 3 - 5

 48423435

- 4 - 8 .67e - 4 - 2 .54e - 3 - 4 .22e - 3 - 5

 48423435

- 4 - 8 .67e - 4 - 2 .54e - 3 - 4 .22e - 3 - 5

 48423435

- 4 - 8 .67e - 4 - 2 .54e - 3 - 4 .22e - 3 - 5 - 3 - 1 .78e - 3 - 5 .33e - 3 - 8 .89e - 3 - 1 .24e - 2 - 1 . 1 2.39e - 1 7.03e - 2 - 9 .80e - 2 - 2 .66e - 1 - 4 .35e - 1 - 6 .03e - 1 - 7 .71e - 1 - 9 .40e - 1 - 1 .11e 0 ‡ xy 1 - 1 - 1 .80e - 1 - 5 .41e - 1 - 9

 48423435313538312111292214161719111111519 60e -2 ‡ xx 4.07e -.01e

- 2 - 1 .76e - 2 - 8 .98e - 2 - 1 - 2 - 8 .53e - 2 - 2 .56e - 1 - 4 .27e - 1 - 5 .97e - 1 - 7

 2128212822141517

- 2 - 1 .76e - 2 - 8 .98e - 2 - 1 - 2 - 8 .53e - 2 - 2 .56e - 1 - 4 .27e - 1 - 5 .97e - 1 - 7 - 2 - 8 .53e - 2 - 2 .56e - 1 - 4 .27e - 1 - 5 .97e - 1 - 7 .

 21282128221415172822141517 68e-1

- 2 - 1 .76e - 2 - 8 .98e - 2 - 1 - 2 - 8 .53e - 2 - 2 .56e - 1 - 4 .27e - 1 - 5 .97e - 1 - 7

 2128212822141517

- 2 - 1 .76e - 2 - 8 .98e - 2 - 1 - 2 - 8 .53e - 2 - 2 .56e - 1 - 4 .27e - 1 - 5 .97e - 1 - 7

 2128212822141517

- 3 - 2 - 3 - 4 .16e - 3 - 1 .25e - 2 - 2 .08e - 2 - 2 .91e - 2 - 3

 323431222223

- 3 - 2 - 3 - 4 .16e - 3 - 1 .25e - 2 - 2 .08e - 2 - 2 .91e - 2 - 3

 323431222223

- 3 - 2 - 3 - 4 .16e - 3 - 1 .25e - 2 - 2 .08e - 2 - 2 .91e - 2 - 3

 323431222223

- 3 - 2 - 3 - 4 .16e - 3 - 1 .25e - 2 - 2 .08e - 2 - 2 .91e - 2 - 3

 323431222223

9.37e - 3 7.13e - 3 4.90e - 3 2.67e - 3 4.33e - 4 - 1 .80e - 3 - 4 .03e - 3 - 6 .27e - 3 - 8 - 3 - 1 .53e - 3 - 4 .60e - 3 - 7 .67e - 3 - 1 .07e - 2 - 1 . 1 3.00e - 1 1.57e - 1 1.40e - 2 - 1 .29e - 1 - 2 .72e - 1 - 4 .15e - 1 - 5 .58e - 1 - 7 .01e - 1 - 8

 3333413436383134373121111211214151718 38e -2 ‡ xx 4.43e -

9.37e - 3 7.13e - 3 4.90e - 3 2.67e - 3 4.33e - 4 - 1 .80e - 3 - 4 .03e - 3 - 6 .27e - 3 - 8 - 3 - 1 .53e - 3 - 4 .60e - 3 - 7 .67e - 3 - 1 .07e - 2 - 1 . 1 3.00e - 1 1.57e - 1 1.40e - 2 - 1 .29e - 1 - 2 .72e - 1 - 4 .15e - 1 - 5 .58e - 1 - 7 .01e - 1 - 8

 3333413436383134373121111211214151718 38e -2 ‡ xx 4.43e -

9.37e - 3 7.13e - 3 4.90e - 3 2.67e - 3 4.33e - 4 - 1 .80e - 3 - 4 .03e - 3 - 6 .27e - 3 - 8 - 3 - 1 .53e - 3 - 4 .60e - 3 - 7 .67e - 3 - 1 .07e - 2 - 1 . 1 3.00e - 1 1.57e - 1 1.40e - 2 - 1 .29e - 1 - 2 .72e - 1 - 4 .15e - 1 - 5 .58e - 1 - 7 .01e - 1 - 8

 3333413436383134373121111211214151718 38e -2 ‡ xx 4.43e -

9.37e - 3 7.13e - 3 4.90e - 3 2.67e - 3 4.33e - 4 - 1 .80e - 3 - 4 .03e - 3 - 6 .27e - 3 - 8 - 3 - 1 .53e - 3 - 4 .60e - 3 - 7 .67e - 3 - 1 .07e - 2 - 1 . 1 3.00e - 1 1.57e - 1 1.40e - 2 - 1 .29e - 1 - 2 .72e - 1 - 4 .15e - 1 - 5 .58e - 1 - 7 .01e - 1 - 8

 3333413436383134373121111211214151718 38e -2 ‡ xx 4.43e -

7.63e - 1 4.87e - 1 2.11e - 1 - 6 .43e - 2 - 3 .40e - 1 - 6 .16e - 1 - 8 .92e - 1 - 1 - 1 - 1 .66e - 1 - 4 .98e - 1 - 8 .30e - 1 - 1

 11162316181111141811

7.63e - 1 4.87e - 1 2.11e - 1 - 6 .43e - 2 - 3 .40e - 1 - 6 .16e - 1 - 8 .92e - 1 - 1 - 1 - 1 .66e - 1 - 4 .98e - 1 - 8 .30e - 1 - 1

 11162316181111141811

7.63e - 1 4.87e - 1 2.11e - 1 - 6 .43e - 2 - 3 .40e - 1 - 6 .16e - 1 - 8 .92e - 1 - 1 - 1 - 1 .66e - 1 - 4 .98e - 1 - 8 .30e - 1 - 1

 11162316181111141811

7.63e - 1 4.87e - 1 2.11e - 1 - 6 .43e - 2 - 3 .40e - 1 - 6 .16e - 1 - 8 .92e - 1 - 1 - 1 - 1 .66e - 1 - 4 .98e - 1 - 8 .30e - 1 - 1

 11162316181111141811

- 3 - 3 .22e - 3 - 9 .67e - 3 - 1 .61e - 2 - 2 .26e - 2 - 2

 3339312222

- 3 - 3 .22e - 3 - 9 .67e - 3 - 1 .61e - 2 - 2 .26e - 2 - 2

 3339312222

- 3 - 3 .22e - 3 - 9 .67e - 3 - 1 .61e - 2 - 2 .26e - 2 - 2 - 3 - 3 .22e - 3 - 9 .67e - 3 - 1 .61e - 2 - 2 .26e - 2 - 2 2 Résumé

 333931222233393122222 Le comportement mécanique des alliages métalliques industriels, notamment ceux utilisés dans le domaine de l'aéronautique, est contrôlé par la présence de différents types de précipités et par la nucléation et propagation de défauts cristallins tels que les dislocations. La compréhension du comportement de ces matériaux nécessite des modèles continus afin d'accéder à l'échelle macroscopique. Cependant, même aujourd'hui, les théories conventionnelles de la plasticité utilisent des variables mésoscopique et des équations d'évolution qui ne reposent pas sur la notion de transport de dislocations. En conséquence, ces théories sont basées sur des lois phénoménologiques qu'il est nécessaire de calibrer pour chaque matériau et chaque application. Il est donc souhaitable d'établir le lien entre les échelles micro et macro afin de générer une théorie continue de la plasticité déduite analytiquement des équations fondamentales de la dynamique des dislocations. L'objet de cette thèse est précisément de contribuer à l'élaboration d'une telle théorie. La première étape a consisté à établir rigoureusement la procédure de changement d'échelle dans une situation simplifiée. Nous avons alors abouti à un système d'équations de transport hyperboliques sur des densités de dislocations contrôlées par des contraintes locales de friction et de backstress qui émergent du changement d'échelle. Nous avons ensuite développé une procédure numérique pour calculer ces termes et analyser leur comportement. Finalement, nous avons développé un schéma numérique efficace pour intégrer les équations de transport ainsi qu'un schéma spectral multi-grille pour résoudre l'équilibre élastique associé à un champ de déformation propre quelconque dans un milieu élastiquement anisotrope of industrial metallic alloys, in particular those used in the aerospace industry, is controlled by the existence of several types of precipitates and by the nucleation and propagation of crystalline defects such as dislocations. The understanding of this behavior requires continuous models to access the use mesoscopic variables and evolution equations that are not based on the transport of dislocations. Therefore, these theories are based on phenomenological laws that must be calibrated for each material, or, for each specific application. It is therefore highly desirable to make link between the micro and macro scales, in order to derive a continuous theory of plasticity from the fundamental equations of the dislocation dynamics. The aim of this thesis is precisely to contribute the elaboration of such a theory. The first step has consisted to rigorously establish a coarse graining procedure in a simplified situation. We have then obtained a set of hyperbolic transport equations on dislocation densities, controlled by a local friction stress and a local back-stress that emerge from the scale change. We have then developed a numerical procedure to compute these local terms and analyze their behavior.

.5 Physical meaning of the local stresses τ s f and τ s b

	where the local stresses • sc (r, t), • s f (r, t) and • s b (r, t) are defined in equations (1.23), (1.31) and
	(1.32). Next, we discuss the physical meaning of the correlation-induced • s f and • s b .
	1.2
	1. We will precise this point in chapter 2, section 2.3.1

Table 3 .

 3 1: Coefficients of the power law fit.

	L Ô fl	10	14.1	20
	-	1.655 1.719 1.765
	B	0.156 0.197 0.246

Table 3 .

 3 71, B =0.16, C ++ =1.09 and C +≠ = ≠0.23 for which correlations are accounted for. The other parameters are reported in Table (3.2). 2: Shear modulus µ, Poisson coefficient ‹ and Burgers vector b

  a cubic symmetry turned by 18 ¶ with respect to our grid. We perform each simulation with 4 different discretizations ordered in column, respectively: 16◊16, 64◊64, 256◊256 and 1024◊1024. Moreover, we precise the characteristic values of the error due to this discretization for ' xx above the columns: 16 2 , 64 2 and 256 2

	c a	133.8 76.2 90.0 0 0	76.2 90.0 133.8 90.00 0 90.0 120.00 0 0 110.00 0 0≠26.9 26.9 0 0 0 0 0 0 110.00	R d d d d d d b d d d	GPa	(F.1)
		26.9 ≠26.9	0	0	0	82.4
	which corresponds to

Recently, M. Geers and coll.[START_REF] Dogge | Extended modelling of dislocation transport-fomulation and finite element implementation[END_REF] have proposed a set of different transport equations based on dislocation densities. However, the derivation does not rely on a coarse-graining procedure and, thus, does not correspond to a real transition to mesoscale.

This is coherent with the results presented in[START_REF] Ispánovity | Submicron Plasticity: Yield Stress, Dislocation Avalanches, and Velocity Distribution[END_REF], where it is observed that, for a similar system, the plastic strain rate as a function of the applied stress is always finite, even-though it does display two different regimes separated by a critical stress.

where L s 1 s 2 0 (f )= f (-, -, ")+f (-+ s 1 , -, ")+f (-, -+ s 2 , ")+f (-+ s 1 , -+ s 2 , ") 4 (4. [START_REF] Cottrell | Strain hardening in Andrade creep[END_REF] with (-, -, ") the grid index and s i = ±1. A similar definition stands for L 0s 2 s 3 and L s 1 0s 3 .U s i n g this operator, the polarization tensor is computed in real space as

and the elastic equilibrium becomes

(4.17

where we have assumed an implicit summation over repeated indices in subscript position. As usual, equation (4.17) is solved for ' ij in Fourier space.

As detailed in appendix C, we get an equation similar to equation (4.4) where the strain Green tensor is now defined as

ik q

(3) j + G

(3) jk q

(3) i 4

• q

(3)

il q

(3)

where the displacement Green tensor G (3) is given by

(q)=e iq• Pij e ≠iq• Pkl ' I (ij,kl) C 0 ijkl q

(3) j q

(3) k (4.19) where Pij is equal to ≠ d 2 (1≠" i,j )(ę i +ę j ) as detailled in appendix B. ' I (ij,kl) is a complex number defined in table C.10 of appendix C. Note finally that in many cases of practical importance, the elastic stiffness tensor only differs from zero when i = j and k = l or when (i, j)=(k, l). This is for example the case for the isotropic case, but also for several anisotropic symmetries (cubic, tetragonal, orthorhombic) provided that the numerical grid is suitably oriented with the crystal lattice.

Method (4): Energy based interpolation

This last method, based on a variational principle, is computed with a direct solving of the mechanical equilibrium. But the fixed point algorithm can also be implemented, following the appendix D.Firstly, we will gives the system of linear equations and secondly the fixed point algorithm.

This method define a physical interpolation with consistency with the definition of the elastic energy. Moreover, in order to express an accurate discrete energy, we split the main grid (circle-grid) by small voxels (d/2 ◊ d/2), as explained by figure 4.2. As mentioned before, it is assumed that materials parameters (C ijkl and ' ú ij ) are known on the circle-grid and interpolated on cross-grid with linear interpolation of the four neighbors.

Appendix A

Fourier Transform

Let f (r) be a function defined in an orthorhombic box with periodic boundary conditions. The Fourier transform is defined as

where V is the volume of the box. Due to the choice of periodic boundary conditions, the wave vectors q are discrete:

wherei are signed integers and where L i are the box dimensions.

The backward Fourier Transform is

If we are interested in continuous functions whose spatial variations are on a distance larger than the characteristic length d, the Fourier space can be limited to the first Brillouin zone ' B 0 defined as the set of vectors q wherei are integers such that ≠N i /2 <i AE N i /2 where N i = L i /d. Then f (r) is fully characterized by these N 1 N 2 N 3 Fourier components. Alternatively f ( (r) can be characterized by its value on a regular grid B 0 defined as the set of vectors r =(u 1 ,u 2 ,u 3 )d + r0 where r0 is an arbitrary shift vector and u i are integers such that 0 AE u i <N i . Introducing the discrete field f d (u 1 ,u 2 ,u 3 )=f ((u 1 ,u 2 ,u 3 )d+ r0 ) and its discrete Fourier Transform

where

This expression is only the usual property that a global shift in real space is equivalent to a multiplication by a phase term in Fourier space.

Appendix B Deformation and displacements with staggered grids

In the simulations with periodic boundary conditions, the total deformation tensor '( R) is decomposed in a uniform macroscopic strain È' ij Í and an heterogeneous part "'( R), whose spatial average is zero. We then introduce local displacement ų( R) which are related to the heterogeneous strain by the usual relation:

Using the staggered grids, the displacement and deformation fields are defined on different grids as shown in figure 4.1. It leads to the following relationships:

where (ę 1 , ę2 , ę3 ) are the unit vectors of the cubic grid, R is a vector of the grid and d is the grid size. Expressions (B.2) can be rewritten in the following convenient form:

where

and " ij stands for the Kronecker symbol equal to 1 if i = j and 0 otherwise. Going to Fourier space, equation (B.3) becomes

where q

(3)

Appendix C Construction of the fixed point scheme of Method (3) for mechanical equilibrium

This appendix builds the terms used for method (3). We start from equation (4.13) which is rewriten as follow:

We precise the grid position of each component of the divergence by a subscript below the central difference. R describes the main discrete grid (circle-grid) and ( ę1 , ę2 , ę3 ) are the unit vectors of our basis as detailed in figure 4

As usual, this equation is solved in Fourier space. The righthand side of equation (C.5) is

where we have used the notations

In Fourier space we obtain

Using the Fourier vectors q(3) and the vector Pij defined in B, we get the close form

Similarly, the Fourier transform of the lefthand side of equation (C.5) is

In Fourier space, the interpolation operator I (ij,kl) is only a multiplication per a complex number ' I (ij,kl) (q) defined by the following table

with s i = ±1. Similar definitions hold for ' L s 1 0s 3 and ' L 0s 2 s 3 . The mechanical equilibrium is therefore

Using equation (B.4) linking the strain to the displacements in equation (C.12) and using the symmetries of C 0 ijkl and ' I (ij,kl) we obtain, for all l,

where the Green tensor is defined for all q (3) " =0by

which is Hermitian.

Introducing the displacement (C.13) in equation (B.4) we get

where the strain Green tensor is

ik q

(3)

il q

(3)

In the above expression, we have enforced the k ¡ l symmetry Γ ijkl = Γ ijlk = Γ jikl in order to treat Γ in Voigt notation.

Appendix D Construction of Method (4) with variational approach for mechanical equilibrium

In this method, we assume that the strain component fields are defined on staggered grids (see figure 4.1) and are related to the local displacement field ų as detailed in B. We assume that all physical quantities are locally constant over a volume equal to a voxel of size d (see figure 4.2). However, due to the shifts between the staggered grids, it only implies that all quantities are constant inside subvoxels of size d/2. Therefore, each voxel centered on a lattice point R of the main ('circle') grid is divided into 2 D subvoxels where D stands for the space dimension. To define a subvoxel, we must supplement the grid vector R with an additional quantiy s to select one of the 2 D subvoxels. A convenient choice is to choose s such as R + s is the center of the subvoxel. All coordinates of s are equal to ± d 4 . The elastic energy can then be written

where C ijkl ( R, s) and ' el ij ( R, s) are respectively the elastic tensor and the elastic strain of the ( R, s) subvoxel. The sign tilde means the averaged field inside the subvoxel. Assuming that the homogeneous strain È'Í is imposed, elastic equilibrium is reached when the derivative of F el with respect to the local displacements ų is zero. 1 Computing these derivatives, we have

where ‡kl ( RÕ , s)= q mn Cklmn ( RÕ , s)' el mn ( RÕ , s) is the local stress in the subvoxel ( RÕ , s). The strain in a subvoxel 'el mn ( RÕ , s) is related to the strain ' el mn ( RÕ ) defined on the staggered grid (figure 4.1) by

where ◊kl (s) is the shift vector defined by

1. The extension to a condition of imposed average stress is straightforward. See e.g. [START_REF] Khachaturyan | Theory of Structural Transformations in Solids[END_REF].

Using equation (B.3), we get

+ " l,i

3

This expression is then introduced in equation (D.2) leading to the following expression

where we have introduced the quantity

is only the average stress over the voxel at point R. ‡ M ij ( R) with i " = j is the average stress over all subvoxels around the point R + d 2 (ę i + ęj ). Note that in equation (D.6), a divergence of this average stress tensor can be recognized. As in methods (1), ( 2) and ( 3), the stress ‡ M il ( R) is then decomposed into an elastically homogeneous term and a polarization

where the • M il ( R) is related to the polarization of the subvoxels •il ( R, s) by

and, as in equation (4.5),

Expressing equation (D.6) in Fourier space and using equation (D.8) we get

where the displacement Green tensor is

ÿ s e iq•( ◊il (s)≠ ◊mn(s)) .

(D.13)

The equilibrium displacement field is then analytically obtained in Fourier space (q " =0)

Finally, inserting equation (D.14) in equation (B.4) we obtain

where

Appendix E

Second derivatives of elastic energy for direct solving of Method (4)

The second derivatives of the elastic energy (D.1) with respect to the displacements are computed in this appendix. We introduce the dynamical matrix where ‡ M ( R) is defined by equation (D.7). The dynamical matrix is obtained by derivating this expression with respect to u j ( RÕ ). First, the derivative of the stress is:

where s i is the component of s along ęi .